Rate Constants (except user defined ones)#

MICM supports a subset of the rate constants defined as part of the OpenAtmos Mechanism Configuration We will be adding more in the future. The links to the micm classes below detail the class and methods. Please check the OpenAtmos standard for specific on the algorithm implemented for each rate constant type. At present, supported rate constants are:

This tutorial covers all but the last one. See the User-Defined Rate Constants tutorial for examples and use cases on that.

We’ll setup and solve a fake chemical system with 7 species and 6 reactions,

\[\begin{split}A &\longrightarrow B, &k_{1, \mathrm{arrhenius}} \\ B &\longrightarrow C (\mathrm{alkoxy\ products}) + D (\mathrm{nitrate\ products}), &k_{2, \mathrm{branched}} \\ C &\longrightarrow E, &k_{3, \mathrm{surface}} \\ D &\longrightarrow 2F, &k_{4, \mathrm{ternary\ chemical\ activation}} \\ 2E &\longrightarrow G, &k_{5, \mathrm{troe}} \\ F &\longrightarrow G, &k_{6, \mathrm{tunneling}} \\\end{split}\]

MICM can be configured in two ways. We can either build the mechanism up by hand with the micm API, or parse a valid mechanism Configuration in the OpenAtmos format. In this tutorial, we will do both.

If you’re looking for a copy and paste, choose the appropriate tab below and be on your way! Otherwise, stick around for a line by line explanation.

// Each rate constant is in its own header file
#include <micm/process/arrhenius_rate_constant.hpp>
#include <micm/process/branched_rate_constant.hpp>
#include <micm/process/surface_rate_constant.hpp>
#include <micm/process/ternary_chemical_activation_rate_constant.hpp>
#include <micm/process/troe_rate_constant.hpp>
#include <micm/process/tunneling_rate_constant.hpp>
#include <micm/solver/rosenbrock.hpp>
#include <micm/solver/solver_builder.hpp>

#include <iomanip>
#include <iostream>
#include <map>

// Use our namespace so that this example is easier to read
using namespace micm;

