Skip to content

Commit

Permalink
Merge pull request #3 from Jenny-Nelson-Group/mohammed
Browse files Browse the repository at this point in the history 
Mohammed
  • Loading branch information
@mohammedazzouzi15
mohammedazzouzi15 authored Jun 17, 2022
2 parents 0d315bf + cbf850b commit 82522ec
Show file tree
Hide file tree
Showing 7 changed files with 1,280 additions and 60 deletions.
2 changes: 1 addition & 1 deletion Example_workflow.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@

Prec = paramsRec; % initiliase the recombination parameters (default values)
offset = 0.25; % eV % energy difference between the excited state and the CT state
Prec.params.tickness = 100 * 1e-9; % m % thickness of the active layer
Prec.params.tickness = 10 * 1e-9; % m % thickness of the active layer
Prec.params.Ex.DG0 = 1.36;
Prec.params.CT.DG0 = Prec.params.Ex.DG0 - offset;
Prec.params.Ex.f = 2.56;
Expand Down
1,037 changes: 1,037 additions & 0 deletions classes/deviceparams.asv
View file

Large diffs are not rendered by default.

176 changes: 166 additions & 10 deletions classes/deviceparams.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -359,8 +359,8 @@
kbT=DP.physical_const.kB*DP.physical_const.T;%in eV
q=DP.physical_const.e;
offset=Prec.params.Ex.DG0-Prec.params.CT.DG0;
krecex=Prec.params.Ex.results.knr;
krecCT=Prec.params.CT.results.knr;
krecex=Prec.params.Ex.results.knr+Prec.params.Ex.results.krTot;
krecCT=Prec.params.CT.results.knr+Prec.params.CT.results.krTot;
RCTE=Prec.params.RCTE;
CT0=Prec.results.R0rad/q/1e3/(Prec.params.CT.results.krTot+Prec.params.Ex.results.krTot*exp(-offset/kbT)/RCTE);

Expand Down Expand Up @@ -470,21 +470,21 @@
DP.light_properties.Gensprofile_pos=Xpos*1e-7+DP.Layers{activelayer}.XL;
DP.light_properties.Gensprofile_signal=Gx;
end
function Dudt=kineticmodel(DP,activelayer,u,Gen,Ginj)
function Dudt=kineticmodel(DP,activelayer,u,Gen_ex,Gen_CT,Ginj)
kk=activelayer;
Dudt=[DP.Layers{kk}.kdis*(u(3))- DP.Layers{kk}.kfor*((u(1)*u(2)))+Ginj;
DP.Layers{kk}.kdis*(u(3))- DP.Layers{kk}.kfor*((u(1)*u(2)))+Ginj;
DP.Layers{kk}.kdisexc*(u(4))+DP.Layers{kk}.kfor*((u(1)*u(2)))-(DP.Layers{kk}.kdis*u(3)+DP.Layers{kk}.krec*(u(3)-DP.Layers{kk}.CT0))-DP.Layers{kk}.kforEx*(u(3));
Gen-DP.Layers{kk}.kdisexc*(u(4))-DP.Layers{kk}.krecexc*(u(4)-DP.Layers{kk}.Ex0)+DP.Layers{kk}.kforEx*(u(3));];%abs
Gen_CT+DP.Layers{kk}.kdisexc*(u(4))+DP.Layers{kk}.kfor*((u(1)*u(2)))-(DP.Layers{kk}.kdis*u(3)+DP.Layers{kk}.krec*(u(3)-DP.Layers{kk}.CT0))-DP.Layers{kk}.kforEx*(u(3));
Gen_ex-DP.Layers{kk}.kdisexc*(u(4))-DP.Layers{kk}.krecexc*(u(4)-DP.Layers{kk}.Ex0)+DP.Layers{kk}.kforEx*(u(3));];%abs
end
function [t,y]=solveKineticmodel(DP,Gen,Ginj)
function [t,y]=solveKineticmodel(DP,Gen_ex,Gen_CT,Ginj)
kk=2;
tspan = [0 5];
u0 = [DP.Layers{kk}.ni,DP.Layers{kk}.ni,DP.Layers{kk}.CT0,DP.Layers{kk}.Ex0];
[t,y] = ode15s(@(t,u) kineticmodel(DP,kk,u,Gen,Ginj), tspan, u0);
% y(end,4)/y(end,3);
% % figure
% % loglog(t,y,'-o')
% Gen_ex=Gen;
% Gen_CT=0;
[t,y] = ode15s(@(t,u) kineticmodel(DP,kk,u,Gen_ex,Gen_CT,Ginj), tspan, u0);

