% 10.34 - Fall 2006
% HW Set #11 - Kinetic Monte Carlo (Gillespie's Algorithm)
% Rob Ashcraft
% 12/10/2006
% This function performs the continuous solution and the kinetic Monte
% Carlo solution for several diameters and k-values. There are no inputs
% to this function, but they can be specified in the body of the function.
% Additional comment may be found in the file "p1_postprocessing.m".
function pset11p1_main
clear; clc;
warning off all;
Nav = 6.022e23;
time_0 = cputime;
line_style_vec = {'o','-','.','--','d',':','^','-.','x','+','s'};
color_vec = {'k','b','r','g','k','b','r','g','k','b','r','g'};
D_vec = [30; 50; 100; 250; 500; 1000; 2500];
k_vec = [1; 10; 100];
% initial conditions
C_ROO_0 = 0; % molar
C_ROOH_0 = 1e-6; % molar
C_RH_0 = 10; % molar, constant value
C_O2_0 = 1e-4; % molar, constant value
Temp = 313; % Kelvin, constant value
C_ROOH_cutoff = 0.010; % concentration at which to end simulation
% see additional comments for the continous part
disp('Initial number of molecules of: ROO, ROOH, RH, and O2 in droplet of size D');
for i=1:length(D_vec)
[N_ROO_0,N_ROOH_0,N_RH_0,N_O2_0] = calc_N(C_ROO_0,C_ROOH_0,C_RH_0,C_O2_0,D_vec(i));
disp(['Droplet size (nm): ',num2str(D_vec(i)),' -- ROO = ',num2str(N_ROO_0),...
', ROOH = ',num2str(N_ROOH_0),', RH = ',num2str(N_RH_0),', O2 = ',num2str(N_O2_0)]);
legend_D{i,1} = ['D = ',num2str(D_vec(i)),' nm'];
end
disp(' ');
% Solve the continuum problem at a number of different Diameters and k's
odesteps_small = 500;
odesteps_large = 10000;
t_small = linspace(0,3600*24,odesteps_small);
t_large = linspace(0,3600*24,odesteps_large);
options = [];
for j=1:length(k_vec)
%clear time C;
%time =
figure; hold on;
for i=1:length(D_vec)
[time(:,i,j),C(:,(i-1)*2+1:(i-1)*2+2,j)] = ode15s(@continuum_solver,t_large,[C_ROO_0,C_ROOH_0],options,Temp,D_vec(i),k_vec(j),C_RH_0,C_O2_0);
[time_plot(:,i,j),C_plot(:,(i-1)*2+1:(i-1)*2+2,j)] = ode15s(@continuum_solver,t_small,[C_ROO_0,C_ROOH_0],options,Temp,D_vec(i),k_vec(j),C_RH_0,C_O2_0);
[loc] = find(C(:,(i-1)*2+2,j) >= C_ROOH_cutoff);
time_to_failure(i,j) = time(loc(1),i,j);
rate_of_failure(i,j) = 1/time_to_failure(i,j);
plot(time_plot(:,i,j),C_plot(:,(i-1)*2+2,j),strcat(line_style_vec{i},color_vec{i}));
end
legend(legend_D,'Location','NorthWest');
axis([0 max(time_to_failure(:,j))*1.1 0 0.015]); %axis 'auto y';
end
figure;
semilogx(D_vec,time_to_failure(:,1),'o',D_vec,time_to_failure(:,2),'.',D_vec,time_to_failure(:,3),'*')
legend(['k = ',num2str(k_vec(1))],['k = ',num2str(k_vec(2))],['k = ',num2str(k_vec(3))]);
%================================================================
% Now work on the Kinetic Monte Carlo part of the problem:
% I have this set somewhat specfically for my computer in terms of timing.
% That is, I can specify a total wall time I want to spend per D, and the
% code will choose the appropriate number of steps. You compute could run
% faster or slower than mine (likely faster).
t_per_D = 60; % wall time you want to spend per droplet size, in seconds
% m gives the average time per MC trajectory run: Time (sec) = m*N_traj
% these values will different on a different computer
m = [1.85e-5;
8.88e-5;
0.0026;
0.113;
0.898;
7.188;
113.1];
% iterate over the different D-values
for i = 1:length(D_vec)
% compute the number of MC steps. This will be a non-round number, so
% the next part will attempt to round it to a nice, round number
MC_steps_temp(i) = t_per_D/m(i);
places(i) = length(num2str(round(MC_steps_temp(i)))); % number of place left of decimal
if(places(i) > 2) % if N >= 100
MC_steps(i) = ceil(MC_steps_temp(i)/(10^(places(i)-2))) * 10^(places(i)-2); % round up to the 2nd place (e.g. 121 -> 130, 4567304 -> 4600000)
elseif(places(i) == 2) % if 10 < N < 100
MC_steps(i) = ceil(MC_steps_temp(i)/10) * 10; % round up to the nearest 10 (94 -> 100, 18 -> 20)
else % N < 10
MC_steps(i) = ceil(MC_steps_temp(i)); % just round up (4.56 -> 5, 9.2123 -> 10)
end
end
j=2; % for k = 10
for i = 1:length(D_vec)
time_in = cputime; % this is for keeping track of wall time
disp(['Now running MC simulation for D = ',num2str(D_vec(i)),' nm']);
% For droplets that start with no ROOH or deplete all ROO + ROOH, just
% keep track of them with a counter (also for droplets go longer than
% my arbitrary upper time limit)
started_dead(i) = 0; % droplet that start with no ROOH
started_but_died(i) = 0; % droplet that start with 1+ ROOH, and get to a condition with no ROO + ROOH
went_too_long(i) = 0; % droplets don't reach critical concentration by upper time limit
[N_ROO,N_ROOH,N_RH_0,N_O2_0] = calc_N(C_ROO_0,C_ROOH_0,C_RH_0,C_O2_0,D_vec(i));
% defined vector sizes to make the code run faster. This is very
% important for large (>30000 elements) vectors where you are storing
% data. If you don't predefine the size, Matlab will grow the vector
% for you, but it becomes slow for large vectors. The idea here is to
% pre-allocate memory for the vector for the number of elements. If
% you start with more than 2 ROOH, basically every step will yield a
% finite time to failure. In this case the vector size is the number
% of MC steps. When you start with an average ROOH less than ~2, you
% can't prespecify the exact length, since some droplet will start dead
% or die along the way. In this case, I start with a vector slightly
% larger than I would expect to get, and fill it with NaN. This will
% easily allow me to extract the meaningful later on. I will alse make
% sure to fill these vectors from the top down, with no spaces between
% meaningful data.
if(N_ROOH > 2)
N_ROOH_temp{i} = zeros(MC_steps(i),1); % vector to store final number of ROOH (a long vector in each of 7 cells, one per D)
time_fail_temp{i} = zeros(MC_steps(i),1); % vector to store time to failure
rate_fail_temp{i} = zeros(MC_steps(i),1); % vector to store rate of failure
else
N_ROOH_temp{i} = zeros(MC_steps(i)*N_ROOH*1.25,1)*NaN;
time_fail_temp{i} = zeros(MC_steps(i)*N_ROOH*1.25,1)*NaN;
rate_fail_temp{i} = zeros(MC_steps(i)*N_ROOH*1.25,1)*NaN;
end
Vol = pi/6*D_vec(i)^3; % volume in nm^3
sizect = 0;
% "sizect" is a counter that will keep track of the droplets that reach
% the critical concentration. This will allow the vectors to be filled
% from the top down, even though the MC_step counter does not
% calculate the average number of molecules for a given drop size
[N_ROO_0,N_ROOH_0,N_RH_0,N_O2_0] = calc_N(C_ROO_0,C_ROOH_0,C_RH_0,C_O2_0,D_vec(i));
for mc = 1:MC_steps(i) % run trajectories for each D
% since the N_0's calculated from the concentrations will be non-integer,
% we need to round them to the nearest integer based on a probability.
% This is accomplished with the following lines:
N_ROO = round(N_ROO_0 + (rand-.5));
N_ROOH = round(N_ROOH_0 + (rand-.5));
N_RH_0 = round(N_RH_0 + (rand-.5));
N_O2_0 = round(N_O2_0 + (rand-.5));
if(N_ROOH == 0) % droplet start with nothing... do nothing
started_dead(i) = started_dead(i)+1; % keep track of the number that starts dead
mctime = NaN;
else % N_ROOH >= 1 ... something will happen, so run MC
% compute the rates of reation in units of molecules/sec-droplet
% I found that writing the rates out explicitly (instead of
% using a function call to a code that returns the rate speeded
% up the code by ~10 times (in the while loop mainly)
rates(1) = (1e15*exp(-15000/Temp)*N_ROOH/Vol)*Vol; % molecules/drop-s ROOH + 2_RH + 2_O2 ==> 2_ROO + ROH + H2O
rates(2) = (1e9*(1e24/6.022e23)*exp(-6000/Temp)*N_ROO/Vol*N_RH_0/Vol)*Vol; % molecules/drop-s ROO + RH + O2 ==> ROOH + ROO
rates(3) = (1e6*(1e24/6.022e23)*(N_ROO/Vol)*((N_ROO-1)/Vol))*Vol; % molecules/drop-s 2_ROO ==> ROH + aldehyde + O2
rates(4) = (0.10/D_vec(i)*k_vec(j)*N_ROO/Vol)*Vol; % D in nm % molecules/drop-s ROO + X(aq) ==> ROOH + X_dot(aq)
tau = 1/(rates(1) + rates(2) + rates(3) + rates(4)); % compute the Tau under these droplet conditions
mctime = 0; % initially, the time in the MC simulation is zero, this variable keeps track of it
C_ROOH = C_ROOH_0; % the concentration starts out at some initial value (not actually this, but this is not really used)
dead = 0; % this keep track of if the drop dies. if so, the while loop will break
% this "while" runs the MC until one of three conditions is
% satisfied:
% (1) The critical concentration is reached
% (2) The upper limit on time is reached
% (3) The droplet dies (all ROO and ROOH has reacted away)
while(C_ROOH <= C_ROOH_cutoff & mctime <= 1e8 & dead ~= 1)
if(N_ROO + N_ROOH == 0) % the droplet has died
dead = 1; % this will make the while loop break
N_ROOH = N_ROOH;
N_ROO = N_ROO;
elseif(N_ROO + N_ROOH ~= 0) % the droplet still has active species present... keep stepping along in time
rates(1) = (1e15*exp(-15000/Temp)*N_ROOH/Vol)*Vol; % molecules/drop-s ROOH + 2_RH + 2_O2 ==> 2_ROO + ROH + H2O
rates(2) = (1e9*(1e24/6.022e23)*exp(-6000/Temp)*N_ROO/Vol*N_RH_0/Vol)*Vol; % molecules/drop-s ROO + RH + O2 ==> ROOH + ROO
rates(3) = (1e6*(1e24/6.022e23)*(N_ROO/Vol)*((N_ROO-1)/Vol))*Vol; % molecules/drop-s 2_ROO ==> ROH + aldehyde + O2
rates(4) = (0.10/D_vec(i)*k_vec(j)*N_ROO/Vol)*Vol; % D in nm % molecules/drop-s ROO + X(aq) ==> ROOH + X_dot(aq)
tau = 1/(rates(1) + rates(2) + rates(3) + rates(4)); % recompute Tau at this condition
mctime = mctime - tau*log(rand); % determine the time until the next reaction event; add to current time to get new time
rand2 = rand;
%determine the fate next reaction that happens, and
%make the reaction happen by update N_x values
if(rand2 <= tau*rates(1)) % ROOH + 2_RH + 2_O2 ==> 2_ROO + ROH + H2O
N_ROOH = N_ROOH - 1;
N_ROO = N_ROO + 2;
elseif(rand2 <= tau*(rates(1) + rates(2))) % ROO + RH + O2 ==> ROOH + ROO
N_ROO = N_ROO - 0;
N_ROOH = N_ROOH + 1;
elseif(rand2 <= tau*(rates(1) + rates(2) + rates(3))) % 2_ROO ==> ROH + aldehyde + O2
N_ROO = N_ROO - 2;
N_ROOH = N_ROOH;
elseif(rand2 <= tau*(rates(1) + rates(2) + rates(3) + rates(4))) % ROO + X(aq) ==> ROOH + X_dot(aq)
N_ROO = N_ROO - 1;
N_ROOH = N_ROOH + 1;
end
end
C_ROOH = N_ROOH/((Nav*pi/6*D_vec(i)^3)/1e24); % compute the new concentration base on the new N_ROOH value after reaction
end % while loop
% catalogue the information from this trajectory into vectors
if(N_ROOH == 0 & mctime > 0)% & mctime <= 1e7) % started by died
started_but_died(i) = started_but_died(i) + 1;
elseif(C_ROOH < 0.01 & mctime > 1e8) % went too long
went_too_long(i) = went_too_long(i) + 1;
else
sizect = sizect + 1; % keeps track of the number of drops that reach critical concentration
time_fail_temp{i}(sizect,1) = mctime;
rate_fail_temp{i}(sizect,1) = 1/mctime;
N_ROOH_temp{i}(sizect,1) = N_ROOH;
end
end % end else N_ROOH ~= 0 statement
% This just prints out the status of the current droplet every
% 1/10th of the simulation... so you know things are still working
if(places(i) == 1)
disp(['Step Number ',num2str(mc),' of ',num2str(MC_steps(i)),'; Current wall-time = ',num2str(cputime-time_in),' seconds']);
else
if(rem(mc,MC_steps(i)/10) == 0)
disp(['Step Number ',num2str(mc),' of ',num2str(MC_steps(i)),'; Current wall-time = ',num2str(cputime-time_in),' seconds']);
end
end
end % end iterations of N_droplets
% Now we need to get rid of all of the NaN that we started with in our
% solution vectors, and did not fill in with real data. This finds the
% largest position that is not a NaN (where the last data point lies)
max_not_nan = max(find(isnan(N_ROOH_temp{i})==0));
% generate the final vectors to be used in the histogram
N_ROOH_final{i} = N_ROOH_temp{i}(1:max_not_nan);
time_fail_final{i} = time_fail_temp{i}(1:max_not_nan);
rate_fail_final{i} = rate_fail_temp{i}(1:max_not_nan);
% generate histograms based on the log-spaced bins
figure;
hist(time_fail_final{i}(:),logspace(3,8,50));
set(gca,'Xscale','log'); axis([1e3 1e8 0 1]); axis 'auto y';
figure;
hist(rate_fail_final{i}(:),logspace(-8,-3,50));
set(gca,'Xscale','log'); axis([1e-8 1e-3 0 1]); axis 'auto y';
% Print out some information on the stats for this diameter
alive(i) = MC_steps(i) - started_dead(i);
disp(' ');
disp(['Number of steps that reached critical C_ROOH: ',num2str(length(N_ROOH_final{i}(:))),' of ',num2str(MC_steps(i)),' steps']);
disp(['Went too long (but did not die): ',num2str(went_too_long(i)),' of ',num2str(alive(i)),' steps']);
disp(['Number of droplets that died (but started alive): ',num2str(started_but_died(i)),' of ',num2str(alive(i))]);
disp(['Fraction of droplets that started dead: ',num2str(started_dead(i)),' of ',num2str(MC_steps(i)),' (',num2str(started_dead(i)/MC_steps(i)*100),'%)']);
disp(' ');
disp(' ');
pause(1);
end
% displat the total time that the code took to run
time_total = cputime - time_0;
disp(['Total wall-time = ',num2str(time_total/3600),' hours']);
save('workspace_whatever');
% Other functions used in the code
%========================================
function dCdt = continuum_solver(time,C,T,D,k,C_RH,C_O2)
% this function returns the differential equations used to solve the
% continuum problem. Inputs are temperature in kelvin, D in nm, k in
% nm/sec, and concentrations in moles/L. Output is the dX/dt values in
% moles/L-sec
C_ROO = C(1);
C_ROOH = C(2);
rate1 = 1e15*exp(-15000/T)*C_ROOH; % mole/L-s ROOH + 2_RH + 2_O2 ==> 2_ROO + ROH + H2O
rate2 = 1e9*exp(-6000/T)*C_ROO*C_RH; % mole/L-s ROO + RH + O2 ==> ROOH + ROO
rate3 = 1e6*C_ROO^2; % mole/L-s 2_ROO ==> ROH + aldehyde + O2
rate4 = 0.10/D*k*C_ROO; % D in nm % mole/L-s ROO + X(aq) ==> ROOH + X_dot(aq)
dROO_dt = 2*rate1 - 2*rate3 - rate4;
dROOH_dt = -rate1 + rate2 + rate4;
dCdt = [dROO_dt; dROOH_dt];
return;
%=========================================
function [N_ROO,N_ROOH,N_RH,N_O2] = calc_N(C_ROO,C_ROOH,C_RH,C_O2,D)
% this function takes in the concentrations in moles/L and the droplet
% diameter in nm, and returns the number of molecules that would be
% expected in a single droplet of size D.
% diameter should be in nanometers and the concentration in moles/liter
Nav = 6.022e23;
% number of molecules in the droplet: N = C*Nav*pi/6*D^3*(1 dm^3/10^24 nm^3)
N_ROO = (C_ROO*Nav*pi/6*D^3)/1e24;
N_ROOH = (C_ROOH*Nav*pi/6*D^3)/1e24;
N_RH = (C_RH*Nav*pi/6*D^3)/1e24;
N_O2 = (C_O2*Nav*pi/6*D^3)/1e24;
return;