Ingredients#

Ingredients are functions that encode the physics of a model. QC Lab is designed to be operational with a minimal set of ingredients that describe the Hamiltonian of the system which is assumed to have the form:

\[H(q,p) = \hat{H}_{\mathrm{q}} + \hat{H}_{\mathrm{q-c}}(q) + H_{\mathrm{c}}(q,p)\]

where \(\hat{H}_\mathrm{q}\) is the quantum Hamiltonian, \(\hat{H}_{\mathrm{q-c}}(q)\) is the quantum-classical coupling Hamiltonian, and \(H_{\mathrm{c}}(q,p)\) is the classical Hamiltonian.

A generic ingredient has the form:

def ingredient_specifier(model, parameters, **kwargs):
    # Get any keyword arguments, with default values if not provided.
    kwarg1 = kwargs.get('kwarg1', default_value1)
    # Compute the ingredient using attributes of the model and parameters objects.
    constants = model.constants
    # Return the computed ingredient.
    return ingredient

where model is a Model object (see Models) which contains all the constants of the model, parameters is a variable object (see Variable Objects) containing time-dependent parameters of the simulation, and **kwargs are any additional keyword arguments that are specific to that ingredient type.

Ingredients in QC lab can come in different variations, for example the quantum Hamiltonian ingredient could describe a two-level system, a nearest-neighbor lattice, or a more complex Hamiltonian. The type and variety of an ingredient is specified in its name which follows the convention <ingredient_type>_<variety>. For example, the quantum Hamiltonian ingredient for a two-level system is named h_q_two_level where h_q indicates that it is a quantum Hamiltonian ingredient and two_level indicates that it describes a two-level system. Likewise the classical Hamiltonian ingredient for a harmonic oscillator is named h_c_harmonic where h_c indicates that it is a classical Hamiltonian ingredient and harmonic indicates that it describes a harmonic oscillator.

Ingredients can be included in a model by appending them to the model’s ingredients attribute, which is a list of tuples where each tuple contains the name of the ingredient and the ingredient function itself. Because the ingredients list is read back to front, appending an ingredient is sufficient to overwrite any existing ingredient with the same name. For example, to include a custom quantum Hamiltonian ingredient in a model, one would do:

def h_q_custom(model, parameters, **kwargs):
    # Custom quantum Hamiltonian implementation.
    h_q = ...
    return h_q

model.ingredients.append(("h_q", h_q_custom))

The minimal set of ingredients required to run a simulation are:

A quantum Hamiltonian ingredient, named h_q.
A classical Hamiltonian ingredient, named h_c.
A quantum-classical coupling Hamiltonian ingredient, named h_qc.

Additional ingredients that make the simulation more efficient or accurate are:

An initialization function for the classical coordinates, named init_classical.
A gradient of the classical Hamiltonian with respect to the conjugate classical coordinates, named dh_c_dzc.
A gradient of the quantum-classical coupling Hamiltonian with respect to the conjugate classical coordinates, named dh_qc_dzc.
A hopping function for surface hopping algorithms, named hop.

Vectorization#

By default, QC Lab assumes that ingredients are implemented in a vectorized fashion. This means that rather than constructing the respective term for an individual trajectory, each ingredient constructs the term for all trajectories in a batch at once. Additionally, any keyword argument associated with an ingredient will have an initial trajectory dimension. For example, the quantum-classical Hamiltonian ingredient h_qc has a keyword argument z which is the complex classical coordinate. If the batch size is 100, then z will have shape (100, model.constants.num_classical_coordinates) where model.constants.num_classical_coordinates is the number of classical coordinates in the model. The output of the h_qc ingredient will then have shape (100, model.constants.num_quantum_states, model.constants.num_quantum_states) where model.constants.num_quantum_states is the number of quantum states in the model.

Rather than implementing this vectorization yourself, you can use the @vectorize_ingredient decorator provided in the qclab.functions module. This decorator will automatically vectorize an ingredient that is implemented for a single trajectory (i.e., without any batch dimension). For example, the following implementation of the quantum Hamiltonian ingredient for a two-level system is not vectorized:

from qclab.functions import vectorize_ingredient

@vectorize_ingredient
def h_q_two_level(model, parameters, **kwargs):
    """
    Quantum Hamiltonian for a two-level system.
    """
    e_0 = model.constants.two_level_00
    e_1 = model.constants.two_level_11
    v = model.constants.two_level_01_re + 1j * model.constants.two_level_01_im
    h_q = np.zeros((2, 2), dtype=complex)
    h_q[0, 0] = e_0
    h_q[1, 1] = e_1
    h_q[0, 1] = v
    h_q[1, 0] = np.conj(v)
    return h_q

