%Model name: model of Memristor based on Theoretical formulas % this model was written by Dmitry Fliter and Keren Talisveyberg % Technion Israel institute of technology EE faculty December 2011 model = 1; % define the model 0 - Linear Ion Drift ; 1 - Simmons Tunnel Barrier; 2 - Team model ; 3 - Nonlinear Ion Drift model win = 0; % define the window type : 0 - No window; 1 - Jogelkar window ; 2 - Biolek window ; 3 - Prodromakis window ; 4- Kvatinsky window (Team model only) iv = 0; % IV_relation=0 means linear V=IR. IV_relation=1 means nonlinear V=I*exp{..} %Genaral parameters num_of_cycles = 20; amp = 1; freq = 2e6; w_init = 0.5; % the initial state condition [0:1] D = 3e-09; V_t = 0; P_coeff = 2; J = 1; Roff = 2e5; Ron = 100; %Linear Ion Drift parameters uV=1e-15; %%dopants mobility %Simmons Tunnel Bariier % Team parameters a_on = 2e-09; a_off = 1.2e-09; c_on = 40e-06; c_off = 3.5e-06; alpha_on = 3; alpha_off = 3; k_on = -8e-5; k_off = 8e-5; i_on = 8.9e-06; i_off = 115e-06; v_on=-0.02; v_off=0.02; x_on = 3e-09; x_off = 0; X_c = 107e-12; b = 500e-06; %Nonlinear Ion Drift parametrs beta = 9; a = 4; c = 0.01; n = 14; q = 13; g = 4; alpha = 7; %% Linear Ion Drift model if (model==0) tspan=[0 num_of_cycles/freq]; %%time length of the simulation points=2e5; %%number of sampling points W0=w_init*D; %define the initial value of W tspan_vector = linspace(tspan(1),tspan(2),points); % Create vector of initial values I = amp*sin(freq*2*pi*tspan_vector); %%can also use square wave generated by : (square(tspan_vector)); W=zeros(size((tspan_vector))); W_dot=zeros(size((tspan_vector))); delta_t=tspan_vector(2)-tspan_vector(1); %%define the step size W(1)=W0; %% initiliaze the first W vetor elemnt to W0 - the initial condition for i=2:length(tspan_vector) % case this is an ideal window if (((win==0) || (win==4)) && ((abs (I(i))) >= (V_t/ (Ron*W(i-1)/D+Roff*(1-W(i-1)/D))))) W_dot(i)=I(i)*(Ron*uV/D); W(i)=W(i-1)+W_dot(i)*delta_t; elseif ((win==0) && ((abs(I(i))) < (V_t/ (Ron*W(i-1)/D+Roff*(1-W(i-1)/D))))) W(i)=W(i-1); W_dot(i)=0; end % case this is Jogelkar window if ((win==1) && ((abs(I(i))) >= (V_t/(Ron*W(i-1)/D+Roff*(1-W(i-1)/D))))) W_dot(i)=I(i)*(Ron*uV/D); W(i)=W(i-1)+W_dot(i)*delta_t*(1-(2*W(i-1)/D-1)^(2*P_coeff));%%+1e-18*sign(I(i)); elseif ((win==1) && ((abs (I(i)) ) < (V_t/ (Ron*W(i-1)/D+Roff*(1-W(i-1)/D))))) W(i)=W(i-1); W_dot(i)=0; end % case this is Biolek window if ((win==2) && ((abs(I(i))) >= (V_t/(Ron*W(i-1)/D+Roff*(1-W(i-1)/D))))) W_dot(i)=I(i)*(Ron*uV/D); W(i)=W(i-1)+W_dot(i)*delta_t*(1-(W(i-1)/D-heaviside(-I(i)))^(2*P_coeff)); elseif ((win==2) && ((abs(I(i))) < (V_t/(Ron*W(i-1)/D+Roff*(1-W(i-1)/D))))) W(i)=W(i-1); W_dot(i)=0; end % case this is Prodromakis window if ((win==3) && ((abs(I(i))) >= (V_t/(Ron*W(i-1)/D+Roff*(1-W(i-1)/D))))) W_dot(i)=I(i)*(Ron*uV/D); W(i)=W(i-1)+W_dot(i)*delta_t*(J*(1-((W(i-1)/D-0.5)^2+0.75)^P_coeff)); elseif ((win==3) && ((abs(I(i))) < (V_t/(Ron*W(i-1)/D+Roff*(1-W(i-1)/D))))) W(i)=W(i-1); W_dot(i)=0; end % correct the w vector according to bounds [0 D] if W(i) < 0 W(i) = 0; W_dot(i)=0; elseif W(i) > D W(i) = D; W_dot(i)=0; end end R=Ron*W/D+Roff*(1-W/D); %%this parameter might be useful for debug V= R.*I; figure(1); plot(V(20e3:end),I(20e3:end)); title('I-V curve'); xlabel('V[volt]'); ylabel('I[amp]'); figure(2); plot(tspan_vector,W/D,'r'); title('W/D as func of time'); xlabel('time[sec]'); legend('W/D'); end %% Simmons Tunnel Barrier model if (model==1) points=2e5; tspan=[0 num_of_cycles/freq]; t = linspace(tspan(1),tspan(2),points); curr = amp*sin(freq*2*pi*t); X=zeros(1,points); X_dot=zeros(1,points); delta_t=t(2)-t(1); X(1)=w_init*D; for i=2:(length(t)) if curr(i)> 0 X_dot(i)=c_off*sinh(curr(i)/i_off)*exp(-exp((X(i-1)-a_off)/X_c-abs(curr(i))/b)-X(i-1)/X_c); else X_dot(i)=c_on*sinh(curr(i)/i_on)*exp(-exp(-(X(i-1)-a_on)/X_c-abs(curr(i))/b)-X(i-1)/X_c); end X(i)=X(i-1)+delta_t*X_dot(i); if (X(i) < 0) X(i) = 0; X_dot(i)=0; elseif (X(i) > D) X(i) = D; X_dot(i)=0; end end R=Roff.*X./D+Ron.*(1-X./D); V=R.*curr; figure(1); plot(V(20e3:end),curr(20e3:end)); title('I-V curve'); xlabel('V[volt]'); ylabel('I[amp]'); figure(2); plot(t,X/D,'r'); title('X/D as func of time'); xlabel('time[sec]'); legend('X/D'); end %% Team model if (model==2) points=2e5; tspan=[0 num_of_cycles/freq]; t = linspace(tspan(1),tspan(2),points); curr = amp*sin(freq*2*pi*t); X=zeros(1,points); X_dot=zeros(1,points); delta_t=t(2)-t(1); X(1)=w_init*D; for i=2:(length(t)) % case this is an ideal window if (win == 0) if (curr(i) > 0) && (curr(i) > i_off) X_dot(i)=k_off*(curr(i)/i_off-1)^alpha_off; X(i)=X(i-1)+delta_t*X_dot(i); elseif (curr(i) <= 0) && (curr(i) < i_on) X_dot(i)=k_on*(curr(i)/i_on-1)^alpha_on; X(i)=X(i-1)+delta_t*X_dot(i); else X(i)=X(i-1); X_dot(i)=0; end end % case this is Jogelkar window if (win == 1) if (curr(i) > 0) && (curr(i) > i_off) X_dot(i)=k_off*(curr(i)/i_off-1)^alpha_off; X(i)=X(i-1)+delta_t*X_dot(i).*(1-(2*X(i-1)/D-1)^(2*P_coeff)); elseif (curr(i) <= 0) && (curr(i) < i_on) X_dot(i)=k_on*(curr(i)/i_on-1)^alpha_on; X(i)=X(i-1)+delta_t*X_dot(i).*(1-(2*X(i-1)/D-1)^(2*P_coeff)); else X(i)=X(i-1); X_dot(i)=0; end end % case this is Biolek window if (win == 2) if (curr(i) > 0) && (curr(i) > i_off) X_dot(i)=k_off*(curr(i)/i_off-1)^alpha_off; X(i)=X(i-1)+delta_t*X_dot(i).*(1-(1-X(i-1)/D-heaviside(curr(i)))^(2*P_coeff)); elseif (curr(i) <= 0) && (curr(i) < i_on) X_dot(i)=k_on*(curr(i)/i_on-1)^alpha_on; X(i)=X(i-1)+delta_t*X_dot(i).*(1-(1-X(i-1)/D-heaviside(curr(i)))^(2*P_coeff)); else X(i)=X(i-1); X_dot(i)=0; end end % case this is Prodromakis window if (win == 3) if (curr(i) > 0) && (curr(i) > i_off) X_dot(i)=k_off*(curr(i)/i_off-1)^alpha_off; X(i)=X(i-1)+delta_t*X_dot(i).*(J*(1-((X(i-1)/D-0.5)^2+0.75)^P_coeff)); elseif (curr(i) <= 0) && (curr(i) < i_on) X_dot(i)=k_on*(curr(i)/i_on-1)^alpha_on; X(i)=X(i-1)+delta_t*X_dot(i).*(J*(1-((X(i-1)/D-0.5)^2+0.75)^P_coeff)); else X(i)=X(i-1); X_dot(i)=0; end end % case this is Kvatinsky window if (win == 4) if (curr(i) > 0) && (curr(i) > i_off) X_dot(i)=k_off*(curr(i)/i_off-1)^alpha_off*exp(-exp((X(i-1)-a_off)/X_c)); X(i)=X(i-1)+delta_t*X_dot(i); elseif (curr(i) <= 0) && (curr(i) < i_on) X_dot(i)=k_on*(curr(i)/i_on-1)^alpha_on*exp(-exp(-(X(i-1)-a_on)/X_c)); X(i)=X(i-1)+delta_t*X_dot(i); else X(i)=X(i-1); X_dot(i)=0; end end if (X(i) < 0) X(i) = 0; X_dot(i)=0; elseif (X(i) > D) X(i) = D; X_dot(i)=0; end end if (iv == 0) %case I-V relation is linear R=Roff.*X./D+Ron.*(1-X./D); V=R.*curr; else %case the I-V relation is nonlinear lambda = log(Roff/Ron); V=Ron*curr.*exp(lambda*(X-x_on)/(x_off-x_on)); end figure(1); plot(V(20e3:end),curr(20e3:end)); title('I-V curve'); xlabel('V[volt]'); ylabel('I[amp]'); figure(2); plot(t,X/D); title('X/D as func of time'); xlabel('time[sec]'); legend('X/D'); end %% Nonlinear Ion Drift model if (model==3) points=40000; tspan=[0 num_of_cycles/freq]; t = linspace(tspan(1),tspan(2),points); delta_t = t(2) - t(1); V = amp*sin(freq*2*pi*t); W=zeros(size((t))); W_dot=zeros(size((t))); curr=zeros(size((t))); W(1) = w_init; for i=2:points % case this is an ideal window if ((win==0) || (win==4)) W_dot(i)=a*V(i)^q; W(i)=W(i-1)+W_dot(i)*delta_t; end % case this is Jogelkar window if (win==1) W_dot(i)=a*V(i)^q; W(i)=W(i-1)+W_dot(i)*delta_t*(1-(2*W(i-1)-1)^(2*P_coeff)); end % case this is Biolek window if (win==2) W_dot(i)=a*V(i)^q; W(i)=W(i-1)+W_dot(i)*delta_t*(1-(W(i-1)-heaviside(-V(i-1)))^(2*P_coeff)); end % case this is Prodromakis window if (win==3) W_dot(i)=a*V(i)^q; W(i)=W(i-1)+W_dot(i)*delta_t*(J*(1-((W(i-1)-0.5)^2+0.75)^P_coeff)); end % correct the w vector according to bounds [0 D] if W(i) < 0 W(i) = 0; W_dot(i)=0; elseif W(i) > 1 W(i) = 1; W_dot(i)=0; end curr(i)=W(i)^n*beta*sinh(alpha*V(i))+c*(exp(g*V(i))-1); end figure(1); plot(V(20e3:end),curr(20e3:end)); title('I-V curve'); xlabel('V[volt]'); ylabel('I[amp]'); figure(2); plot(t,W,'r'); title('W/D'); xlabel('time[sec]'); legend('W/D'); end %% VTEAM model if (model==4) points=2e5; tspan=[0 num_of_cycles/freq]; t = linspace(tspan(1),tspan(2),points); v = amp*sin(freq*2*pi*t); X=zeros(1,points); X_dot=zeros(1,points); delta_t=t(2)-t(1); X(1)=w_init*D; for i=2:(length(t)) % case this is an ideal window if (win == 0) if (v(i) > 0) && (v(i) > v_off) X_dot(i)=k_off*(v(i)/v_off-1)^alpha_off; X(i)=X(i-1)+delta_t*X_dot(i); elseif (v(i) <= 0) && (v(i) < v_on) X_dot(i)=k_on*(v(i)/v_on-1)^alpha_on; X(i)=X(i-1)+delta_t*X_dot(i); else X(i)=X(i-1); X_dot(i)=0; end end % case this is Jogelkar window if (win == 1) if (v(i) > 0) && (v(i) > v_off) X_dot(i)=k_off*(v(i)/v_off-1)^alpha_off; X(i)=X(i-1)+delta_t*X_dot(i).*(1-(2*X(i-1)/D-1)^(2*P_coeff)); elseif (v(i) <= 0) && (v(i) < v_on) X_dot(i)=k_on*(v(i)/v_on-1)^alpha_on; X(i)=X(i-1)+delta_t*X_dot(i).*(1-(2*X(i-1)/D-1)^(2*P_coeff)); else X(i)=X(i-1); X_dot(i)=0; end end % case this is Biolek window if (win == 2) if (v(i) > 0) && (v(i) > v_off) X_dot(i)=k_off*(v(i)/v_off-1)^alpha_off; X(i)=X(i-1)+delta_t*X_dot(i).*(1-(1-X(i-1)/D-heaviside(v(i)))^(2*P_coeff)); elseif (v(i) <= 0) && (v(i) < v_on) X_dot(i)=k_on*(v(i)/v_on-1)^alpha_on; X(i)=X(i-1)+delta_t*X_dot(i).*(1-(1-X(i-1)/D-heaviside(v(i)))^(2*P_coeff)); else X(i)=X(i-1); X_dot(i)=0; end end % case this is Prodromakis window if (win == 3) if (v(i) > 0) && (v(i) > v_off) X_dot(i)=k_off*(v(i)/v_off-1)^alpha_off; X(i)=X(i-1)+delta_t*X_dot(i).*(J*(1-((X(i-1)/D-0.5)^2+0.75)^P_coeff)); elseif (v(i) <= 0) && (v(i) < v_on) X_dot(i)=k_on*(v(i)/v_on-1)^alpha_on; X(i)=X(i-1)+delta_t*X_dot(i).*(J*(1-((X(i-1)/D-0.5)^2+0.75)^P_coeff)); else X(i)=X(i-1); X_dot(i)=0; end end % case this is Kvatinsky window if (win == 4) if (v(i) > 0) && (v(i) > v_off) X_dot(i)=k_off*(v(i)/v_off-1)^alpha_off*exp(-exp((X(i-1)-a_off)/X_c)); X(i)=X(i-1)+delta_t*X_dot(i); elseif (v(i) <= 0) && (v(i) < v_on) X_dot(i)=k_on*(v(i)/v_on-1)^alpha_on*exp(-exp(-(X(i-1)-a_on)/X_c)); X(i)=X(i-1)+delta_t*X_dot(i); else X(i)=X(i-1); X_dot(i)=0; end end if (X(i) < 0) X(i) = 0; X_dot(i)=0; elseif (X(i) > D) X(i) = D; X_dot(i)=0; end end if (iv == 0) %case I-V relation is linear R=Roff.*X./D+Ron.*(1-X./D); curr=v./R; else %case the I-V relation is nonlinear lambda = log(Roff/Ron); curr=v./(Ron.*exp(lambda*(X-x_on)/(x_off-x_on))); end figure(1); plot(v(20e3:end),curr(20e3:end)); title('I-V curve'); xlabel('V[volt]'); ylabel('I[amp]'); figure(2); plot(t,X/D); title('X/D as func of time'); xlabel('time[sec]'); legend('X/D'); end