A <-> B elementary unimolecular reaction¶with 1st-order kinetics in both directions, taken to equilibrium, using a simple, coarse fixed-timestep simulation.
In Part 2, below, we perform some analysis of the results: in particular, we examine the reaction rates.
(See also the experiment "1D/reactions/reaction_1" for a multi-compartment version)
LAST_REVISED = "Feb. 16, 2026"
LIFE123_VERSION = "1.0.0rc7" # Library version this experiment is based on
#import set_path # Using MyBinder? Uncomment this before running the next cell!
#import sys
#sys.path.append("C:/some_path/my_env_or_install") # CHANGE to the folder containing your venv or libraries installation!
# NOTE: If any of the imports below can't find a module, uncomment the lines above, or try: import set_path
import numpy as np
import ipynbname
from life123 import check_version, UniformCompartment, ThermoDynamics, PlotlyHelper
check_version(LIFE123_VERSION)
OK
# Initialize the HTML logging (for the graphics)
log_file = ipynbname.name() + ".log.htm" # Use the notebook base filename for the log file
# IN CASE OF PROBLEMS, set manually to any desired name
# Instantiate the simulator and specify the chemicals
uc = UniformCompartment()
# Reaction A <-> B , with 1st-order kinetics in both directions
uc.add_reaction(reactants="A", products="B", kF=3., kR=2.)
print("Number of reactions: ", uc.number_of_reactions())
add_reaction(): detected reaction type `ReactionUnimolecular` Number of reactions: 1
uc.describe_reactions()
Number of reactions: 1 0: A <-> B Elementary Unimolecular reaction (kF = 3 / kR = 2 / delta_G = -1,005.1 / K = 1.5 / Temp = 25 C) Chemicals involved in the above reactions: ['A', 'B']
# Send a plot of the network of reactions to the HTML log file
uc.plot_reaction_network(log_file=log_file)
[GRAPHIC ELEMENT SENT TO FILE `react_2_a.log.htm`]
# Initial concentrations of all the chemicals
uc.set_conc({"A": 10., "B": 50.})
uc.describe_state()
SYSTEM STATE at Time t = 0: 2 species: Species 0 (A). Conc: 10.0 Species 1 (B). Conc: 50.0 Chemicals involved in reactions: ['B', 'A']
uc.get_history()
| SYSTEM TIME | A | B | step | caption | |
|---|---|---|---|---|---|
| 0 | 0.0 | 10.0 | 50.0 | Set concentration |
uc.find_equilibrium_conc(rxn_index=0) # This is an EXACT equilibrium solution,
# for 1 reaction (the only reaction)
{'A': 24.0, 'B': 36.0}
That's consistent with the 3/2 ratio of forward/reverse rate constants (and the 1st order of the reaction)
B, relative to the small initial concentration of A¶More precisely, it's because the reaction quotient Q, at the current initial concentrations, is larger than our equilibrium constant K, which is 1.5 :
ThermoDynamics.compute_reaction_quotient(reactant_data=["A"], product_data=["B"],
conc={"A": 10., "B": 50.})
5.0
Now, let's see the reaction in action!
# First step of reaction
uc.single_compartment_react(initial_step=0.1, n_steps=1, variable_steps=False) # NOT using variable steps!
1 total fixed step(s) taken in 0.005 sec
uc.get_history()
| SYSTEM TIME | A | B | step | caption | |
|---|---|---|---|---|---|
| 0 | 0.0 | 10.0 | 50.0 | Set concentration | |
| 1 | 0.1 | 17.0 | 43.0 | 1 | last reaction step |
We can already see the reaction proceeding in reverse...
# Numerous more fixed steps
uc.single_compartment_react(initial_step=0.1, n_steps=10, variable_steps=False) # Again, fixed steps used
10 total fixed step(s) taken in 0.020 sec
uc.get_history()
| SYSTEM TIME | A | B | step | caption | |
|---|---|---|---|---|---|
| 0 | 0.0 | 10.000000 | 50.000000 | Set concentration | |
| 1 | 0.1 | 17.000000 | 43.000000 | 1 | last reaction step |
| 2 | 0.2 | 20.500000 | 39.500000 | 1 | 1st reaction step |
| 3 | 0.3 | 22.250000 | 37.750000 | 2 | |
| 4 | 0.4 | 23.125000 | 36.875000 | 3 | |
| 5 | 0.5 | 23.562500 | 36.437500 | 4 | |
| 6 | 0.6 | 23.781250 | 36.218750 | 5 | |
| 7 | 0.7 | 23.890625 | 36.109375 | 6 | |
| 8 | 0.8 | 23.945312 | 36.054688 | 7 | |
| 9 | 0.9 | 23.972656 | 36.027344 | 8 | |
| 10 | 1.0 | 23.986328 | 36.013672 | 9 | |
| 11 | 1.1 | 23.993164 | 36.006836 | 10 | last reaction step |
uc.get_system_conc() # The current concentrations, in the order the chemicals were added
array([23.99316406, 36.00683594])
That's very close to the values we previewed earlier: {'A': 24.0, 'B': 36.0}
# Verify that the reaction has reached equilibrium
uc.is_in_equilibrium()
0: A <-> B
Current concentrations: [A] = 23.99 ; [B] = 36.01
1. Ratio of reactant/product concentrations, adjusted for reaction orders: 1.50071
Formula used: [B] / [A]
2. Ratio of forward/reverse reaction rates: 1.5
Discrepancy between the two values: 0.04749 %
Reaction IS in equilibrium (within 1% tolerance)
True
uc.plot_history(colors=['darkturquoise', 'green'], show_intervals=True)
react_2_b this simulation gets repeated with an adaptive variable time resolution that takes smaller steps at the beginning, when the reaction is proceeding faster¶
There's no need to compute this; it is automatically saved during the reaction simulations
rates_df = uc.get_rate_history()
rates_df
| SYSTEM TIME | rxn0_rate | step | |
|---|---|---|---|
| 0 | 0.0 | -70.000000 | 0 |
| 1 | 0.1 | -35.000000 | 0 |
| 2 | 0.2 | -17.500000 | 1 |
| 3 | 0.3 | -8.750000 | 2 |
| 4 | 0.4 | -4.375000 | 3 |
| 5 | 0.5 | -2.187500 | 4 |
| 6 | 0.6 | -1.093750 | 5 |
| 7 | 0.7 | -0.546875 | 6 |
| 8 | 0.8 | -0.273438 | 7 |
| 9 | 0.9 | -0.136719 | 8 |
| 10 | 1.0 | -0.068359 | 9 |
Note that reaction rates are defined for the reaction products; since A is a reactant (in reaction 0, our only reaction), we must flip its sign; since the stoichiometry of A is simply 1, no further adjustment needed.
rates_df['A_dot'] = - rates_df['rxn0_rate'] # Add a column to the Pandas dataframe
rates_df
| SYSTEM TIME | rxn0_rate | step | A_dot | |
|---|---|---|---|---|
| 0 | 0.0 | -70.000000 | 0 | 70.000000 |
| 1 | 0.1 | -35.000000 | 0 | 35.000000 |
| 2 | 0.2 | -17.500000 | 1 | 17.500000 |
| 3 | 0.3 | -8.750000 | 2 | 8.750000 |
| 4 | 0.4 | -4.375000 | 3 | 4.375000 |
| 5 | 0.5 | -2.187500 | 4 | 2.187500 |
| 6 | 0.6 | -1.093750 | 5 | 1.093750 |
| 7 | 0.7 | -0.546875 | 6 | 0.546875 |
| 8 | 0.8 | -0.273438 | 7 | 0.273438 |
| 9 | 0.9 | -0.136719 | 8 | 0.136719 |
| 10 | 1.0 | -0.068359 | 9 | 0.068359 |
uc.get_history() # Revisited from earlier
| SYSTEM TIME | A | B | step | caption | |
|---|---|---|---|---|---|
| 0 | 0.0 | 10.000000 | 50.000000 | Set concentration | |
| 1 | 0.1 | 17.000000 | 43.000000 | 1 | last reaction step |
| 2 | 0.2 | 20.500000 | 39.500000 | 1 | 1st reaction step |
| 3 | 0.3 | 22.250000 | 37.750000 | 2 | |
| 4 | 0.4 | 23.125000 | 36.875000 | 3 | |
| 5 | 0.5 | 23.562500 | 36.437500 | 4 | |
| 6 | 0.6 | 23.781250 | 36.218750 | 5 | |
| 7 | 0.7 | 23.890625 | 36.109375 | 6 | |
| 8 | 0.8 | 23.945312 | 36.054688 | 7 | |
| 9 | 0.9 | 23.972656 | 36.027344 | 8 | |
| 10 | 1.0 | 23.986328 | 36.013672 | 9 | |
| 11 | 1.1 | 23.993164 | 36.006836 | 10 | last reaction step |
Notice that we lack a rate for the last time value, in the above table, because no reaction simulation starting at that time has been performed
p1 = uc.plot_history(chemicals="A", colors="darkturquoise") # The plot of [A] from the system history
p2 = PlotlyHelper.plot_pandas(df=rates_df, x_var="SYSTEM TIME", fields="A_dot", colors="brown") # The plot of A_dot, from rates_df
PlotlyHelper.combine_plots([p1, p2],
title="Concentration of A with time, and its rate of change (A_dot)",
y_label="[A] (turquoise) /<br> A_dot (brown)",
legend_title="Plot",
curve_labels=["A", "A_dot"])