int main(const int argc, const char* argv[])
{
  auto a = Species("A");
  auto b = Species("B");
  auto c = Species(
      "C",
      std::map<std::string, double>{ { "molecular weight [kg mol-1]", 0.025 },
                                     { "diffusion coefficient [m2 s-1]", 2.3e2 } });
  auto d = Species("D");
  auto e = Species("E");
  auto f = Species("F");
  auto g = Species("G");

  Phase gas_phase{ std::vector<Species>{ a, b, c, d, e, f, g } };

  Process r1 = Process::Create()
                   .SetReactants({ a })
                   .SetProducts({ Yields(b, 1) })
                   .SetRateConstant(ArrheniusRateConstant({ .A_ = 2.15e-1, .B_ = 0, .C_ = 110 }))
                   .SetPhase(gas_phase);

  // a branched reaction has two output pathways
  // this is represnted internal to micm as two different reactions
  auto branched_params = BranchedRateConstantParameters{ .X_ = 1.2, .Y_ = 204.3, .a0_ = 1.0e-3, .n_ = 2 };
  branched_params.branch_ = BranchedRateConstantParameters::Branch::Alkoxy;

  Process r2 = Process::Create()
                   .SetReactants({ b })
                   .SetProducts({ Yields(c, 1) })
                   .SetRateConstant(BranchedRateConstant(branched_params))
                   .SetPhase(gas_phase);

  branched_params.branch_ = BranchedRateConstantParameters::Branch::Nitrate;
  Process r3 = Process::Create()
                   .SetReactants({ b })
                   .SetProducts({ Yields(d, 1) })
                   .SetRateConstant(BranchedRateConstant(branched_params))
                   .SetPhase(gas_phase);

  // A surface rate constant also needs to know the effective radius and particle number concentration
  // we will set those later
  Process r4 = Process::Create()
                   .SetReactants({ c })
                   .SetProducts({ Yields(e, 1) })
                   .SetRateConstant(SurfaceRateConstant({ .label_ = "C", .species_ = c, .reaction_probability_ = 0.90 }))
                   .SetPhase(gas_phase);

  Process r5 = Process::Create()
                   .SetReactants({ d })
                   .SetProducts({ Yields(f, 2) })
                   .SetRateConstant(TernaryChemicalActivationRateConstant({ .k0_A_ = 1.2,
                                                                          .k0_B_ = 2.3,
                                                                          .k0_C_ = 302.3,
                                                                          .kinf_A_ = 2.6,
                                                                          .kinf_B_ = -3.1,
                                                                          .kinf_C_ = 402.1,
                                                                          .Fc_ = 0.9,
                                                                          .N_ = 1.2 }))
                   .SetPhase(gas_phase);

  // to have a stoichiemetric coefficient of more than one for reactants,
  // list the reactant that many times
  Process r6 = Process::Create()
                   .SetReactants({ e, e })
                   .SetProducts({ Yields(g, 1) })
                   .SetRateConstant(TroeRateConstant({ .k0_A_ = 1.2e4,
                                                     .k0_B_ = 167.0,
                                                     .k0_C_ = 3.0,
                                                     .kinf_A_ = 136.0,
                                                     .kinf_B_ = 5.0,
                                                     .kinf_C_ = 24.0,
                                                     .Fc_ = 0.9,
                                                     .N_ = 0.8 }))
                   .SetPhase(gas_phase);

  Process r7 = Process::Create()
                   .SetReactants({ f })
                   .SetProducts({ Yields(g, 1) })
                   .SetRateConstant(TunnelingRateConstant({ .A_ = 1.2, .B_ = 2.3, .C_ = 302.3 }))
                   .SetPhase(gas_phase);

  auto chemical_system = System(micm::SystemParameters{ .gas_phase_ = gas_phase });
  auto reactions = std::vector<micm::Process>{ r1, r2, r3, r4, r5, r6, r7 };

  auto solver = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
                    .SetSystem(chemical_system)
                    .SetReactions(reactions)
                    .Build();
  State state = solver.GetState();

  state.conditions_[0].temperature_ = 287.45;  // K
  state.conditions_[0].pressure_ = 101319.9;   // Pa

  state.SetConcentration(a, 1.0);  // mol m-3
  state.SetConcentration(b, 0.0);  // mol m-3
  state.SetConcentration(c, 0.0);  // mol m-3
  state.SetConcentration(d, 0.0);  // mol m-3
  state.SetConcentration(e, 0.0);  // mol m-3
  state.SetConcentration(f, 0.0);  // mol m-3
  state.SetConcentration(g, 0.0);  // mol m-3

  state.SetCustomRateParameter("C.effective radius [m]", 1e-7);
  state.SetCustomRateParameter("C.particle number concentration [# m-3]", 2.5e6);

  // choose a timestep and print the initial state
  double time_step = 500;  // s

  state.PrintHeader();
  state.PrintState(0);

  // solve for ten iterations
  for (int i = 0; i < 10; ++i)
  {
    // Depending on how stiff the system is
    // the solver integration step may not be able to solve for the full time step
    // so we need to track how much time the solver was able to integrate for and continue
    // solving until we finish
    double elapsed_solve_time = 0;
    solver.CalculateRateConstants(state);

    while (elapsed_solve_time < time_step)
    {
      auto result = solver.Solve(time_step - elapsed_solve_time, state);
      elapsed_solve_time = result.final_time_;
    }

    state.PrintState(time_step * (i + 1));
  }

  return 0;
}

// Each rate constant is in its own header file
#include <micm/configure/solver_config.hpp>
#include <micm/process/arrhenius_rate_constant.hpp>
#include <micm/process/branched_rate_constant.hpp>
#include <micm/process/surface_rate_constant.hpp>
#include <micm/process/ternary_chemical_activation_rate_constant.hpp>
#include <micm/process/troe_rate_constant.hpp>
#include <micm/process/tunneling_rate_constant.hpp>
#include <micm/solver/rosenbrock.hpp>
#include <micm/solver/solver_builder.hpp>

#include <iomanip>
#include <iostream>
#include <map>

