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 the initialization Tasks, update Tasks, and collect Tasks 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. Algorithm objects define the transient quantities of an algorithm in the State object, which is a Python dictionary.

Algorithms in QC Lab 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 Model objects implemented without such a basis (and vice versa). As an example, the FewestSwitchesSurfaceHoppingAbInitio and MeanFieldAbInitio Algorithm objects can only be used with Model objects defined in an adiabatic basis. Ab initio Models can only be used with ab initio Algorithms, and vice versa. 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 a Constants object algorithm.settings which contains the settings specific to the Algorithm object. 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, users 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 update time step of the simulation, and the collect Recipe is executed at each collect time step 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.

Mean-field 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
 1class MeanField(Algorithm):
 2    """
 3    Mean-field dynamics algorithm class.
 4    """
 5
 6    def __init__(self, settings=None):
 7        if settings is None:
 8            settings = {}
 9        self.default_settings = {}
10        super().__init__(self.default_settings, settings)
11
12    initialization_recipe = [
13        tasks.initialize_variable_objects,
14        tasks.initialize_norm_factor,
15        tasks.initialize_z,
16        tasks.update_h_q_tot,
17    ]
18
19    update_recipe = [
20        # Begin RK4 integration steps.
21        partial(tasks.update_classical_force, z_name="z"),
22        tasks.update_quantum_classical_force,
23        tasks.update_z_rk4_k123,
24        partial(tasks.update_classical_force, z_name="z_1"),
25        partial(
26            tasks.update_quantum_classical_force,
27            z_name="z_1",
28            wf_changed=False,
29        ),
30        partial(tasks.update_z_rk4_k123, z_name="z", z_k_name="z_2", k_name="z_rk4_k2"),
31        partial(tasks.update_classical_force, z_name="z_2"),
32        partial(
33            tasks.update_quantum_classical_force,
34            z_name="z_2",
35            wf_changed=False,
36        ),
37        partial(
38            tasks.update_z_rk4_k123,
39            z_name="z",
40            z_k_name="z_3",
41            k_name="z_rk4_k3",
42            dt_factor=1.0,
43        ),
44        partial(tasks.update_classical_force, z_name="z_3"),
45        partial(
46            tasks.update_quantum_classical_force,
47            z_name="z_3",
48            wf_changed=False,
49        ),
50        tasks.update_z_rk4_k4,
51        # End RK4 integration steps.
52        tasks.update_wf_db_rk4,
53        tasks.update_h_q_tot,
54    ]
55
56    collect_recipe = [
57        tasks.update_t,
58        tasks.update_dm_db_wf,
59        tasks.update_quantum_energy_wf,
60        tasks.update_classical_energy,
61        tasks.collect_t,
62        tasks.collect_dm_db,
63        tasks.collect_classical_energy,
64        tasks.collect_quantum_energy,
65    ]

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
  1class FewestSwitchesSurfaceHopping(Algorithm):
  2    """
  3    Fewest switches surface hopping algorithm class.
  4    """
  5
  6    def __init__(self, settings=None):
  7        if settings is None:
  8            settings = {}
  9        self.default_settings = {
 10            "fssh_deterministic": False,
 11            "gauge_fixing": "sign_overlap",
 12            "use_gauge_field_force": False,
 13        }
 14        super().__init__(self.default_settings, settings)
 15
 16    initialization_recipe = [
 17        tasks.initialize_variable_objects,
 18        tasks.initialize_norm_factor,
 19        tasks.initialize_branch_seeds,
 20        tasks.initialize_z,
 21        tasks.update_h_q_tot,
 22        partial(
 23            tasks.diagonalize_matrix,
 24            matrix_name="h_q_tot",
 25            eigvals_name="eigvals",
 26            eigvecs_name="eigvecs",
 27        ),
 28        partial(
 29            tasks.update_eigvecs_gauge,
 30            gauge_fixing="phase_der_couple",
 31            eigvecs_previous_name="eigvecs",
 32        ),
 33        partial(tasks.copy_in_state, copy_name="eigvecs_previous", orig_name="eigvecs"),
 34        partial(
 35            tasks.update_vector_basis,
 36            input_vec_name="wf_db",
 37            basis_name="eigvecs",
 38            output_vec_name="wf_adb",
 39            adb_to_db=False,
 40        ),
 41        tasks.initialize_random_values_fssh,
 42        tasks.initialize_active_surface,
 43        tasks.initialize_dm_adb_0_fssh,
 44        tasks.update_act_surf_wf,
 45    ]
 46
 47    update_recipe = [
 48        partial(tasks.copy_in_state, copy_name="eigvecs_previous", orig_name="eigvecs"),
 49        # Begin RK4 integration steps.
 50        tasks.update_classical_force,
 51        partial(
 52            tasks.update_quantum_classical_force,
 53            wf_db_name="act_surf_wf",
 54            wf_changed=True,
 55        ),
 56        tasks.update_z_rk4_k123,
 57        partial(tasks.update_classical_force, z_name="z_1"),
 58        partial(
 59            tasks.update_quantum_classical_force,
 60            wf_db_name="act_surf_wf",
 61            z_name="z_1",
 62            wf_changed=False,
 63        ),
 64        partial(tasks.update_z_rk4_k123, z_name="z", z_k_name="z_2", k_name="z_rk4_k2"),
 65        partial(tasks.update_classical_force, z_name="z_2"),
 66        partial(
 67            tasks.update_quantum_classical_force,
 68            wf_db_name="act_surf_wf",
 69            z_name="z_2",
 70            wf_changed=False,
 71        ),
 72        partial(
 73            tasks.update_z_rk4_k123,
 74            z_name="z",
 75            z_k_name="z_3",
 76            k_name="z_rk4_k3",
 77            dt_factor=1.0,
 78        ),
 79        partial(tasks.update_classical_force, z_name="z_3"),
 80        partial(
 81            tasks.update_quantum_classical_force,
 82            wf_db_name="act_surf_wf",
 83            z_name="z_3",
 84            wf_changed=False,
 85        ),
 86        tasks.update_z_rk4_k4,
 87        # End RK4 integration steps.
 88        tasks.update_wf_db_propagator,
 89        tasks.update_h_q_tot,
 90        partial(
 91            tasks.diagonalize_matrix,
 92            matrix_name="h_q_tot",
 93            eigvals_name="eigvals",
 94            eigvecs_name="eigvecs",
 95        ),
 96        tasks.update_eigvecs_gauge,
 97        partial(
 98            tasks.update_vector_basis,
 99            input_vec_name="wf_db",
100            basis_name="eigvecs",
101            output_vec_name="wf_adb",
102            adb_to_db=False,
103        ),
104        tasks.update_hop_prob_fssh,
105        tasks.update_hop_inds_fssh,
106        tasks.update_hop_vals_fssh,
107        tasks.update_z_hop,
108        tasks.update_act_surf_hop,
109        tasks.update_act_surf_wf,
110    ]
111
112    collect_recipe = [
113        tasks.update_t,
114        tasks.update_dm_db_fssh,
115        tasks.update_quantum_energy_act_surf,
116        tasks.update_classical_energy_fssh,
117        tasks.collect_t,
118        tasks.collect_dm_db,
119        tasks.collect_quantum_energy,
120        tasks.collect_classical_energy,
121    ]

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.fewest_switches_surface_hopping (class FewestSwitchesSurfaceHoppingAbInitio).

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