The ingredient can then be included in a model as:

model.ingredients.append(("h_q", h_q_two_level))

Alternatively, vectorization can be implemented manually. For example, the following implementation of the quantum Hamiltonian ingredient for a two-level system is manually vectorized:

def h_q_two_level(model, parameters, **kwargs):
    """
    Quantum Hamiltonian for a two-level system.
    """
    batch_size = kwargs["batch_size"]
    e_0 = model.constants.two_level_00
    e_1 = model.constants.two_level_11
    v = model.constants.two_level_01_re + 1j * model.constants.two_level_01_im
    h_q = np.zeros((2, 2), dtype=complex)
    h_q[0, 0] = e_0
    h_q[1, 1] = e_1
    h_q[0, 1] = v
    h_q[1, 0] = np.conj(v)
    # Here we use np.broadcast_to to add the batch dimension to the output.
    # This is only possible because the output does not depend on any trajectory-specific
    # information.
    return np.broadcast_to(h_q, (batch_size, 2, 2))

Sparse Quantum-Classical Gradients#

If left unspecified, the gradient of the quantum-classical Hamiltonian will be calculated numerically using finite differences. This can be computationally expensive and introduce inaccuracies. This can be avoided by providing an analytical implementation of the gradient ingredient dh_qc_dzc. In QC Lab, we implement this gradient in a sparse manner, meaning that we only compute the non-zero elements of the gradient matrix. This reduces the potentially large gradient tensor with shape (sim.settings.batch_size, model.constants.num_classical_coordinates, model.constants.num_quantum_states, model.constants.num_quantum_states) to a list of non-zero elements, mels, their indices inds, and the shape of the full tensor shape.

In practice, the ordering of the matrix elements can have a dramatic impact on performance due to memory access patterns. Therefore, we recommend using numpy.where to determine the indices of the non-zero elements. For example, the following implementation of the gradient of the quantum-classical Hamiltonian for a spin-boson model computes only the non-zero elements of the gradient:

@vectorize_ingredient
def dh_qc_dzc_spinboson(model, parameters, **kwargs):
    z = kwargs["z"]
    A = model.constants.get("A")
    diag_coupling = model.constants.diagonal_linear_coupling
    out = np.zeros((A, 2, 2), dtype=complex)
    for i in range(A):
        out[i, 0, 0] = diag_coupling[0, i]
        out[i, 1, 1] = diag_coupling[1, i]
    inds = np.where(out != 0)
    mels = out[inds]
    shape = (len(z), A, 2, 2)
    return inds, mels, shape

A dense ingredient can be automatically made sparse by invoking the @make_ingredient_sparse decorator from the qclab.functions module. This decorator will convert a dense ingredient into a sparse one by determining the non-zero elements, their indices, and the shape of the full tensor. For example, the following implementation of the gradient of the quantum-classical Hamiltonian for a spin-boson model is dense, but made sparse by the decorator:

from qclab.functions import make_ingredient_sparse

@make_ingredient_sparse
@vectorize_ingredient
def dh_qc_dzc_spinboson(model, parameters, **kwargs):
    z = kwargs["z"]
    A = model.constants.get("A")
    diag_coupling = model.constants.diagonal_linear_coupling
    out = np.zeros((A, 2, 2), dtype=complex)
    for i in range(A):
        out[i, 0, 0] = diag_coupling[0, i]
        out[i, 1, 1] = diag_coupling[1, i]
    return out

Of course the most efficient implementation is one that is both analytical and sparse without invoking the decorator. This is what is implemented in the ingredient ingredients.dh_qc_dzc_diagonal_linear which is included in the Spin-Boson model by default.

Ingredients in QC Lab#

The built-in ingredients in QC Lab can be found in the qclab.ingredients module.

Note

All ingredients assume that the model object has a minimal set of constants including num_quantum_states (the number of quantum states) and num_classical_coordinates (the number of classical coordinates), classical_coordinate_mass (the mass of the classical coordinates), and classical_coordinate_weight (the weight of the classical coordinates). These constants are discussed in Models. For brevity we exclude explicit mention of these constants in the task documentation.

This module contains ingredients for use in Model classes.

qclab.ingredients.h_c_harmonic(model, parameters, **kwargs)#

Harmonic oscillator classical Hamiltonian function.

\(H_{\mathrm{c}} = \frac{1}{2}\sum_{n} \left( \frac{p_n^2}{m_n} + m_n \omega_n^2 q_n^2 \right)\)

Keyword Args

zndarray: Complex classical coordinate.

Model Constants