// Use our namespace so that this example is easier to read
using namespace micm;

int main(const int argc, const char* argv[])
{
  SolverConfig solverConfig;

  std::string config_path = "./configs/rate_constants_no_user_defined";
  try
  {
    solverConfig.ReadAndParse(config_path);
  }
  catch (const std::system_error& e)
  {
    std::cerr << "Error reading and parsing config file: " << e.what() << std::endl;
    return 1;
  }

  micm::SolverParameters solver_params = solverConfig.GetSolverParams();

  auto chemical_system = solver_params.system_;
  auto reactions = solver_params.processes_;

  auto solver = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
                    .SetSystem(chemical_system)
                    .SetReactions(reactions)
                    .Build();

  State state = solver.GetState();

  state.conditions_[0].temperature_ = 287.45;  // K
  state.conditions_[0].pressure_ = 101319.9;   // Pa

  std::unordered_map<std::string, std::vector<double>> intial_concentration = {
    { "A", { 1.0 } },  // mol m-3
    { "B", { 0.0 } },  // mol m-3
    { "C", { 0.0 } },  // mol m-3
    { "D", { 0.0 } },  // mol m-3
    { "E", { 0.0 } },  // mol m-3
    { "F", { 0.0 } },  // mol m-3
    { "G", { 0.0 } },  // mol m-3
  };

  state.SetConcentrations(intial_concentration);

  state.SetCustomRateParameter("SURF.C surface.effective radius [m]", 1e-7);
  state.SetCustomRateParameter("SURF.C surface.particle number concentration [# m-3]", 2.5e6);

  // choose a timestep and print the initial state
  double time_step = 500;  // s

  state.PrintHeader();
  state.PrintState(0);

  // solve for ten iterations
  for (int i = 0; i < 10; ++i)
  {
    // Depending on how stiff the system is
    // the solver integration step may not be able to solve for the full time step
    // so we need to track how much time the solver was able to integrate for and continue
    // solving until we finish
    double elapsed_solve_time = 0;
    solver.CalculateRateConstants(state);

    while (elapsed_solve_time < time_step)
    {
      auto result = solver.Solve(time_step - elapsed_solve_time, state);
      elapsed_solve_time = result.final_time_;
    }

    state.PrintState(time_step * (i + 1));
  }

  return 0;
}

Line-by-line explanation#

To get started, we’ll need to include each rate constant type and the rosenbrock solver at the top of the file.

// Each rate constant is in its own header file
#include <micm/process/arrhenius_rate_constant.hpp>
#include <micm/process/branched_rate_constant.hpp>
#include <micm/process/surface_rate_constant.hpp>
#include <micm/process/ternary_chemical_activation_rate_constant.hpp>
#include <micm/process/troe_rate_constant.hpp>
#include <micm/process/tunneling_rate_constant.hpp>
#include <micm/solver/rosenbrock.hpp>
#include <micm/solver/solver_builder.hpp>

#include <iomanip>
#include <iostream>
#include <map>

After that, we’ll use the micm namespace so that we don’t have to repeat it everywhere we need it.

// Use our namespace so that this example is easier to read

To create a micm::RosenbrockSolver, we have to define a chemical system (micm::System) and our reactions, which will be a vector of micm::Process We will use the species to define these.

To do this by hand, we have to define all of the chemical species in the system. This allows us to set any properties of the species that may be necessary for rate constanta calculations, like molecular weights and diffusion coefficients for the surface reaction. We will also put these species into the gas phase.