end
function [X,Y]=simulateTAS(DP,G,Gpulse,fighandle)
kk=2;
Expand Down Expand Up @@ -725,6 +725,162 @@ function simulate_TCSPC(DP,Prec,Background_Ex_Gen_rate,Laser_Ex_gen_rate,laser_w
xlabel('Time [ns]')
ylabel('Charge carrier density /pulse intensity [a.u]')
legend()


end
function [time,PL_emission,wavelength]=simulate_TCSPC_wavelength(DP,Prec,Background_Ex_Gen_rate,Laser_Ex_gen_rate,laser_width,exp_length,legendname)
%laser_width and exp_length in ns
%

kk=2;% active layer number
% Get the equilibrium solution first
tspan = [0 1e-1];
u0 = [DP.Layers{kk}.ni,DP.Layers{kk}.ni,DP.Layers{kk}.CT0,DP.Layers{kk}.Ex0];
[t,yeq] = ode15s(@(t,u) kineticmodel(DP,kk,u,Background_Ex_Gen_rate,0), tspan, u0);
u0 = yeq(end,:);
% run the solution with the laser on
tspan_laser = [0 laser_width*1e-9];
laser_ex_gen_func= @(t) Laser_Ex_gen_rate*exp(-t/1e-6);
[t_Laser,y_Laser] = ode15s(@(t,u) kineticmodel(DP,kk,u,Background_Ex_Gen_rate+laser_ex_gen_func(t),0), tspan_laser, u0);
u0 = y_Laser(end,:);
CTsum_Laser=y_Laser(:,3);
Exsum_Laser=y_Laser(:,4);
krE_Laser=Prec.params.CT.results.krE.*CTsum_Laser+Prec.params.Ex.results.krE.*Exsum_Laser;
% run the solution after the laser pulse
tspan = [0 exp_length*1e-9];
[t,y] = ode15s(@(t,u) kineticmodel(DP,kk,u,Background_Ex_Gen_rate,0), tspan, u0);
CTsum=y(:,3);
Exsum=y(:,4);
krE=Prec.params.CT.results.krE.*CTsum+Prec.params.Ex.results.krE.*Exsum;
[Maxlum,Peakpos]=max(krE');
[Maxlum_Laser,Peakpos_Laser]=max(krE_Laser');
ccdens_Laser=y_Laser(:,1);%/Laser_Ex_gen_rate/laser_width;
ccdens=y(:,1);%/Laser_Ex_gen_rate/laser_width;
t=t_Laser(end)+t;
PL_emission=[krE_Laser;krE(2:end,:)];
Enegy_distr=Prec.const.Edistribution;
time=[t_Laser*1e9;t(2:end)*1e9];
% figure
subplot(1,3,1)
title('PL intensity')
for wavelength=900:500:1200
energy=1240./wavelength;
Ydata=PL_emission(:,find(Enegy_distr>energy,1));
loglog(time,Ydata,"DisplayName",[num2str(wavelength) ' nm'])
hold on
end
hold on
% semilogy(t*1e9,(Maxlum)./max(Maxlum))
xlabel('Time [ns]')
ylabel('TCPSC signal [a.u]')
% ylim([Maxlum(end) 1.1*max(max([(Maxlum_Laser),(Maxlum)]))])
legend()
subplot(1,3,2)
title('PL intensity')
for wavelength=900:500:1200
energy=1240./wavelength;
Ydata=PL_emission(:,find(Enegy_distr>energy,1));
plot(time(time>50),Ydata(time>50)./max(Ydata(time>80)),"DisplayName",[num2str(wavelength) ' nm'])
hold on
end
hold on
% semilogy(t*1e9,(Maxlum)./max(Maxlum))
xlabel('Time [ns]')
ylabel('TCPSC signal [a.u]')
% ylim([Maxlum(end) 1.1*max(max([(Maxlum_Laser),(Maxlum)]))])
legend()

subplot(1,3,3)
title('normalised PL against wavelength')
colormap cool
for tt=logspace(0,log(max(time)*0.99)/log(10),6)
semilogy(1240./Prec.const.Edistribution,PL_emission(find(time>tt,1),:)/max(PL_emission(find(time>tt,1),:)),"DisplayName",[legendname ' @ ' num2str(tt,'%1.1e') ' ns'])
hold on
end
hold off
% semilogy(t*1e9,(Maxlum)./max(Maxlum))
xlabel('wavelength in [nm]')
ylabel('TCPSC signal [a.u]')
ylim([1e-4 1.1])
legend()

end
function [time,PL_emission,wavelength]=simulate_TCSPC_gaussianlaserpulse(DP,Prec,Background_Ex_Gen_rate,Laser_Ex_gen_rate_gauss,Laser_Ex_gen_rate_exp,...
gaussianwidth,exponentialdecay,laser_peak,exp_length,legendname)
%laser_width and exp_length in ns
%

kk=2;% active layer number
% Get the equilibrium solution first
tspan = [0 1e-1];
u0 = [DP.Layers{kk}.ni,DP.Layers{kk}.ni,DP.Layers{kk}.CT0,DP.Layers{kk}.Ex0];
[t,yeq] = ode15s(@(t,u) kineticmodel(DP,kk,u,Background_Ex_Gen_rate,0), tspan, u0);
u0 = yeq(end,:);
% run the solution with the laser on
tspan_laser = [0 exp_length*1e-9];
laser_ex_gen_func= @(t) Laser_Ex_gen_rate_exp*exp(-t/exponentialdecay)+Laser_Ex_gen_rate_gauss*exp(-(t-laser_peak*1e-9).^2./(gaussianwidth*1e-9)^2);
[t_Laser,y_Laser] = ode15s(@(t,u) kineticmodel(DP,kk,u,Background_Ex_Gen_rate+laser_ex_gen_func(t),0), tspan_laser, u0);
u0 = y_Laser(end,:);
CTsum_Laser=y_Laser(:,3);
Exsum_Laser=y_Laser(:,4);
krE_Laser=Prec.params.CT.results.krE.*CTsum_Laser+Prec.params.Ex.results.krE.*Exsum_Laser;
% run the solution after the laser pulse
[Maxlum_Laser,Peakpos_Laser]=max(krE_Laser');
ccdens_Laser=y_Laser(:,1);%/Laser_Ex_gen_rate/laser_width;
PL_emission=[krE_Laser];
Enegy_distr=Prec.const.Edistribution;
time=[t_Laser*1e9];
% figure
subplot(1,3,1)
title('PL intensity')
for wavelength=900:500:1200
energy=1240./wavelength;
Ydata=PL_emission(:,find(Enegy_distr>energy,1));
loglog(time,Ydata,"DisplayName",[num2str(wavelength) ' nm'])
hold on
end
hold on
% semilogy(t*1e9,(Maxlum)./max(Maxlum))
xlabel('Time [ns]')
ylabel('TCPSC signal [a.u]')
xlim([1e-1,exp_length])
% ylim([Maxlum(end) 1.1*max(max([(Maxlum_Laser),(Maxlum)]))])
legend()
subplot(1,3,2)
title('PL intensity')
for wavelength=900:500:1200
energy=1240./wavelength;
Ydata=PL_emission(:,find(Enegy_distr>energy,1));
loglog(time,Ydata./max(Ydata),"DisplayName",[num2str(wavelength) ' nm'])
hold on
end
hold on
loglog(time,laser_ex_gen_func(time*1e-9)./max(laser_ex_gen_func(time*1e-9)),'*-',"DisplayName",'normalised laser pulse intensity')

% semilogy(t*1e9,(Maxlum)./max(Maxlum))
xlabel('Time [ns]')
ylabel('TCPSC signal normalised [a.u]')
xlim([1e-1,exp_length])

ylim([1e-7 1.1])
legend()

subplot(1,3,3)
title('normalised PL against wavelength')
clist = colormap(jet(6));
jj=0;
for tt=logspace(0,log(max(time)*0.99)/log(10),6)
jj=jj+1;
semilogy(1240./Prec.const.Edistribution,PL_emission(find(time>tt,1),:)/max(PL_emission(find(time>tt,1),:)),"DisplayName",[legendname ' @ ' num2str(tt,'%1.1e') ' ns']...
,'color',clist(jj,:,:))
hold on
end
hold off
% semilogy(t*1e9,(Maxlum)./max(Maxlum))
xlabel('wavelength in [nm]')
ylabel('TCPSC signal [a.u]')
ylim([1e-4 1.1])
legend()

end
function [X,Y,Z]=simulate_PL(DP,Prec,fighandle)
Expand Down
20 changes: 17 additions & 3 deletions classes/dfplot.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -905,11 +905,25 @@ function Electroluminescence(DV,layernum,Voltage,Currentdensity)
disp("the JV did not run until the desired voltage")
end
else
if max(Current_density>Currentdensity*0.9 & Current_density<Currentdensity*1.1)
CTsum=mean(trapz(x(x>XL & x<XR), CT(Current_density>Currentdensity*0.9 & Current_density<Currentdensity*1.1,x>XL & x<XR), 2));
Exsum=mean(trapz(x(x>XL & x<XR), Ex(Current_density>Currentdensity*0.9 & Current_density<Currentdensity*1.1,x>XL & x<XR), 2));
if max(Current_density>Currentdensity )
try
CTsum=mean(trapz(x(x>XL & x<XR), CT(Current_density>Currentdensity*0.9 & Current_density<Currentdensity*1.1,x>XL & x<XR), 2));
Exsum=mean(trapz(x(x>XL & x<XR), Ex(Current_density>Currentdensity*0.9 & Current_density<Currentdensity*1.1,x>XL & x<XR), 2));
if isnan(CTsum)
error('could not find the right current density')
end
catch ME
disp( ME.message+"\n using an interpolation method ")
CTsum_h=mean(trapz(x(x>XL & x<XR), CT(find(Current_density>Currentdensity,1),x>XL & x<XR), 2));
Exsum_h=mean(trapz(x(x>XL & x<XR), Ex(find(Current_density>Currentdensity,1),x>XL & x<XR), 2));
CTsum_l=mean(trapz(x(x>XL & x<XR), CT(find(Current_density<Currentdensity,end),x>XL & x<XR), 2));
Exsum_l=mean(trapz(x(x>XL & x<XR), Ex(find(Current_density<Currentdensity,end),x>XL & x<XR), 2));
CTsum=(CTsum_h+CTsum_l)/2;
Exsum=(Exsum_l+Exsum_h)/2;
end
else
disp("the JV did not run until the desired current density")

end
end

Expand Down
21 changes: 12 additions & 9 deletions classes/paramsRec.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@
params.Ex.hW=0.15;%main vibronic energy considered in eV
params.Ex.sigma=0.001;%gaussian distribution for CT state
params.Ex.Dmus=3*3.33e-30/1.6e-19;%difference in static dipole moment (10 in DEbye ) then th
params.Ex.numbrestate=1;
params.Ex.numbrestate=2;
%% CT properties
params.CT.f=0.001;%oscillator strength of the CT state
params.CT.L0=0.14;%outer reorganisation energy(low frequency) in eV
Expand All @@ -49,7 +49,7 @@
params.CT.hW=0.15;%main vibronic energy considered in eV
params.CT.sigma=0.001;%gaussian distribution for CT state
params.CT.Dmus=10*3.33e-30/1.6e-19;%difference in static dipole moment (10 in DEbye ) then th
params.CT.numbrestate=1;
params.CT.numbrestate=2;
%%%%%%% add this to account for the effect of Hybredisation
params.Vstar=0.020; % Coupling between S1 and CT in eV
Prec.params=params;
Expand All @@ -63,7 +63,7 @@
for m=0:1:params.Number_Vibronic_Mode_final
for n=0:1:params.Number_Vibronic_Mode_initial

params.funlaguerre{m+1,n+1} = @(y) laguerreL(n,(m-n),y);
params.funlaguerre(m+1,n+1) = laguerreL(n,(m-n),params.S);
end
end
end
Expand Down Expand Up @@ -134,12 +134,15 @@
function params=FCWD(params,const)

%%%%%%%%%%%%%%%%Laguerre poly%%%%%%%%%%%%%%%%
laguerrecalc=zeros(params.Number_Vibronic_Mode_initial+1,params.Number_Vibronic_Mode_final+1);
for m=0:1:params.Number_Vibronic_Mode_final
for n=0:1:params.Number_Vibronic_Mode_initial
laguerrecalc(m+1,n+1)=params.funlaguerre{m+1,n+1}(params.S);
end
end

laguerrecalc=params.funlaguerre;
%zeros(params.Number_Vibronic_Mode_initial+1,params.Number_Vibronic_Mode_final+1);
% for m=0:1:params.Number_Vibronic_Mode_final
% for n=0:1:params.Number_Vibronic_Mode_initial
% laguerrecalc(m+1,n+1)=params.funlaguerre{m+1,n+1}(params.S);
% end
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
istate=0;

Expand Down
37 changes: 24 additions & 13 deletions scripts/Simulate_transient_PL.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,24 +8,25 @@
Prec.params.Ex.f = 2.56; % [1e-1, 10]%controls the strength of the optical transition from the Ex to ground
Prec.params.CT.f = 2e-2;%[1e-4, 1e-1] %controls the strength of the optical transition from the CT to ground

Prec.params.Ex.L0 = 0.06; %[0.5,0.2] % low frequency reorganisation energy % strongly affects the spectral width of the emission
Prec.params.Ex.Li = 0.045; %[0.5,0.2] % high frequency reorganisation energy % strongly affects the ratio between the vibronic peak emission
Prec.params.CT.L0 = 0.07; %[0.5,0.2]
Prec.params.CT.Li = 0.1; %[0.5,0.2]
Prec.params.Ex.L0 = 0.06; %[0.05,0.2] % low frequency reorganisation energy % strongly affects the spectral width of the emission
Prec.params.Ex.Li = 0.045; %[0.05,0.2] % high frequency reorganisation energy % strongly affects the ratio between the vibronic peak emission
Prec.params.CT.L0 = 0.05; %[0.05,0.2]
Prec.params.CT.Li = 0.04; %[0.05,0.2]
Prec.params.RCTE = 1e-2; % [1e-2,10] % ratio of CT to Exiton density of state


%%
figure
for kdis_exc=[1e11,5e10,1e10]
for Grate=5e25

% Generate a device with the defined parameters
% Parameters are from Prec which is defined above and from the PINDevice file, which is loaded below

activelayer = 2; % Active Layer Index % integer ( no need to change)
%NC = 2e19; % Number of Charge Carriers % cm^-3
Kfor = 1e-10; %[1e-8,1e-12] % Rate Constant CS to CT % cm^3 / s
kdis = 1e11; %[1e9,1e12] % Rate Constant CT dissociation % 1 / s
kdisex = kdis_exc; %[1e9,1e12] % Rate Constatn Ex dissociation % 1 / s
Kfor = 1e-11; %[1e-8,1e-12] % Rate Constant CS to CT % cm^3 / s
kdis = 1e9; %[1e9,1e12] % Rate Constant CT dissociation % 1 / s
kdisex = 1e9; %[1e9,1e12] % Rate Constatn Ex dissociation % 1 / s
%mobility = 5e-4; % Charge Carrier Mobility % cm^2 / V / s
%Generate the deviceparams class
deviceParameterFile = 'DeviceParameters_Default.xlsx';
Expand All @@ -35,13 +36,19 @@
%experiment parameter

Background_Ex_Gen_rate=1e10; %1sun would be around 1e22 % in cm-3 s-1
Laser_Ex_gen_rate=1e26; % in cm-3 s-1
laser_width= 5; % in ns
exp_length= 500; % in ns
%laser pulse properties
Laser_Ex_gen_rate_gauss=Grate; % in cm-3 s-1
Laser_Ex_gen_rate_exp=Grate*1e-6; % in cm-3 s-1

laser_peak= 2.5; % in ns
gaussianwidth=5;% in ns width of the gaussian
exponentialdecay=1e-9;%decay of the exponential component of the laser pulse in s

exp_length= 5001; % in ns
Temperature=300;%[30,350] % in Kelvin


exp_name=['Ex diss = ' num2str(kdis_exc,'%1.1e')];%,'%1.1e')];%
exp_name=['G rate = ' num2str(Laser_Ex_gen_rate,'%1.1e')];%,'%1.1e')];%
%doing the calculations

Prec.const.T = Temperature; % temperature in Kelvin
Expand All @@ -51,5 +58,9 @@
krecex = Prec.params.Ex.results.knr +Prec.params.Ex.results.krTot;
DP = DP.generateDeviceparams(2e19, activelayer, 5e-4, kdis, kdisex, Prec, Kfor, 0);

DP.simulate_TCSPC(Prec,Background_Ex_Gen_rate,Laser_Ex_gen_rate,laser_width,exp_length,exp_name)
%DP.simulate_TCSPC(Prec,Background_Ex_Gen_rate,Laser_Ex_gen_rate,laser_width,exp_length,exp_name)
[time,PL_emission,wavelength]=DP.simulate_TCSPC_gaussianlaserpulse(Prec,Background_Ex_Gen_rate,Laser_Ex_gen_rate_gauss,Laser_Ex_gen_rate_exp,...
gaussianwidth,exponentialdecay,laser_peak,exp_length,exp_name);
%simulate_TCSPC_gaussianlaserpulse(DP,Prec,Background_Ex_Gen_rate,Laser_Ex_gen_rate,...
% gaussianwidth,exponentialdecay,laserwidth,exp_length,legendname)
end
Loading

0 comments on commit 82522ec

Please sign in to comment.