// Author of the Code : Mohit D. Ganeriwala, E. G. Marin, F. Pasadas, F. G. Ruiz and Andres Godoy - University of Granada, Spain
// In case of using the code please cite the following paper 
// Cruces, S., Ganeriwala, M.D., Lee, J., Völkel, L., Braun, D., Grundmann, A., Ran, K., González Marín, E., Kalisch, H., Heuken, M., Vescan, A., Mayer, J., Godoy, A., Daus, A. and Lemme, M.C. (2025), 
// Volatile MoS2 Memristors with Lateral Silver Ion Migration for Artificial Neuron Applications. Small Sci., 5: 2400523. https://doi.org/10.1002/smsc.202400523


`include "disciplines.vams"
`include "constants.vams"
module Diff_Mem(l, r);
	inout l, r;
	electrical l, r, n_d;
	parameter real g0 = 1.524 from [0:inf);
	parameter real alpha0 = 1.941 from [0:inf);  //1/nm
	parameter real V0 = 84.33 from (0:inf);    //1/Vnm   //epsilon_s/q
	parameter real Vel0 = 18.84 from (0:inf);
	parameter real I0 = 14.79 from (0:inf);  // Ampere.nm^3/V  //q*u_0*A
	parameter real Beta = 0.2547 from (0:inf);
	parameter real gamma0 = 11.07 from (0:inf);
	parameter real Ea = 0.405 from (0:inf);  //eV
	parameter real a0 = 8.690 from (0:inf); // nm
	parameter real L = 1.393e3 from (0:inf);  // nm
	parameter real shift = 2.472 from [0:inf);
	parameter real Nt = 369.8e-9 from[0:inf);     // 1/nm3
	parameter real Ir = 36.56e-6 from(0:inf); // Ampere

	parameter real max_d = 1.849 from (0:inf);
	parameter real min_d = 0.2547 from (0:inf);

	parameter real smoothing = 1e-8 from (0:inf);
	parameter real GMIN = 1e-12 from (0:inf);
	parameter real Kclip = 50 from (0:inf);

	real d, ddt_d, Gamma, Fw1, Fw2, clip_max_d, clip_min_d, T1;

	analog function real smoothstep;
		input x, smoothing;
		real x, smoothing;
		begin
			smoothstep = 0.5*(x/sqrt(x*x + smoothing)+1);
		end
	endfunction // smoothstep
	
	analog begin
		d = V(n_d, r);
		T1 = ((V0*V(l,r)*exp(-alpha0*d)/pow(L,2))+Nt)/(1+exp(-alpha0*d));
		I(l,r) <+ I0*T1*exp(g0*sqrt(V(l,r)/L)/$vt)*(V(l,r)/L) +  Ir*pow(V(l,r),2)/pow(L,3)  +  GMIN*V(l, r);
		
		Gamma = gamma0 - Beta * pow(d, 3);
		ddt_d = -Vel0*exp(-Ea/$vt)*sinh((V(l, r)-shift)*Gamma*a0/L/$vt);

		Fw1 = smoothstep(min_d-d, smoothing);
		Fw2 = smoothstep(d-max_d, smoothing);
		clip_min_d = (limexp(Kclip*(min_d-d)) - ddt_d) * Fw1;
		clip_max_d = (-limexp(Kclip*(d-max_d)) - ddt_d) * Fw2;

		I(n_d, r) <+ ddt_d + clip_min_d + clip_max_d;
		I(n_d, r) <+ ddt(-1e-9*d);
		
	end
endmodule