Algorithms#

Algorithms in QC Lab define the sequence of operations that evolve the system defined by the model object (see Models) in time. They are composed of three recipes which define initialization steps, update steps, and collect steps that together define the desired algorithm. Each recipe is a list of “tasks” (see Tasks) which are executed in the order specified by the recipe list. Algorithms define the transient quantities of an algorithm in the state object, which is an instance of the Variable class (see Variable Objects).

Algorithm Objects#

Algorithm objects in QC Lab are instances of the qclab.Algorithm class. Each algorithm object is composed of three recipes: an initialization recipe algorithm.initialization_recipe, an update recipe algorithm.update_recipe, and a collect recipe algorithm.collect_recipe. Like a model object, an algorithm object has an instance of the Constants class algorithm.settings which contains the settings specific to the algorithm. Unlike the model object, algorithm objects do not have internal constants and so there is no initialization method as there is for model objects (see Models). Instead, the settings of the algorithm object are set directly by the user during or after instantiation of the algorithm object.

The empty Algorithm class is:

class Algorithm:
    """
    Algorithm class for defining and executing algorithm recipes.
    """

    def __init__(self, default_settings=None, settings=None):
        if settings is None:
            settings = {}
        if default_settings is None:
            default_settings = {}
        # Merge default settings with user-provided settings.
        settings = {**default_settings, **settings}
        # Construct a Constants object to hold settings.
        self.settings = Constants()
        # Put settings from the dictionary into the Constants object.
        for key, val in settings.items():
            setattr(self.settings, key, val)
        # Copy the recipes and output variables to ensure they are not shared
        # across instances.
        self.initialization_recipe = copy.deepcopy(self.initialization_recipe)
        self.update_recipe = copy.deepcopy(self.update_recipe)
        self.collect_recipe = copy.deepcopy(self.collect_recipe)

    initialization_recipe = []
    update_recipe = []
    collect_recipe = []

    def execute_recipe(self, sim, state, parameters, recipe):
        """
        Carry out the given recipe for the simulation by running
        each task in the recipe.
        """
        for func in recipe:
            state, parameters = func(sim, state, parameters)
        return state, parameters

After instantiating an algorithm object, the user can populate its recipes by assigning tasks to each recipe. For example, the mean-field algorithm can be defined from an empty Algorithm object as:

from qclab import Algorithm
import qclab.tasks as tasks
from functools import partial

# Create an empty algorithm object.
algorithm = Algorithm()
# Populate the initialization recipe.
algorithm.initialization_recipe = [
        tasks.initialize_variable_objects,
        tasks.initialize_norm_factor,
        tasks.initialize_z,
        tasks.update_h_quantum,
]
# Populate the update recipe.
algorithm.update_recipe = [
    # Begin RK4 integration steps.
    # RK4 steps excluded for brevity.
    # End RK4 integration steps.
    tasks.update_wf_db_rk4,
    tasks.update_h_q
]
# Populate the collect recipe.
algorithm.collect_recipe = [
    tasks.update_t,
    tasks.update_dm_db_mf,
    tasks.update_quantum_energy,
    tasks.update_classical_energy,
    tasks.collect_t,
    tasks.collect_dm_db,
    tasks.collect_classical_energy,
    tasks.collect_quantum_energy,
]

Each recipe is executed by the method algorithm.execute_recipe. The initialization recipe is executed once at the beginning of the simulation, the update recipe is executed at each time step of the simulation, and the collect recipe is executed once at the end of the simulation to gather and process results.

Variable Objects#

The tasks that make up an algorithm operate on the attributes of the state and parameters objects that are passed to each task. These objects are instances of the qclab.Variable class. The variable class is a simple container for attributes that are created and updated by tasks during the simulation. Importantly, when accessing an attribute that does not exist, a variable class returns None instead of raising an error. This allows tasks to check for the existence of attributes in the variable objects and create them if they do not exist.

Additionally, the contents of a variable class can be passed to it when instantiating it:

from qclab import Variable

# Create a variable object with some initial attributes.
state = Variable({"wf_db": None, "z": None, "p": None})

but can also be created empty.

By default, each variable object has an attribute output_dict which is an empty dictionary. This dictionary is used to gather results during the collect recipe of the algorithm. At the end of the simulation, the contents of this dictionary are returned as the output of the simulation. For example, if during the collect recipe a task adds an entry to the output dictionary as:

