% Program to simulate a simple model of a unimolecular reaction with 5
% equally spaced energy levels, and reaction from the last level.
% Make sure that we don't start with any predefined values.
%clear
% Transition rates in s^{-1}
wdown = 1e9
wup = 0.8*wdown
wr = 1e13
% Number of molecules and initial populations:
t = 0;
Ntot = 10000
N = zeros(10,1);
N(1) = Ntot;
% A counter for the number of simulation steps
nsteps = 0;
% Main simulation loop: continue until we run out of molecules.
while (Ntot > 0)
% Calculate propensities.
a(1:4) = wup*N(1:4); % You can only go up from levels 1 to 4.
a(5:8) = wdown*N(2:5); % You can only go down from levels 2 to 5.
a(9) = wr*N(5); % Reaction from level 5.
% KMC selection of time and next jump.
a_sums = cumsum(a);
atot = a_sums(end);
r = rand(2,1);
Delta_t = -log(r(1))/atot;
k = find(a_sums > r(2)*atot,1,'first');
% What reaction did we pick? Change the populations accordingly.
if (k<5)
% This is an upward jump.
start_state = k;
N(start_state) = N(start_state) - 1;
N(start_state+1) = N(start_state+1) + 1;
elseif (k<9)
% Downward jump.
start_state = k - 3;
N(start_state) = N(start_state) - 1;
N(start_state-1) = N(start_state-1) + 1;
else % k = 9: reaction
N(5) = N(5) - 1;
Ntot = Ntot - 1;
end
% Increment t, then store the current total number of molecules.
t = t + Delta_t;
nsteps = nsteps + 1;
tval(nsteps) = t;
Ntotval(nsteps) = Ntot;
end
% Simulation done. Plot ln(Ntot) vs t.
% The last value of Ntotval will always be zero so we can't take its log.
% Get rid of this point.
nsteps = nsteps - 1;
tval = tval(1:nsteps);
Ntotval = Ntotval(1:nsteps);
plot(tval,log(Ntotval))
xlabel('t/s')
ylabel('ln(N_{total})','interpreter','tex')
hold % This will allow us to add a line to the plot.
% Fit the data to get a rate constant.
% line_coeffs(1) will be the slope and line_coeffs(2) the intercept.
line_coeffs = polyfit(tval,log(Ntotval),1);
% The first time point is at t=0. The last time point is stored in t.
t_endpoints = [0,t];
% Calculate the corresponding N values from the equation of the line:
N_endpoints = line_coeffs(1)*t_endpoints + line_coeffs(2);
plot(t_endpoints,N_endpoints)
hold off
% The rate constant is -(slope).
k = -line_coeffs(1)