% 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;