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 a dictionary.

Algorithms in QC Lab can are tailored to Model objects defined in adiabatic or diabatic bases (see Models) in order to optimize their performance. Such tailoring breaks the compatibility between an algorithm implemented assuming a diabatic basis and those models implemented without such a basis (and vice versa). As an example, the AbInitioFewestSwitchesSurfaceHopping and AbInitioMeanField algorithms can only be used with models defined in an adiabatic basis. In most cases, model problems are defined in a diabatic basis and so we tailor the present adiabatic algorithms towards ab initio simulations which are the most common use case for an adiabatic basis.

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_q_tot,
]
# 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.

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_q_tot,
26    ]
27
28    update_recipe = [
29        # Begin RK4 integration steps.
30        partial(tasks.update_classical_force, z_name="z"),
31        tasks.update_quantum_classical_force,
32        tasks.update_z_rk4_k123,
33        partial(tasks.update_classical_force, z_name="z_1"),
34        partial(
35            tasks.update_quantum_classical_force,
36            z_name="z_1",
37            wf_changed=False,
38        ),
39        partial(tasks.update_z_rk4_k123, z_name="z", z_k_name="z_2", k_name="z_rk4_k2"),
40        partial(tasks.update_classical_force, z_name="z_2"),
41        partial(
42            tasks.update_quantum_classical_force,
43            z_name="z_2",
44            wf_changed=False,
45        ),
46        partial(
47            tasks.update_z_rk4_k123,
48            z_name="z",
49            z_k_name="z_3",
50            k_name="z_rk4_k3",
51            dt_factor=1.0,
52        ),
53        partial(tasks.update_classical_force, z_name="z_3"),
54        partial(
55            tasks.update_quantum_classical_force,
56            z_name="z_3",
57            wf_changed=False,
58        ),
59        tasks.update_z_rk4_k4,
60        # End RK4 integration steps.
61        tasks.update_wf_db_rk4,
62        tasks.update_h_q_tot,
63    ]
64
65    collect_recipe = [
66        tasks.update_t,
67        tasks.update_dm_db_wf,
68        tasks.update_quantum_energy_wf,
69        tasks.update_classical_energy,
70        tasks.collect_t,
71        tasks.collect_dm_db,
72        tasks.collect_classical_energy,
73        tasks.collect_quantum_energy,
74    ]
75
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_q_tot,
 31        partial(
 32            tasks.diagonalize_matrix,
 33            matrix_name="h_q_tot",
 34            eigvals_name="eigvals",
 35            eigvecs_name="eigvecs",
 36        ),
 37        partial(
 38            tasks.update_eigvecs_gauge,
 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.update_vector_basis,
 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_force,
 60        partial(
 61            tasks.update_quantum_classical_force,
 62            wf_db_name="act_surf_wf",
 63            wf_changed=True,
 64        ),
 65        tasks.update_z_rk4_k123,
 66        partial(tasks.update_classical_force, z_name="z_1"),
 67        partial(
 68            tasks.update_quantum_classical_force,
 69            wf_db_name="act_surf_wf",
 70            z_name="z_1",
 71            wf_changed=False,
 72        ),
 73        partial(tasks.update_z_rk4_k123, z_name="z", z_k_name="z_2", k_name="z_rk4_k2"),
 74        partial(tasks.update_classical_force, z_name="z_2"),
 75        partial(
 76            tasks.update_quantum_classical_force,
 77            wf_db_name="act_surf_wf",
 78            z_name="z_2",
 79            wf_changed=False,
 80        ),
 81        partial(
 82            tasks.update_z_rk4_k123,
 83            z_name="z",
 84            z_k_name="z_3",
 85            k_name="z_rk4_k3",
 86            dt_factor=1.0,
 87        ),
 88        partial(tasks.update_classical_force, z_name="z_3"),
 89        partial(
 90            tasks.update_quantum_classical_force,
 91            wf_db_name="act_surf_wf",
 92            z_name="z_3",
 93            wf_changed=False,
 94        ),
 95        tasks.update_z_rk4_k4,
 96        # End RK4 integration steps.
 97        tasks.update_wf_db_propagator,
 98        tasks.update_h_q_tot,
 99        partial(
100            tasks.diagonalize_matrix,
101            matrix_name="h_q_tot",
102            eigvals_name="eigvals",
103            eigvecs_name="eigvecs",
104        ),
105        tasks.update_eigvecs_gauge,
106        partial(
107            tasks.update_vector_basis,
108            input_vec_name="wf_db",
109            basis_name="eigvecs",
110            output_vec_name="wf_adb",
111            adb_to_db=False,
112        ),
113        tasks.update_hop_prob_fssh,
114        tasks.update_hop_inds_fssh,
115        tasks.update_hop_vals_fssh,
116        tasks.update_z_hop,
117        tasks.update_act_surf_hop,
118        tasks.update_act_surf_wf,
119    ]
120
121    collect_recipe = [
122        tasks.update_t,
123        tasks.update_dm_db_fssh,
124        tasks.update_quantum_energy_act_surf,
125        tasks.update_classical_energy_fssh,
126        tasks.collect_t,
127        tasks.collect_dm_db,
128        tasks.collect_quantum_energy,
129        tasks.collect_classical_energy,
130    ]
131