harmonic_frequencyndarray: Harmonic frequency of each classical coordinate.

Returns

h_cndarray: Classical Hamiltonian. (batch_size,).

qclab.ingredients.h_c_free(model, parameters, **kwargs)#

Free particle classical Hamiltonian function.

\(H_{\mathrm{c}} = \sum_{n} \left( \frac{p_n^2}{2m_n} \right)\)

Keyword Args

zndarray: Complex classical coordinate.

Model Constants

None

Returns

h_cndarray: Classical Hamiltonian. (batch_size,).

qclab.ingredients.dh_c_dzc_harmonic(model, parameters, **kwargs)#

Derivative of the harmonic oscillator classical Hamiltonian function with respect to the conjugate z coordinate. This is an ingredient that calls the low-level function dh_c_dzc_harmonic_jit.

Keyword Args

zndarray: Complex classical coordinate.

Model Constants

harmonic_frequencyndarray: Harmonic frequency of each classical coordinate.

Returns

dh_c_dzcndarray: Gradient of the classical Hamiltonian with respect to the conjugate classical coordinate. (batch_size, num_classical_coordinates).

qclab.ingredients.dh_c_dzc_free(model, parameters, **kwargs)#

Derivative of the free particle classical Hamiltonian function with respect to the conjugate z coordinate.

Keyword Args

zndarray: Complex classical coordinate.

Model Constants

None

Returns

dh_c_dzcndarray: Gradient of the classical Hamiltonian with respect to the conjugate classical coordinate. (batch_size, num_classical_coordinates).

qclab.ingredients.h_q_two_level(model, parameters, **kwargs)#

Quantum Hamiltonian for a two-level system.

\(H_{nm} = \delta_{nm}\mathrm{two\_level\_nn}+(1-\delta_{nm})(\mathrm{two\_level\_nm\_re} + i \mathrm{two\_level\_nm\_im})\)

Keyword Args

batch_sizeint: Number of trajectories in a batch.

Model Constants

two_level_00float: Energy of the first level.
two_level_11float: Energy of the second level.
two_level_01_refloat: Real part of the off-diagonal coupling.
two_level_01_imfloat: Imaginary part of the off-diagonal coupling.

Returns

h_qndarray: Quantum Hamiltonian. (batch_size, num_states, num_states).

qclab.ingredients.h_q_nearest_neighbor(model, parameters, **kwargs)#

Quantum Hamiltonian for a nearest-neighbor lattice.

\(H_{nm} = -t (\delta_{n,m+1} + \delta_{n,m-1})\)

Keyword Args

batch_sizeint: Number of trajectories in a batch.

Model Constants

nearest_neighbor_hopping_energyfloat: Hopping energy between sites \(t\).
nearest_neighbor_periodicbool: Whether to apply periodic boundary conditions.

Returns

h_qndarray: Quantum Hamiltonian. (batch_size, num_states, num_states).

qclab.ingredients.h_qc_diagonal_linear(model, parameters, **kwargs)#

Diagonal linear quantum-classical Hamiltonian.

\(H_{nm} = \delta_{nm}\sum_{j} \gamma_{nj} (z_{j} + z_{j}^*)\)

Keyword Args

zndarray: Complex classical coordinate.

Model Constants

diagonal_linear_couplingndarray: Coupling constants \(\gamma\).

Returns

h_qcndarray: Quantum-classical coupling Hamiltonian. (batch_size, num_states, num_states).

qclab.ingredients.dh_qc_dzc_diagonal_linear(model, parameters, **kwargs)#

Gradient of the diagonal linear quantum-classical coupling Hamiltonian in sparse format.

\([\partial_{z} H_{qc}]_{ijkl} = \delta_{kl}\gamma_{kj}\)

Keyword Args

zndarray: Complex classical coordinate.

Model Constants

diagonal_linear_couplingndarray: Coupling constants \(\gamma\).

Returns

indstuple of ndarray: Indices of the non-zero elements of the gradient. (batch_index, coordinate_index, row_index, column_index).
melsndarray: Values of the non-zero elements of the gradient.
shapetuple: Shape of the full gradient array. (batch_size, num_classical_coordinates, num_states, num_states).

qclab.ingredients.hop_harmonic(model, parameters, **kwargs)#

FSSH hop function for taking the classical coordinates to represent harmonic oscillators.

Determines the shift in the classical coordinates required to conserve energy following a hop between quantum states.

If enough energy is available, the function returns the shift in the classical coordinates such that the new classical coordinate is z + shift and a Boolean equaling True if the hop has occurred. If not enough energy is available, the shift becomes zero and the Boolean is False.