state.output_dict['dm_db'] = state.dm_db

then the output of the simulation will contain an entry with key 'dm_db' and value equal to the contents of the attribute state.dm_db averaged over its batch dimension.

Importantly, the parameters object is not used for the generation of outputs in QC Lab so its output_dict attribute is not used. The parameters object is instead used to interface with the ingredients of a model (see Models and Ingredients).

Mean Field Example#

As an example of a complete algorithm we include the source code for the mean-field algorithm below. This algorithm is defined in the qclab.algorithms.MeanField module and uses tasks from the qclab.tasks module to populate its recipes.

FSSH Collected Observables#

Key

Description

quantum_energy

The quantum energy of the system.

classical_energy

The classical energy of the system.

dm_db

The diabatic density matrix of the quantum subsystem.

t

The time points of the simulation.
View full source
 1"""
 2This module contains the MeanField algorithm class.
 3"""
 4
 5from functools import partial
 6from qclab.algorithm import Algorithm
 7from qclab import tasks
 8
 9
10class MeanField(Algorithm):
11    """
12    Mean-field dynamics algorithm class.
13    """
14
15    def __init__(self, settings=None):
16        if settings is None:
17            settings = {}
18        self.default_settings = {}
19        super().__init__(self.default_settings, settings)
20
21    initialization_recipe = [
22        tasks.initialize_variable_objects,
23        tasks.initialize_norm_factor,
24        tasks.initialize_z,
25        tasks.update_h_quantum,
26    ]
27
28    update_recipe = [
29        # Begin RK4 integration steps.
30        partial(tasks.update_classical_forces, z_name="z"),
31        tasks.update_quantum_classical_forces,
32        tasks.update_z_rk4_k123,
33        partial(tasks.update_classical_forces, z_name="z_1"),
34        partial(
35            tasks.update_quantum_classical_forces,
36            z_name="z_1",
37            wf_changed=False,
38        ),
39        partial(
40            tasks.update_z_rk4_k123, z_name="z", z_output_name="z_2", k_name="z_rk4_k2"
41        ),
42        partial(tasks.update_classical_forces, z_name="z_2"),
43        partial(
44            tasks.update_quantum_classical_forces,
45            z_name="z_2",
46            wf_changed=False,
47        ),
48        partial(
49            tasks.update_z_rk4_k123,
50            z_name="z",
51            z_output_name="z_3",
52            k_name="z_rk4_k3",
53            dt_factor=1.0,
54        ),
55        partial(tasks.update_classical_forces, z_name="z_3"),
56        partial(
57            tasks.update_quantum_classical_forces,
58            z_name="z_3",
59            wf_changed=False,
60        ),
61        tasks.update_z_rk4_k4,
62        # End RK4 integration steps.
63        tasks.update_wf_db_rk4,
64        tasks.update_h_quantum,
65    ]
66
67    collect_recipe = [
68        tasks.update_t,
69        tasks.update_dm_db_mf,
70        tasks.update_quantum_energy,
71        tasks.update_classical_energy,
72        tasks.collect_t,
73        tasks.collect_dm_db,
74        tasks.collect_classical_energy,
75        tasks.collect_quantum_energy,
76    ]

Surface Hopping Example#

As an additional example of a complete algorithm we include the source code for the fewest-switches surface hopping algorithm below. This algorithm is defined in the qclab.algorithms.FewestSwitchesSurfaceHopping module and uses tasks from the qclab.tasks module to populate its recipes.

FSSH Collected Observables#

Key

Description

quantum_energy

The quantum energy of the system.

classical_energy

The classical energy of the system.

dm_db

The diabatic density matrix of the quantum subsystem.

t