Ab Initio Surface Hopping Example#

As an example of an algorithm customized to Model objects defined in an adiabatic basis for compatibility with ab initio calculations, here we include the source code for the ab initio fewest-switches surface hopping algorithm implemented in the module qclab.algorithms.FewestSwitchesSurfaceHopping.

Ab initio FSSH Collected Observables#

Key

Description

quantum_energy

The quantum energy of the system.

classical_energy

The classical energy of the system.

dm_adb

The adiabatic density matrix of the quantum subsystem.

t

The time points of the simulation.
View full source
  1
  2class FewestSwitchesSurfaceHoppingAbInitio(Algorithm):
  3    """
  4    Fewest switches surface hopping algorithm class implemented in the adiabatic basis
  5    for compatibility with ab initio calculations.
  6    """
  7
  8    def __init__(self, settings=None):
  9        if settings is None:
 10            settings = {}
 11        self.default_settings = {
 12            "fssh_deterministic": False,
 13            "use_gauge_field_force": False,
 14            "update_wf_adb_eig_num_substeps": 10,
 15            "use_wf_overlaps_for_adb_connection": True,
 16        }
 17        super().__init__(self.default_settings, settings)
 18
 19    initialization_recipe = [
 20        tasks.initialize_variable_objects,
 21        partial(tasks.copy_to_parameters, state_name="seed", parameters_name="seed"),
 22        tasks.initialize_norm_factor,
 23        tasks.initialize_branch_seeds,
 24        tasks.initialize_z,
 25        partial(
 26            tasks.update_ab_initio_property,
 27            property_dict={
 28                "energy": {"z": "z", "excited_amplitudes": True},
 29                "gradient": {"z": "z", "state_inds_gradient": None},
 30                "derivative_coupling": {
 31                    "z": "z",
 32                    "state_inds_derivative_coupling": None,
 33                },
 34            },
 35        ),
 36        tasks.update_h_q_tot,
 37        tasks.update_classical_force,
 38        tasks.update_derivative_coupling_dzc,
 39        partial(tasks.update_quantum_classical_force, wf_db_name="wf_adb"),
 40        partial(tasks.update_adb_connection, update_derivative_coupling=False),
 41        tasks.initialize_random_values_fssh,
 42        tasks.initialize_active_surface,
 43        tasks.initialize_dm_adb_0_fssh,
 44        partial(
 45            tasks.diagonalize_matrix,
 46            matrix_name="h_q_tot",
 47            eigvals_name="eigvals",
 48            eigvecs_name="eigvecs",
 49        ),
 50        tasks.update_act_surf_wf,
 51        tasks.update_quantum_energy_act_surf,
 52        tasks.update_classical_energy_fssh,
 53    ]
 54
 55    update_recipe = [
 56        partial(
 57            tasks.copy_in_state,
 58            copy_name="aip_excited_amplitudes_previous",
 59            orig_name="aip_excited_amplitudes",
 60        ),
 61        partial(
 62            tasks.copy_in_state,
 63            copy_name="eigvecs_previous",
 64            orig_name="eigvecs",
 65        ),
 66        partial(
 67            tasks.copy_in_state,
 68            copy_name="adb_connection_previous",
 69            orig_name="adb_connection",
 70        ),
 71        partial(tasks.copy_in_state, copy_name="h_q_tot_previous", orig_name="h_q_tot"),
 72        partial(
 73            tasks.copy_in_state,
 74            copy_name="quantum_classical_force_previous",
 75            orig_name="quantum_classical_force",
 76        ),
 77        partial(
 78            tasks.copy_in_state,
 79            copy_name="classical_force_previous",
 80            orig_name="classical_force",
 81        ),
 82        partial(
 83            tasks.copy_in_state,
 84            copy_name="z_previous",
 85            orig_name="z",
 86        ),
 87        tasks.update_q_velocity_verlet,
 88        partial(
 89            tasks.update_ab_initio_property,
 90            property_dict={
 91                "energy": {"z": "z", "excited_amplitudes": True},
 92                "wf_overlaps": {
 93                    "z": "z",
 94                    "z_previous": "z_previous",
 95                    "amplitudes_previous": "aip_excited_amplitudes_previous",
 96                    "amplitudes_current": "aip_excited_amplitudes",
 97                },
 98            },
 99        ),
100        # tasks.update_adb_connection,
101        partial(tasks.update_adb_connection, update_derivative_coupling=True),
102        tasks.update_h_q_tot,
103        partial(
104            tasks.diagonalize_matrix,
105            matrix_name="h_q_tot",
106            eigvals_name="eigvals",
107            eigvecs_name="eigvecs",
108        ),
109        partial(
110            tasks.update_wf_adb_hop_prob,
111            update_hopping_probabilities=True,
112        ),
113        partial(
114            tasks.update_hop_inds_fssh,
115            hop_bool_name="hop_bool",
116            hop_pairs_name="hop_pairs",
117        ),
118        partial(
119            tasks.update_ab_initio_property,
120            property_dict={
121                "derivative_coupling": {
122                    "calc_property": "hop_bool",
123                    "z": "z",
124                    "state_inds_derivative_coupling": "hop_pairs",
125                },
126            },
127        ),
128        tasks.update_derivative_coupling_dzc,
129        partial(
130            tasks.update_hop_vals_fssh,
131            derivative_coupling_dzc_name="derivative_coupling_dzc",
132        ),
133        tasks.update_z_hop,
134        tasks.update_act_surf_hop,
135        tasks.update_act_surf_wf,
136        partial(
137            tasks.update_ab_initio_property,
138            property_dict={
139                "gradient": {"z": "z", "state_inds_gradient": "act_surf_ind"},
140            },
141        ),
142        partial(tasks.update_quantum_classical_force, wf_db_name="act_surf_wf"),
143        tasks.update_p_velocity_verlet,
144        tasks.update_classical_force,
145    ]
146
147    collect_recipe = [
148        tasks.update_t,
149        tasks.update_dm_db_fssh,
150        tasks.update_quantum_energy_act_surf,
151        tasks.update_classical_energy_fssh,
152        tasks.collect_t,
153        partial(tasks.collect_dm_db, dm_db_name="dm_adb", dm_db_output_name="dm_adb"),
154        tasks.collect_quantum_energy,
155        tasks.collect_classical_energy,
156    ]