{
  auto a = Species("A");
  auto b = Species("B");
  auto c = Species(
      "C",
      std::map<std::string, double>{ { "molecular weight [kg mol-1]", 0.025 },
                                     { "diffusion coefficient [m2 s-1]", 2.3e2 } });
  auto d = Species("D");
  auto e = Species("E");
  auto f = Species("F");
  auto g = Species("G");

Now that we have a gas phase and our species, we can start building the reactions. Two things to note are that stoichiemtric coefficients for reactants are represented by repeating that product as many times as you need. To specify the yield of a product, we’ve created a typedef micm::Yield and a function micm::Yields() that produces these.

  Process r1 = Process::Create()
                   .SetReactants({ a })
                   .SetProducts({ Yields(b, 1) })
                   .SetRateConstant(ArrheniusRateConstant({ .A_ = 2.15e-1, .B_ = 0, .C_ = 110 }))
                   .SetPhase(gas_phase);

  // a branched reaction has two output pathways
  // this is represnted internal to micm as two different reactions
  auto branched_params = BranchedRateConstantParameters{ .X_ = 1.2, .Y_ = 204.3, .a0_ = 1.0e-3, .n_ = 2 };
  branched_params.branch_ = BranchedRateConstantParameters::Branch::Alkoxy;

  Process r2 = Process::Create()
                   .SetReactants({ b })
                   .SetProducts({ Yields(c, 1) })
                   .SetRateConstant(BranchedRateConstant(branched_params))
                   .SetPhase(gas_phase);

  branched_params.branch_ = BranchedRateConstantParameters::Branch::Nitrate;
  Process r3 = Process::Create()
                   .SetReactants({ b })
                   .SetProducts({ Yields(d, 1) })
                   .SetRateConstant(BranchedRateConstant(branched_params))
                   .SetPhase(gas_phase);

  // A surface rate constant also needs to know the effective radius and particle number concentration
  // we will set those later
  Process r4 = Process::Create()
                   .SetReactants({ c })
                   .SetProducts({ Yields(e, 1) })
                   .SetRateConstant(SurfaceRateConstant({ .label_ = "C", .species_ = c, .reaction_probability_ = 0.90 }))
                   .SetPhase(gas_phase);

  Process r5 = Process::Create()
                   .SetReactants({ d })
                   .SetProducts({ Yields(f, 2) })
                   .SetRateConstant(TernaryChemicalActivationRateConstant({ .k0_A_ = 1.2,
                                                                          .k0_B_ = 2.3,
                                                                          .k0_C_ = 302.3,
                                                                          .kinf_A_ = 2.6,
                                                                          .kinf_B_ = -3.1,
                                                                          .kinf_C_ = 402.1,
                                                                          .Fc_ = 0.9,
                                                                          .N_ = 1.2 }))
                   .SetPhase(gas_phase);

  // to have a stoichiemetric coefficient of more than one for reactants,
  // list the reactant that many times
  Process r6 = Process::Create()
                   .SetReactants({ e, e })
                   .SetProducts({ Yields(g, 1) })
                   .SetRateConstant(TroeRateConstant({ .k0_A_ = 1.2e4,
                                                     .k0_B_ = 167.0,
                                                     .k0_C_ = 3.0,
                                                     .kinf_A_ = 136.0,
                                                     .kinf_B_ = 5.0,
                                                     .kinf_C_ = 24.0,
                                                     .Fc_ = 0.9,
                                                     .N_ = 0.8 }))
                   .SetPhase(gas_phase);

  Process r7 = Process::Create()
                   .SetReactants({ f })
                   .SetProducts({ Yields(g, 1) })
                   .SetRateConstant(TunnelingRateConstant({ .A_ = 1.2, .B_ = 2.3, .C_ = 302.3 }))

And finally we define our chemical system and reactions

  auto chemical_system = System(micm::SystemParameters{ .gas_phase_ = gas_phase });

After defining a valid OpenAtmos configuration with reactions that micm supports, configuring the chemical system and the processes is as simple as using the micm::SolverConfig class

{
  SolverConfig solverConfig;

  std::string config_path = "./configs/rate_constants_no_user_defined";
  try
  {
    solverConfig.ReadAndParse(config_path);
  }
  catch (const std::system_error& e)
  {
    std::cerr << "Error reading and parsing config file: " << e.what() << std::endl;
    return 1;
  }

Now that we have a chemical system and a list of reactions, we can create the RosenbrockSolver. There are several ways to configure the solver. Here we are using a three stage solver. More options can be found in the micm::RosenbrockSolverParameters and in the Solver configurations tutorial.

The rosenbrock solver will provide us a state, which we can use to set the concentrations, custom rate parameters, and temperature and pressure. Note that setting the custom rate paramters is different depending on if you define the configuration by hand or read it in. The parser has defaults for the names of the custom parameters and when defined by hand we choose these.

Initializing the state#

  auto solver = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
                    .SetSystem(chemical_system)
                    .SetReactions(reactions)
                    .Build();
  State state = solver.GetState();

  state.conditions_[0].temperature_ = 287.45;  // K
  state.conditions_[0].pressure_ = 101319.9;   // Pa

  state.SetConcentration(a, 1.0);  // mol m-3
  state.SetConcentration(b, 0.0);  // mol m-3
  state.SetConcentration(c, 0.0);  // mol m-3
  state.SetConcentration(d, 0.0);  // mol m-3
  state.SetConcentration(e, 0.0);  // mol m-3
  state.SetConcentration(f, 0.0);  // mol m-3

  auto chemical_system = solver_params.system_;
  auto reactions = solver_params.processes_;

  auto solver = micm::CpuSolverBuilder<micm::RosenbrockSolverParameters>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
                    .SetSystem(chemical_system)
                    .SetReactions(reactions)
                    .Build();

  State state = solver.GetState();

  state.conditions_[0].temperature_ = 287.45;  // K
  state.conditions_[0].pressure_ = 101319.9;   // Pa

  std::unordered_map<std::string, std::vector<double>> intial_concentration = {
    { "A", { 1.0 } },  // mol m-3
    { "B", { 0.0 } },  // mol m-3
    { "C", { 0.0 } },  // mol m-3
    { "D", { 0.0 } },  // mol m-3
    { "E", { 0.0 } },  // mol m-3

Finally, we are ready to pick a timestep and solve the system.

  state.SetCustomRateParameter("C.effective radius [m]", 1e-7);
  state.SetCustomRateParameter("C.particle number concentration [# m-3]", 2.5e6);

  // choose a timestep and print the initial state
  double time_step = 500;  // s

  state.PrintHeader();
  state.PrintState(0);

  // solve for ten iterations
  for (int i = 0; i < 10; ++i)
  {
    // Depending on how stiff the system is
    // the solver integration step may not be able to solve for the full time step
    // so we need to track how much time the solver was able to integrate for and continue
    // solving until we finish
    double elapsed_solve_time = 0;
    solver.CalculateRateConstants(state);

    while (elapsed_solve_time < time_step)
    {
      auto result = solver.Solve(time_step - elapsed_solve_time, state);
      elapsed_solve_time = result.final_time_;
    }

This is the output:

The Change of Concentration with Time#

time

A

B

C

D

E

F

G

0

1.00e+00

0.00e+00

0.00e+00

0.00e+00

0.00e+00

0.00e+00

0.00e+00

500

3.18e-09

3.66e-09

9.83e-01

3.88e-14

1.41e-03

2.02e-13

7.92e-03

1000

1.14e-14

1.31e-14

9.66e-01

1.39e-19

1.40e-03

7.24e-19

1.64e-02

1500

4.09e-20

4.71e-20

9.49e-01

4.98e-25

1.39e-03

2.59e-24

2.48e-02

2000

1.47e-25

1.69e-25

9.33e-01

1.79e-30

1.38e-03

9.30e-30

3.30e-02

2500

5.26e-31

6.05e-31

9.17e-01

6.40e-36

1.37e-03

3.33e-35

4.11e-02

3000

1.89e-36

2.17e-36

9.01e-01

2.30e-41

1.36e-03

1.20e-40

4.90e-02

3500

6.77e-42

7.78e-42

8.85e-01

8.23e-47

1.34e-03

4.29e-46

5.68e-02

4000

2.43e-47

2.79e-47

8.70e-01

2.95e-52

1.33e-03

1.54e-51

6.44e-02

4500

8.70e-53

1.00e-52

8.55e-01

1.06e-57

1.32e-03

5.51e-57

7.20e-02

5000

3.12e-58

3.59e-58

8.40e-01

3.80e-63

1.31e-03

1.98e-62

7.94e-02