The time points of the simulation.
View full source
  1"""
  2This module contains the FewestSwitchesSurfaceHopping algorithm class.
  3"""
  4
  5from functools import partial
  6from qclab.algorithm import Algorithm
  7from qclab import tasks
  8
  9
 10class FewestSwitchesSurfaceHopping(Algorithm):
 11    """
 12    Fewest switches surface hopping algorithm class.
 13    """
 14
 15    def __init__(self, settings=None):
 16        if settings is None:
 17            settings = {}
 18        self.default_settings = {
 19            "fssh_deterministic": False,
 20            "gauge_fixing": "sign_overlap",
 21            "use_gauge_field_force": False,
 22        }
 23        super().__init__(self.default_settings, settings)
 24
 25    initialization_recipe = [
 26        tasks.initialize_variable_objects,
 27        tasks.initialize_norm_factor,
 28        tasks.initialize_branch_seeds,
 29        tasks.initialize_z,
 30        tasks.update_h_quantum,
 31        partial(
 32            tasks.diagonalize_matrix,
 33            matrix_name="h_quantum",
 34            eigvals_name="eigvals",
 35            eigvecs_name="eigvecs",
 36        ),
 37        partial(
 38            tasks.gauge_fix_eigs,
 39            gauge_fixing="phase_der_couple",
 40            eigvecs_previous_name="eigvecs",
 41        ),
 42        partial(tasks.copy_in_state, copy_name="eigvecs_previous", orig_name="eigvecs"),
 43        partial(
 44            tasks.basis_transform_vec,
 45            input_vec_name="wf_db",
 46            basis_name="eigvecs",
 47            output_vec_name="wf_adb",
 48            adb_to_db=False,
 49        ),
 50        tasks.initialize_random_values_fssh,
 51        tasks.initialize_active_surface,
 52        tasks.initialize_dm_adb_0_fssh,
 53        tasks.update_act_surf_wf,
 54    ]
 55
 56    update_recipe = [
 57        partial(tasks.copy_in_state, copy_name="eigvecs_previous", orig_name="eigvecs"),
 58        # Begin RK4 integration steps.
 59        tasks.update_classical_forces,
 60        partial(
 61            tasks.update_quantum_classical_forces,
 62            wf_db_name="act_surf_wf",
 63            wf_changed=True,
 64        ),
 65        tasks.update_z_rk4_k123,
 66        partial(tasks.update_classical_forces, z_name="z_1"),
 67        partial(
 68            tasks.update_quantum_classical_forces,
 69            wf_db_name="act_surf_wf",
 70            z_name="z_1",
 71            wf_changed=False,
 72        ),
 73        partial(
 74            tasks.update_z_rk4_k123, z_name="z", z_output_name="z_2", k_name="z_rk4_k2"
 75        ),
 76        partial(tasks.update_classical_forces, z_name="z_2"),
 77        partial(
 78            tasks.update_quantum_classical_forces,
 79            wf_db_name="act_surf_wf",
 80            z_name="z_2",
 81            wf_changed=False,
 82        ),
 83        partial(
 84            tasks.update_z_rk4_k123,
 85            z_name="z",
 86            z_output_name="z_3",
 87            k_name="z_rk4_k3",
 88            dt_factor=1.0,
 89        ),
 90        partial(tasks.update_classical_forces, z_name="z_3"),
 91        partial(
 92            tasks.update_quantum_classical_forces,
 93            wf_db_name="act_surf_wf",
 94            z_name="z_3",
 95            wf_changed=False,
 96        ),
 97        tasks.update_z_rk4_k4,
 98        # End RK4 integration steps.
 99        tasks.update_wf_db_eigs,
100        tasks.update_h_quantum,
101        partial(
102            tasks.diagonalize_matrix,
103            matrix_name="h_quantum",
104            eigvals_name="eigvals",
105            eigvecs_name="eigvecs",
106        ),
107        tasks.gauge_fix_eigs,
108        partial(
109            tasks.basis_transform_vec,
110            input_vec_name="wf_db",
111            basis_name="eigvecs",
112            output_vec_name="wf_adb",
113            adb_to_db=False,
114        ),
115        tasks.update_hop_probs_fssh,
116        tasks.update_hop_inds_fssh,
117        tasks.update_hop_vals_fssh,
118        tasks.update_z_hop_fssh,
119        tasks.update_act_surf_hop_fssh,
120        tasks.update_act_surf_wf,
121    ]
122
123    collect_recipe = [
124        tasks.update_t,
125        tasks.update_dm_db_fssh,
126        tasks.update_quantum_energy_fssh,
127        tasks.update_classical_energy_fssh,
128        tasks.collect_t,
129        tasks.collect_dm_db,
130        tasks.collect_quantum_energy,
131        tasks.collect_classical_energy,
132    ]