import numpy as np
import numba
from scipy.integrate import solve_ivp
from chemistrylab import material,vessel
from typing import NamedTuple, Tuple, Callable, Optional, List
from chemistrylab.reactions.reaction_info import ReactInfo
def _get_amounts(materials: Tuple[str], vessel: vessel.Vessel):
n=np.zeros(len(materials))
for i,key in enumerate(materials):
if key in vessel.material_dict:
n[i] = vessel.material_dict[key].mol
else:
n[i] = 0
return n
def _set_amounts(materials, solvents, material_classes, n, vessel):
for i,key in enumerate(materials):
amount = n[i]
if key in vessel.material_dict:
vessel.material_dict[key].mol=amount
elif n[i]>0:
mat = material_classes[i]()
mat.mol=amount
vessel.material_dict[key] = mat
vessel.validate_solvents()
vessel.validate_solutes()
#####################################################################################################################################
[docs]
@numba.jit(nopython=True)
def get_rates(stoich_coeff_arr, pre_exp_arr, activ_energy_arr, conc_coeff_arr, num_reagents, temp, conc):
"""
Finds the rate of reaction :math:`\\frac{dy}{dt}`
Args:
num_reagents (int): The number of reactants involved in the reaction
temp (float): The temperature of the reactions
conc (float): The initial concentrations of the materials
*_arr (np.array): See :class:`~chemistrylab.reactions.reaction_info.ReactInfo`
Returns:
np.array: Rates of change in concentration :math:`\\frac{dy}{dt}`.
"""
R = 8.314462619
conc = np.clip(conc, 0, None)
#k are the reaction constants
k = pre_exp_arr * np.exp((-1.0 * activ_energy_arr) / (R * temp))
rates = k*1
for i in range(len(rates)):
for j in range(num_reagents):
rates[i] *= conc[j] ** stoich_coeff_arr[i][j]
conc_change = np.zeros(conc_coeff_arr.shape[0])
for i in range(conc_change.shape[0]):
for j in range(rates.shape[0]):
conc_change[i] += conc_coeff_arr[i][j]*rates[j]
return conc_change
[docs]
@numba.njit
def newton_solve(stoich_coeff_arr, pre_exp_arr, activ_energy_arr, conc_coeff_arr, num_reagents, temp, conc, dt, N):
"""
Args:
num_reagents (int): The number of reactants involved in the reaction
temp (float): The temperature of the reactions
conc (float): The initial concentrations of the materials
dt (float): The amount of time to pass
N (int): The minimum number of time-steps to break dt into
*_arr (np.array): See :class:`~chemistrylab.reactions.reaction_info.ReactInfo`
Returns:
np.array: The final concentrations y(dt)
Solves the initial value problem :math:`\\frac{dy}{dt} = f(y)` specifically in the case where f(y) is chemical a rate calculation
(The problem is that we have :math:`y(t_0)` and need :math:`y(t_0+dt))`
This solver uses newton's method:
set :math:`T = \\frac{dt}{N}`
set :math:`y_0 = y(t_0)`
for n = 1,2,3,...N:
:math:`y_n = y_{n-1} + f(y_{n-1})*T`
:math:`y(dt) \\leftarrow y_N`
Intuitively, it is like taking a Riemann sum of dy/dt (but you get dy/dt by bootstrapping your current sum for y(t))
This implementation uses a variable step size in order to account for super fast-changing concentrations (wurtz distill)
"""
R = 8.314462619
#if your updates are below 5e-4 you can increase factor (I decided this is a good number)
targ = 5e-4
ddt=dt/N
k = (ddt*pre_exp_arr) * np.exp((-1.0 * activ_energy_arr) / (R * temp))
factor=1
count=0
d_conc=conc*0
while dt>0:
conc = np.clip(conc, 0, None)
#k are the reaction constants
rates = k*1
for i in range(len(rates)):
for j in range(num_reagents):
rates[i] *= conc[j] ** stoich_coeff_arr[i][j]
#calculate concentration changes
for i in range(conc.shape[0]):
d_conc[i]=0
for j in range(rates.shape[0]):
d_conc[i] += conc_coeff_arr[i][j]*rates[j]
#maximum proportion of material reduction
ratio=np.max(-d_conc/(conc+1e-6))
#mess with the step size to make sure you don't get any super huge concentration decreases
while ratio*factor<targ and factor<10:
factor*=2
while ratio*factor>0.1:
factor*=0.5
if factor*ddt>=dt:
factor = dt/ddt
dt=0
dt-=factor*ddt
count+=1
# Add concentration changes
conc+=d_conc*factor
return conc
[docs]
class Reaction():
def __init__(self,react_info: ReactInfo, solver: str = 'RK45', newton_steps: int = 100):
"""
A class to update concentrations of the materials in a vessel according to a reaction.
Args:
react_info (ReactInfo): Named Tuple containing all necessary reaction information
solver (str): Which solver to use
newton_steps (int): How many steps to use when the solver is 'newton'
"""
if not solver in {'RK45', 'RK23', 'DOP853', 'DBF', 'LSODA','newton'}:
solver='RK45'
self.solver=solver
self.newton_steps=newton_steps
#has to be set somewhere
self.threshold=1e-12
#materials we need for the reaction
self.reactants=react_info.REACTANTS
self.products=react_info.PRODUCTS
self.solvents=react_info.SOLVENTS
#Concatenate all of the materials (this should realistically be done in the reaction file since it has a direct
#Impact on the conc_coeff_arr
self.materials=react_info.MATERIALS
self.material_classes = tuple(material.REGISTRY[key] for key in self.materials)
#Necessary for calculating rates
self.stoich_coeff_arr = react_info.stoich_coeff_arr
self.pre_exp_arr = react_info.pre_exp_arr
self.activ_energy_arr = react_info.activ_energy_arr
self.conc_coeff_arr = react_info.conc_coeff_arr
self.num_reagents = len(self.reactants)
[docs]
def update_concentrations(self,vessel: vessel.Vessel, dt: float = 0):
"""
Takes in a vessel and applies the reaction to it, updating the material and solvent dicts in the process
Args:
vessel (Vessel): The vessel to perform the reaction on
dt (float): The amount of time passed during the reaction
"""
n = _get_amounts(self.materials, vessel)
if n.sum() < 1e-12:return
temperature = vessel.temperature
current_volume = vessel.filled_volume()
if dt==0:
dt = vessel.default_dt
#update concentrations
new_n = self.react(n, temperature, current_volume, dt)
#set the updated concentrations
_set_amounts(self.materials, self.solvents, self.material_classes, new_n, vessel)
[docs]
def react(self, n: np.array, temp: float, volume: float, dt: float):
"""
Args:
n (np.array): An array containing the amounts of each material in the vessel.
temp (float): The temperature of the system in Kelvin.
volume (float): The volume of the system in Litres.
dt (float): The time-step demarcating separate steps in seconds.
Returns:
np.array: The new amounts of each material in the vessel
"""
# set the intended vessel temperature in the differential equation module
self.temp = temp
conc=n/volume
if self.solver=='newton':
#newton solver should be faster but less accurate
new_conc = newton_solve(self.stoich_coeff_arr, self.pre_exp_arr,
self.activ_energy_arr, self.conc_coeff_arr,
self.num_reagents, self.temp, conc, dt, self.newton_steps)
else:
new_conc = solve_ivp(self, (0, dt), conc, method=self.solver).y[:, -1]
new_n = new_conc * volume
#set negligible amounts to 0
new_n *= (new_n > self.threshold)
return new_n
def __call__(self, t, conc):
"""
a function that calculates the change in concentration given a current concentration
remember to set the temperature before you call this function
This function is mainly used with the scipy ODE solvers
"""
return get_rates(self.stoich_coeff_arr, self.pre_exp_arr,
self.activ_energy_arr, self.conc_coeff_arr,
self.num_reagents, self.temp, conc)
NoneType = type(None)
[docs]
def react(vessel: vessel.Vessel, dt: float, other_vessel: NoneType, reaction: Reaction):
"""
:class:`~chemistrylab.vessel.Event` function to perform a reaction
"""
reaction.update_concentrations(vessel , dt)
return 0
vessel.Vessel.register(func = react, name = 'react')