Solves the equation:

\[H_{\mathrm{c}}(z) + \epsilon_{\mathrm{initial}} = H_{\mathrm{c}}(z + \mathrm{shift}) + \epsilon_{\mathrm{final}}\]

Keyword Args

zndarray: Current classical coordinate.
delta_zndarray: Rescaling direction of z.
eigval_difffloat: Energy difference between final and initial states.

Model Constants

harmonic_frequencyndarray: Harmonic frequency of each classical coordinate.

Returns

shiftndarray: Shift in the classical coordinate.
hopbool: Whether the hop has occurred.

qclab.ingredients.hop_free(model, parameters, **kwargs)#

FSSH hop function taking the classical coordinates to represent free particles.

Determines the shift in the classical coordinates required to conserve energy following a hop between quantum states.

If enough energy is available, the function returns the shift in the classical coordinates such that the new classical coordinate is z + shift and a Boolean equaling True if the hop has occurred. If not enough energy is available, the shift becomes zero and the Boolean is False.

Solves the equation:

\[H_{\mathrm{c}}(z) + \epsilon_{\mathrm{initial}} = H_{\mathrm{c}}(z + \mathrm{shift}) + \epsilon_{\mathrm{final}}\]

Keyword Args

zndarray: Current classical coordinate.
delta_zndarray: Rescaling direction.
eigval_difffloat: Energy difference between final and initial states.

Model Constants

None

Returns

shiftndarray: Shift in the classical coordinate.
hopbool: Whether the hop has occurred.

qclab.ingredients.init_classical_boltzmann_harmonic(model, parameters, **kwargs)#

Initialize classical coordinates according to Boltzmann statistics taking the classical coordinates to represent harmonic oscillators.

\(P(z) \propto \exp(-H_{\mathrm{c}}/k_{\mathrm{B}}T)\)

Keyword Args

seedndarray, int: Random seeds for each trajectory.

Model Constants

kBTfloat: Thermal quantum.
harmonic_frequencyndarray: Harmonic frequency of each classical coordinate.

Returns

zndarray: Complex classical coordinate.

qclab.ingredients.init_classical_wigner_harmonic(model, parameters, **kwargs)#

Initialize classical coordinates according to the Wigner distribution of the ground state of a harmonic oscillator.

Keyword Args

seedndarray, int: Random seeds for each trajectory.

Model Constants

kBTfloat: Thermal quantum.
harmonic_frequencyndarray: Harmonic frequency of each classical coordinate.

Returns

zndarray: Complex classical coordinate.

qclab.ingredients.init_classical_definite_position_momentum(model, parameters, **kwargs)#

Initialize classical coordinates with definite position and momentum.

Keyword Args

seedndarray, int: Random seeds for each trajectory.

Model Constants

init_positionndarray: Initial position of the classical coordinates.
init_momentumndarray: Initial momentum of the classical coordinates.

Returns

zndarray: Complex classical coordinate.

qclab.ingredients.init_classical_wigner_coherent_state(model, parameters, **kwargs)#

Initialize classical coordinates according to the Wigner distribution of a coherent state of a harmonic oscillator.

\(\vert a\rangle = \exp(a\hat{b}^{\dagger} - a^{*}\hat{b})\vert 0\rangle\)

where \(a\) is the complex displacement parameter of the coherent state. The Wigner distribution is given by

\(P(z) \propto \exp(-2\vert(z - a)\vert^{2})\).

Keyword Args

seedndarray, int: Random seeds for each trajectory.

Model Constants

coherent_state_displacementndarray: Complex displacement parameter of the coherent state.
harmonic_frequencyndarray: Harmonic frequency of each classical coordinate.

Returns

zndarray: Complex classical coordinate.

qclab.ingredients.rescaling_direction_random(model, parameters, **kwargs)#

Random rescaling direction function.

This function returns a random array for the rescaling direction. It is only included for documentation purposes.

Note that this is not a vectorized ingredient, as it is only called on the per-trajectory level.

Keyword Args

z_trajndarray: Current classical coordinate in the trajectory being rescaled.
init_state_indint: Index of the initial quantum state.
final_state_indint: Index of the final quantum state.

Model Constants

None

Returns

delta_zndarray: Direction in which to rescale coordinates.

qclab.ingredients.gauge_field_force_zero(model, parameters, **kwargs)#

Gauge field force function.

This function returns a zero array for the gauge field force. It is only included for documentation purposes.

Keyword Args

zndarray: Current classical coordinate.
state_indint: Index of the state for which the gauge field force is calculated.

Model Constants

None

Returns

gauge_forcendarray: Gauge field force.