Models#

Models in QC Lab define the physics of the quantum-classical system under study. A model object is an instance of the qclab.Model class and is equipped with a set of constants and ingredients that specify the properties of the system in a manner that is agnostic to the quantum-classical algorithm being used.

At a minimum, the model object defines the Hamiltonian of the system:

\[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. These ingredients are discussed in detail in Ingredients.

The model object also contains a mandatory set of constants that define properties of the system:

  • num_quantum_states: the number of quantum states in the system,

  • num_classical_coordinates: the number of classical coordinates in the system,

  • classical_coordinate_mass: the mass of the classical coordinates,

  • classical_coordinate_weight: the weight of the classical coordinates (\(h\) in the complex-coordinate formalism).

The Model Class#

The model class is defined in the qclab.Model module. It is equipped with a constants object model.constants, an ingredients list model.ingredients, and a dictionary of default constants model.default_constants.

Constants#

Often, a model’s properties can be captured by a set of high-level constants that are suitable for user input. For example, the spin-boson model is defined by the following user-defined constants:

  • kBT: the thermal energy at a given temperature,

  • l_reorg: the reorganization energy of the bath,

  • E: the energy bias between the two quantum states,

  • V: the diabatic coupling between the two quantum states,

  • A: the number of bosonic modes in the bath,

  • W: the characteristic frequency of the bosonic modes in the bath.

  • boson_mass: the mass of each bosonic mode in the bath.

Each of these constants have a default value stored in the dictionary model.default_constants. At initialization, these defaults can be overwritten by passing a dictionary to the model constructor, as in:

from qclab.models import SpinBoson

# Create a dictionary of input constants to overwrite the defaults.
input_constants = {
    "kBT": 1.0,
    "l_reorg": 0.005,
    "E": 0.5,
    "V": 0.5,
    "A": 100,
    "W": 0.1,
    "boson_mass": 1.0
}
# Initialize the spin-boson model with the input constants.
model = SpinBoson(input_constants)

These input constants are first read into the model’s constants object model.constants which is an instance of the qclab.Constants class. Any input constants that are not specified will take on their default values. The input constants are then used to compute the mandatory constants required by QC Lab (specified above), as well as any additional constants that may be needed by the ingredients of the model. This computation is performed by a set of initialization ingredients that are typically unique to each model. The resulting “internal” constants are stored in the model’s constants object.

For example, the spin-boson model class uses the following initialization ingredients to compute its constants:

Note

When included within a class, the first argument of the ingredient is self instead of model. Here, specifying them outside of the class, we use model to refer to the instance of the model class.

def _init_h_q(model, parameters, **kwargs):
    """
    Initializes the constants required for the two-level quantum Hamiltonian.
    """
    model.constants.two_level_00 = model.constants.get("E")
    model.constants.two_level_11 = -model.constants.get("E")
    model.constants.two_level_01_re = model.constants.get("V")
    model.constants.two_level_01_im = 0
    return

def _init_h_qc(model, parameters, **kwargs):
    """
    Initializes the constants required for the diagonal linear quantum-classical Hamiltonian.
    """
    A = model.constants.get("A")
    l_reorg = model.constants.get("l_reorg")
    boson_mass = model.constants.get("boson_mass")
    h = model.constants.classical_coordinate_weight
    w = model.constants.harmonic_frequency
    model.constants.diagonal_linear_coupling = np.zeros((2, A))
    model.constants.diagonal_linear_coupling[0] = (
        w * np.sqrt(2.0 * l_reorg / A) * (1.0 / np.sqrt(2.0 * boson_mass * h))
    )
    model.constants.diagonal_linear_coupling[1] = (
        -w * np.sqrt(2.0 * l_reorg / A) * (1.0 / np.sqrt(2.0 * boson_mass * h))
    )
    return

def _init_h_c(model, parameters, **kwargs):
    """
    Initializes the constants required for the harmonic classical Hamiltonian.
    """
    A = model.constants.get("A")
    W = model.constants.get("W")
    model.constants.harmonic_frequency = W * np.tan(
        np.arange(0.5, A + 0.5, 1.0) * np.pi * 0.5 / A
    )
    return

def _init_model(model, parameters, **kwargs):
    """
    Initializes the mandatory constants required by QC Lab.
    """
    A = model.constants.get("A")
    boson_mass = model.constants.get("boson_mass")
    model.constants.num_classical_coordinates = A
    model.constants.num_quantum_states = 2
    model.constants.classical_coordinate_weight = model.constants.harmonic_frequency
    model.constants.classical_coordinate_mass = boson_mass * np.ones(A)
    return

For more information on the formatting of an ingredient, please refer to Ingredients. In the subsequent section we will discuss how these ingredients are included in a model class.

Ingredients List#

The ingredients in a model are contained in a list of tuples model.ingredients. Each tuple contains the name of the ingredient as a string and the ingredient function itself. For example, the spin-boson model includes the following ingredients:

ingredients = [
    ("h_q", ingredients.h_q_two_level),
    ("h_qc", ingredients.h_qc_diagonal_linear),
    ("h_c", ingredients.h_c_harmonic),
    ("dh_qc_dzc", ingredients.dh_qc_dzc_diagonal_linear),
    ("dh_c_dzc", ingredients.dh_c_dzc_harmonic),
    ("init_classical", ingredients.init_classical_boltzmann_harmonic),
    ("hop", ingredients.hop_harmonic),
    ("_init_h_q", _init_h_q),
    ("_init_h_qc", _init_h_qc),
    ("_init_model", _init_model),
    ("_init_h_c", _init_h_c),
]

As you can see, the ingredients list includes both the Hamiltonian ingredients (h_q, h_qc, h_c), their gradients (dh_qc_dzc, dh_c_dzc), as well as other ingredients used in the dynamics (init_classical, hop). Other ingredients define initialization steps that compute the model’s constants (_init_h_q, _init_h_qc, _init_h_c, _init_model). These are distinguished by their leading underscore, which indicates that they are to be run when the model is initialized.

To initialize the model’s constants manually one can run

model.initialize_constants()

which will execute all the ingredients in the list that begin with an underscore. After doing so, all the internal constants will be available in the model’s constants object model.constants. By default, this is done whenever a model object is initialized and whenever a constant is changed.

Importantly, a model’s ingredients list is executed from back to front. This means that one can add or overwrite an existing ingredient by appending a new tuple to the ingredients list. For example, if we wanted to change the quantum-classical coupling from diagonal to off-diagonal coupling, we could define a new ingredient and append it to the ingredients list:

from qclab import ingredients

def h_qc_off_diagonal(model, parameters, **kwargs):
    """
    A vectorized ingredient that couples the boson coordinates
    to the off-diagonal elements of the quantum Hamiltonian.
    """
    z = kwargs['z']
    A = model.constants.get("A")
    m = model.constants.classical_coordinate_mass
    h = model.constants.classical_coordinate_weight
    g = model.constants.diagonal_linear_coupling[0]
    N = model.constants.num_quantum_states
    batch_size = len(z)
    h_qc = np.zeros((batch_size, N, N), dtype=complex)
    h_qc[:, 0, 1] = g[np.newaxis, :] * (z + np.conj(z))
    h_qc[:, 1, 0] = np.conj(h_qc[:, 0, 1])
    return h_qc

# Overwrite the quantum-classical coupling ingredient.
model.ingredients.append(("h_qc", h_qc_off_diagonal))
# Overwrite the gradient of the quantum-classical coupling ingredient.
model.ingredients.append(("dh_qc_dzc", None))  # No analytical gradient available.

Spin-Boson Model#

Spin-Boson Model Constants#

Constant

Description

Default

kBT

Thermal energy.

1.0

V

Onsite energy.

0.5

E

Diabatic coupling.

0.5

A

Number of bosonic modes.

100

W

Characteristic frequency.

0.1

l_reorg

Reorganization energy.

0.005

boson_mass

Boson mass.

1.0
View full source
 1"""
 2This module contains the spin-boson Model class.
 3"""
 4
 5import numpy as np
 6from qclab.model import Model
 7from qclab import ingredients
 8
 9
10class SpinBoson(Model):
11    """
12    Spin-Boson model class.
13
14    Reference publication:
15    Tempelaar & Reichman. J. Chem. Phys. 148, 102309 (2018); https://doi.org/10.1063/1.5000843
16    """
17
18    def __init__(self, constants=None):
19        if constants is None:
20            constants = {}
21        self.default_constants = {
22            "kBT": 1.0,
23            "V": 0.5,
24            "E": 0.5,
25            "A": 100,
26            "W": 0.1,
27            "l_reorg": 0.005,
28            "boson_mass": 1.0,
29        }
30        super().__init__(self.default_constants, constants)
31        self.update_dh_qc_dzc = False
32        self.update_h_q = False
33
34    def _init_h_q(self, parameters, **kwargs):
35        self.constants.two_level_00 = self.constants.get("E")
36        self.constants.two_level_11 = -self.constants.get("E")
37        self.constants.two_level_01_re = self.constants.get("V")
38        self.constants.two_level_01_im = 0
39        return
40
41    def _init_h_qc(self, parameters, **kwargs):
42        A = self.constants.get("A")
43        l_reorg = self.constants.get("l_reorg")
44        boson_mass = self.constants.get("boson_mass")
45        h = self.constants.classical_coordinate_weight
46        w = self.constants.harmonic_frequency
47        self.constants.diagonal_linear_coupling = np.zeros((2, A))
48        self.constants.diagonal_linear_coupling[0] = (
49            w * np.sqrt(2.0 * l_reorg / A) * (1.0 / np.sqrt(2.0 * boson_mass * h))
50        )
51        self.constants.diagonal_linear_coupling[1] = (
52            -w * np.sqrt(2.0 * l_reorg / A) * (1.0 / np.sqrt(2.0 * boson_mass * h))
53        )
54        return
55
56    def _init_h_c(self, parameters, **kwargs):
57        A = self.constants.get("A")
58        W = self.constants.get("W")
59        self.constants.harmonic_frequency = W * np.tan(
60            np.arange(0.5, A + 0.5, 1.0) * np.pi * 0.5 / A
61        )
62        return
63
64    def _init_model(self, parameters, **kwargs):
65        A = self.constants.get("A")
66        boson_mass = self.constants.get("boson_mass")
67        self.constants.num_classical_coordinates = A
68        self.constants.num_quantum_states = 2
69        self.constants.classical_coordinate_weight = self.constants.harmonic_frequency
70        self.constants.classical_coordinate_mass = boson_mass * np.ones(A)
71        return
72
73    ingredients = [
74        ("h_q", ingredients.h_q_two_level),
75        ("h_qc", ingredients.h_qc_diagonal_linear),
76        ("h_c", ingredients.h_c_harmonic),
77        ("dh_qc_dzc", ingredients.dh_qc_dzc_diagonal_linear),
78        ("dh_c_dzc", ingredients.dh_c_dzc_harmonic),
79        ("init_classical", ingredients.init_classical_boltzmann_harmonic),
80        ("hop", ingredients.hop_harmonic),
81        ("_init_h_q", _init_h_q),
82        ("_init_h_qc", _init_h_qc),
83        ("_init_model", _init_model),
84        ("_init_h_c", _init_h_c),
85    ]

FMO Complex Model#

FMO Model Constants#

Constant

Description

Default

kBT

Thermal energy.

1

mass

Coordinate mass.

1

l_reorg

Reorganization energy.

35 cm -1

w_c

Characteristic frequency.

106.14 cm -1

N

Number of bosonic modes.

200
View full source
  1"""
  2This module contains the Model class for the Fenna-Matthews-Olson (FMO) complex.
  3"""
  4
  5import numpy as np
  6from qclab.model import Model
  7from qclab import ingredients
  8from qclab.numerical_constants import INVCM_TO_300K
  9
 10
 11class FMOComplex(Model):
 12    """
 13    A model representing the Fenna-Matthews-Olson (FMO) complex.
 14
 15    All quantities in this model are taken to be in units of kBT at 300 K.
 16
 17    At 300 K, kBT = 0.025852 eV = 208.521 cm^-1. Any quantity in wavenumbers
 18    is made unitless by dividing by kBT.
 19
 20    One unit of time is given by \hbar/kBT = 25.4595 fs.
 21
 22    Reference publication:
 23    Mulvihill et al. J. Chem. Phys. 154, 204109 (2021); https://doi.org/10.1063/5.0051101
 24    """
 25
 26    def __init__(self, constants=None):
 27        if constants is None:
 28            constants = {}
 29        self.default_constants = {
 30            "kBT": 1.0,
 31            "mass": 1.0,
 32            "l_reorg": 35.0 * INVCM_TO_300K,  # reorganization energy
 33            "w_c": 106.14 * INVCM_TO_300K,  # characteristic frequency
 34            "N": 200,
 35        }
 36        super().__init__(self.default_constants, constants)
 37        self.update_dh_qc_dzc = False
 38        self.update_h_q = False
 39
 40    def _init_model(self, parameters, **kwargs):
 41        """
 42        Initialize the model-specific constants.
 43        """
 44        self.constants.num_quantum_states = 7
 45        N = self.constants.get("N")
 46        w_c = self.constants.get("w_c")
 47        mass = self.constants.get("mass")
 48        self.constants.w = (
 49            w_c
 50            * np.tan(np.arange(0.5, N + 0.5, 1.0) * np.pi * 0.5 / (N))[np.newaxis, :]
 51            * np.ones((self.constants.num_quantum_states, N))
 52        ).flatten()
 53        self.constants.num_classical_coordinates = self.constants.num_quantum_states * N
 54
 55        self.constants.classical_coordinate_weight = self.constants.w
 56        self.constants.classical_coordinate_mass = mass * np.ones(
 57            self.constants.num_classical_coordinates
 58        )
 59        return
 60
 61    def _init_h_c(self, parameters, **kwargs):
 62        self.constants.harmonic_frequency = self.constants.w
 63        return
 64
 65    def _init_h_qc(self, parameters, **kwargs):
 66        N = self.constants.get("N")
 67        l_reorg = self.constants.get("l_reorg")
 68        m = self.constants.classical_coordinate_mass
 69        h = self.constants.classical_coordinate_weight
 70        w = self.constants.w
 71        self.constants.diagonal_linear_coupling = np.zeros(
 72            (
 73                self.constants.num_quantum_states,
 74                self.constants.num_classical_coordinates,
 75            )
 76        )
 77        for n in range(self.constants.num_quantum_states):
 78            self.constants.diagonal_linear_coupling[n, n * N : (n + 1) * N] = (
 79                w * np.sqrt(2.0 * l_reorg / N) * (1.0 / np.sqrt(2.0 * m * h))
 80            )[n * N : (n + 1) * N]
 81        return
 82
 83    def h_q(self, parameters, **kwargs):
 84        batch_size = kwargs["batch_size"]
 85        matrix_elements = (
 86            np.array(
 87                [
 88                    [12410, -87.7, 5.5, -5.9, 6.7, -13.7, -9.9],
 89                    [-87.7, 12530, 30.8, 8.2, 0.7, 11.8, 4.3],
 90                    [5.5, 30.8, 12210.0, -53.5, -2.2, -9.6, 6.0],
 91                    [-5.9, 8.2, -53.5, 12320, -70.7, -17.0, -63.3],
 92                    [6.7, 0.7, -2.2, -70.7, 12480, 81.1, -1.3],
 93                    [-13.7, 11.8, -9.6, -17.0, 81.1, 12630, 39.7],
 94                    [-9.9, 4.3, 6.0, -63.3, -1.3, 39.7, 12440],
 95                ],
 96                dtype=complex,
 97            )
 98            * INVCM_TO_300K
 99        )
100        # To reduce numerical errors we offset the diagonal elements by
101        # their minimum value.
102        matrix_elements = matrix_elements - np.min(
103            np.diag(matrix_elements)
104        ) * np.identity(7)
105        # Finally we broadcast the array to the desired shape
106        return np.broadcast_to(matrix_elements, (batch_size, 7, 7))
107
108    ingredients = [
109        ("h_q", h_q),
110        ("h_qc", ingredients.h_qc_diagonal_linear),
111        ("h_c", ingredients.h_c_harmonic),
112        ("dh_qc_dzc", ingredients.dh_qc_dzc_diagonal_linear),
113        ("dh_c_dzc", ingredients.dh_c_dzc_harmonic),
114        ("init_classical", ingredients.init_classical_boltzmann_harmonic),
115        ("hop", ingredients.hop_harmonic),
116        ("_init_h_qc", _init_h_qc),
117        ("_init_h_c", _init_h_c),
118        ("_init_model", _init_model),
119    ]

Tully Problem One#

Tully Problem One Model Constants#

Constant

Description

Default

init_momentum

Initial momentum.

10.0

init_position

Initial position.

-25.0

mass

Coordinate mass.

2000.0

A

See reference publication.

0.01

B

See reference publication.

1.6

C

See reference publication.

0.005

D

See reference publication.

1.0
View full source
  1"""
  2The module defines the Model class for Tully's first problem, a simple avoided crossing.
  3"""
  4
  5import numpy as np
  6from qclab.model import Model
  7from qclab import ingredients
  8from qclab import functions
  9
 10
 11class TullyProblemOne(Model):
 12    """
 13    Tully's first problem: a simple avoided crossing.
 14
 15    Reference publication:
 16    Tully. J. Chem. Phys. 93, 1061 (1990); https://doi.org/10.1063/1.459170
 17    """
 18
 19    def __init__(self, constants=None):
 20        if constants is None:
 21            constants = {}
 22        self.default_constants = {
 23            "init_momentum": 10.0,
 24            "init_position": -25.0,
 25            "mass": 2000.0,
 26            "A": 0.01,
 27            "B": 1.6,
 28            "C": 0.005,
 29            "D": 1.0,
 30        }
 31        super().__init__(self.default_constants, constants)
 32        self.update_dh_qc_dzc = True
 33        self.update_h_q = False
 34
 35    def _init_model(self, parameters, **kwargs):
 36        self.constants.num_quantum_states = 2
 37        self.constants.num_classical_coordinates = 1
 38        self.constants.classical_coordinate_mass = np.array(
 39            [self.constants.get("mass")]
 40        )
 41        self.constants.classical_coordinate_weight = np.array([1.0])
 42        self.constants.init_position = np.array([self.constants.get("init_position")])
 43        self.constants.init_momentum = np.array([self.constants.get("init_momentum")])
 44        return
 45
 46    def h_qc(self, parameters, **kwargs):
 47        """
 48        Quantum-classical Hamiltonian for Tully's first problem.
 49        """
 50        z = kwargs["z"]
 51        batch_size = len(z)
 52        num_quantum_states = self.constants.num_quantum_states
 53        A = self.constants.get("A")
 54        B = self.constants.get("B")
 55        C = self.constants.get("C")
 56        D = self.constants.get("D")
 57        # Calculate q.
 58        m = self.constants.classical_coordinate_mass[np.newaxis, :]
 59        h = self.constants.classical_coordinate_weight[np.newaxis, :]
 60        q = functions.z_to_q(z, m, h)[:, 0]
 61        # Calculate matrix elements.
 62        v_11 = np.zeros(batch_size, dtype=complex)
 63        v_11[q >= 0.0] = A * (1.0 - np.exp(-B * q))[q >= 0.0]
 64        v_11[q < 0.0] = -A * (1.0 - np.exp(B * q))[q < 0.0]
 65        v_12 = C * np.exp(-D * (q**2))
 66        # Assemble Hamiltonian.
 67        h_qc = np.zeros(
 68            (batch_size, num_quantum_states, num_quantum_states), dtype=complex
 69        )
 70        h_qc[:, 0, 0] = v_11
 71        h_qc[:, 0, 1] = v_12
 72        h_qc[:, 1, 0] = v_12
 73        h_qc[:, 1, 1] = -v_11
 74        return h_qc
 75
 76    def dh_qc_dzc(self, parameters, **kwargs):
 77        """
 78        Gradient w.r.t. to the conjugate z coordinate of the quantum-classical Hamiltonian
 79        for Tully's first problem.
 80        """
 81        z = kwargs["z"]
 82        batch_size = len(z)
 83        num_quantum_states = self.constants.num_quantum_states
 84        num_classical_coordinates = self.constants.num_classical_coordinates
 85        A = self.constants.get("A")
 86        B = self.constants.get("B")
 87        C = self.constants.get("C")
 88        D = self.constants.get("D")
 89        # Calculate q.
 90        m = self.constants.classical_coordinate_mass[np.newaxis, :]
 91        h = self.constants.classical_coordinate_weight[np.newaxis, :]
 92        q = functions.z_to_q(z, m, h)[:, 0]
 93        # Calculate phase-space gradients.
 94        dv_11_dq = np.zeros(batch_size, dtype=complex)
 95        dv_11_dq[q >= 0] = A * B * np.exp(-B * q)[q >= 0]
 96        dv_11_dq[q < 0] = A * B * np.exp(B * q)[q < 0]
 97        dv_12_dq = -2 * C * D * q * np.exp(-D * q**2)
 98        # Convert to complex gradients.
 99        dv_11_dzc = functions.dqdp_to_dzc(dv_11_dq, None, m[0], h[0])
100        dv_12_dzc = functions.dqdp_to_dzc(dv_12_dq, None, m[0], h[0])
101        # Assemble indices.
102        batch_idx = np.repeat(np.arange(batch_size, dtype=int), 4)
103        coord_idx = np.zeros(4 * batch_size, dtype=int)
104        state_i_idx = np.tile(np.array([0, 0, 1, 1], dtype=int), batch_size)
105        state_j_idx = np.tile(np.array([0, 1, 0, 1], dtype=int), batch_size)
106        inds = (batch_idx, coord_idx, state_i_idx, state_j_idx)
107        # Assemble matrix elements.
108        mels = np.empty(4 * batch_size, dtype=complex)
109        mels[0::4] = dv_11_dzc
110        mels[1::4] = dv_12_dzc
111        mels[2::4] = dv_12_dzc
112        mels[3::4] = -dv_11_dzc
113        # Assemble shape.
114        shape = (
115            batch_size,
116            num_classical_coordinates,
117            num_quantum_states,
118            num_quantum_states,
119        )
120        return inds, mels, shape
121
122    ingredients = [
123        ("h_q", ingredients.h_q_two_level),
124        ("h_qc", h_qc),
125        ("h_c", ingredients.h_c_free),
126        ("dh_qc_dzc", dh_qc_dzc),
127        ("dh_c_dzc", ingredients.dh_c_dzc_free),
128        ("init_classical", ingredients.init_classical_definite_position_momentum),
129        ("hop", ingredients.hop_free),
130        ("_init_model", _init_model),
131    ]