////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//Copyright @ 2013 Massachusetts Institute of Technology (MIT)
 
//The terms under which the software and associated documentation (the Software) is provided are as the following:
 
//The Software is provided "as is", without warranty of any kind, express or implied, including but not limited to the warranties of merchantability, fitness for a particular purpose and noninfringement. In no event shall the authors or copyright holders be liable for any claim, damages or other liability, whether in an action of contract, tort or otherwise, arising from, out of or in connection with the Software or the use or other dealings in the Software.
 
//MIT grants, free of charge, to any users the right to modify, copy, and redistribute the Software, both within the user's organization and externally, subject to the following restrictions:
 
//1. The users agree not to charge for the MIT code itself but may charge for additions, extensions, or support.
 
//2. In any product based on the Software, the users agree to acknowledge the MIT VS Model Research Group that developed the software. This acknowledgment shall appear in the product documentation.
 
//3. The users agree to obey all U.S. Government restrictions governing redistribution or export of the software.
 
//4. The users agree to reproduce any copyright notice which appears on the software on any copy or modification of such made available to others.
 
//Agreed to by 
//Dimitri A. Antoniadis, MIT
//May 27 2013
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// VerilogA for ambipolar virtual-source (AVS) based self-consistent transport/capacitance model for graphene FETs
// transport and charge model:  S. Rakheja, et al T-NANO 2014.
// Implemented on April 23, 2014 by S. Rakheja


`include "constants.vams"
`include "disciplines.vams"

module avs_1_0_0(d, g, s);
inout d, g, s;
electrical d, g, s;
electrical di, si;

// AVS parameters
parameter real version = 1.00;
parameter real W = 20e-6 from (0:inf); // Device width [m]
parameter real Lg = 650e-9 from (0:inf); // Device length [m]
parameter real Cg = 3.63e-3 from (0:inf); // Gate capacitance [F/m^2]
parameter real n0 = 2 from (1:inf); // Non-ideality factor [unit-less]
parameter real alpha = 6 from (1:inf); // Shift in threshold voltage in subthreshold and strong inversion regimes [unit-less]
parameter real zeta = 0.2 from (0:1]; // Channel ballisticity parameter. invoked in dynamic model.
parameter real CTM_select = 1 from [1:inf); //Parameter to select charge-transport model. If CTM_select = 1, then classical DD-NVSAT model is used, else blended QB and DD-NVSAT model is used.
parameter real Cofs = 0 from [0:inf); // Outer-fringing capacitance for the source-terminal [F/m]
parameter real Cofd = 0 from [0:inf); // Outer-fringing capacitance for the source-terminal [F/m]
parameter real Rc_elec = 100e-6 from (0:inf); // Electron-branch channel-access resistance [Ohm-meter]
parameter real Rc_hole = 100e-6 from (0:inf); // Hole-branch channel-access resistance [Ohm-meter]
parameter real beta = 1.8 from (0:inf); // Empirical parameter for Fsat
parameter real Qmin = 0 from [0:inf); // Minimum background doping [C/m^2]
parameter real mu = 600e-4 from (0:inf); // Carrier mobility [m^2/Vs]
parameter real vxo = 3e5 from (0:inf); // VS velocity [m/s]
parameter real delta = 0 from [0:inf); // Shift in threshold voltage for charge trapping [V].
parameter real Vmin0 = 0 from (-inf:inf); // Minimum conduction point/Dirac voltage at Vds = 0V. [V]
parameter real Tjun = 298 from [173:inf); // Junction temperature [K] 
parameter real mc = 0.06 from [0.005:10]; // Relative effective mass of electrons 

`define  SMALL_VALUE  (1e-10)
`define  LARGE_VALUE  (40) 

real Relec, Rhole, phit0, me, Qref, Qelec, Qhole, etae, etah;
real Vgsraw, Vgdraw, Vds, Vdsi, Vgs, Vgsi, Vdg, Vdgi, Vtn, Vtp, Vdsat,Vdint, Vsint, dir; 
real FFSe, FFSh,Id_elec, Id_hole, Id, Fsat, Vdratio, Vdbetabeta, Vdbeta; 
real Vdsatqe, Vdsatqh, Fsatqe, Fsatqh, te, th, qs_nvsat, qd_nvsat;

real kq, kq2,lhs, rhs, qse_qb, qsh_qb, qde_qb, qdh_qb, qs_net_qb, qd_net_qb, lhse, lhsh, rhse, rhsh; 
real Fsatc,m0, qx, qy, Qs_ch, Qd_ch, Qofx, Qofy, Qs, Qd, Qg;

analog begin
	// Voltage definitions
	Vgsraw = V(g)-V(si);
	Vgdraw = V(g)-V(di);
	if (Vgsraw >= Vgdraw) begin
		Vds = V(d)-V(s);
		Vgs = V(g)-V(s);
		Vdg = V(d)-V(g);
		Vdsi = V(di)-V(si);
		Vdgi = V(di)-V(g);
		Vgsi = Vgsraw;
		dir = 1;
	end
	else begin
		Vds = V(s)-V(d);
		Vgs = V(g)-V(d);
		Vdg = V(s)-V(g);
		Vdsi = V(si)-V(di);
		Vdgi = V(si)-V(g);	
		Vgsi = Vgdraw;
		dir = -1;
	end


	Relec = Rc_elec/W; // Electron-branch resistance [Ohms]
	Rhole = Rc_hole/W; // Hole-branch resistance [Ohms]
	phit0 = $vt(Tjun); // Thermal voltage [V]
	m0 = 9.1e-31; // Mass of free electron [Kg]
	me = m0*mc; // Effective mass of electron [Kg]

	Vtn = Vmin0+delta; // Threshold voltage for electron branch [V]
	Vtp = Vmin0-delta; // Threshold voltage for hole branch [V]

	Qref = Cg*n0*phit0;
	
	FFSe = 1/(1+exp((Vgsi-(Vtn-alpha*phit0/2.0))/(alpha*phit0)));// Fermi-like function for electrons
	FFSh = 1/(1+exp((Vdgi+(Vtp+alpha*phit0/2.0))/(alpha*phit0))); // Fermi-like function for holes

	etae = (Vgsi-(Vtn-alpha*phit0*FFSe))/(n0*phit0);
	etah = (Vdgi+(Vtp+alpha*phit0*FFSh))/(n0*phit0);

	if (etae <= `LARGE_VALUE) begin
		Qelec = Qref*ln(1.0+exp(etae))+Qmin; // Net electron concentration at the VS for electrons [C/m^2]
	end
	else begin
		Qelec = Qref*etae+Qmin; 
	end

	if (etah <= `LARGE_VALUE) begin
		Qhole = Qref*ln(1.0+exp(etah))+Qmin; // Net electron concentration at the VS for holes [C/m^2]
	end
	else begin
		Qhole = Qref*etah+Qmin; 
	end

	Vdsat = vxo*Lg/mu; // Saturation voltage for transport
	//Fsat = Vdsi/Vdsat/(pow(1.0+pow(Vdsi/Vdsat,beta), 1.0/beta); 
	Vdratio	= abs( Vdsi/ Vdsat);
	Vdbeta	= pow( Vdratio, beta);
	Vdbetabeta = pow( 1.0 + Vdbeta, 1.0/ beta);

	Fsat = Vdratio / Vdbetabeta;		
	Id_elec = Qelec*Fsat*vxo*W; // Electron-branch current
	Id_hole = Qhole*Fsat*vxo*W; // Hole-branch current


	Id = Id_elec+Id_hole;

// Begin charge model

//quasi-ballistic charge computation

	kq = sqrt(2.0*zeta*`P_Q/me*(Vdsi))/vxo;
	kq2=kq*kq;
    	lhs = 2.0/(3.0*kq2*kq2)*((2.0+2.0*kq2)*pow((1+kq2),0.5)-(2.0+3.0*kq2));
    	rhs = 2.0/(3.0*kq2*kq2)*(2.0+(kq2-2.0)*pow((kq2+1),0.5));

	qse_qb= -(Qelec*lhs)*W*Lg;
    	qde_qb = -(Qelec*rhs)*W*Lg;
    	qsh_qb = (Qhole*rhs*W*Lg);
    	qdh_qb = (Qhole*lhs*W*Lg);

	if (Vds < 1e-3) begin
		qs_net_qb = 0.5*(Qhole-Qelec)*W*Lg;
    		qd_net_qb = 0.5*(Qhole-Qelec)*W*Lg;
	end
	else begin
		qs_net_qb = (qse_qb+qsh_qb);
		qd_net_qb = (qde_qb+qdh_qb);
	end

// DD-NVSAT model begin (Check Tsividis, chapter #8)
	Vdsatqe = abs(Vgsi-Vtn);
	Fsatqe = Vdsi/Vdsatqe/(pow(1.0+pow(Vdsi/Vdsatqe, beta),1.0/beta));
	te = 1-Fsatqe;
	
	Vdsatqh = abs(Vdgi+Vtp);
	Fsatqh = Vdsi/Vdsatqh/(pow(1.0+pow(Vdsi/Vdsatqh, beta),1.0/beta));
	th = 1-Fsatqh;

	lhse = (6 + 12 * te + 8 * te * te + 4 * te * te * te)/(15.0*(1+te)*(1+te));
	lhsh = (4 + 8 * th + 12 * th * th + 6 * th * th * th)/(15.0*(1+th)*(1+th));
	rhse = (4 + 8 * te + 12 * te * th + 6 * te * te * te)/(15.0*(1+te)*(1+te));
	rhsh = (6 + 12 * th + 8 * th * th + 4 * th * th * th)/(15.0*(1+th)*(1+th));

	qs_nvsat = -Qelec*W*Lg*lhse+Qhole*W*Lg*lhsh;
	qd_nvsat = -Qelec*W*Lg*rhse+Qhole*W*Lg*rhsh;

	if (CTM_select==1) begin
        	Fsatc =0;
	end
    	else begin
        	Fsatc = 1;
	end
    
    	qx = qs_nvsat*(1-Fsatc)+qs_net_qb*(Fsatc);
    	qy = qd_nvsat*(1-Fsatc)+qd_net_qb*(Fsatc);

    	Qs_ch = ((1+dir)*qx+(1-dir)*qy)/2.0;
    	Qd_ch = ((1-dir)*qx+(1+dir)*qy)/2.0;

	Vsint=V(s)+dir*(Id_elec*Relec+Id_hole*Rhole);
    	Vdint=V(d)-dir*(Id_elec*Relec+Id_hole*Rhole);

	Qofx=-Cofs*(Vgsraw)*W;
	Qofy=-Cofd*(Vgdraw)*W;

	Qs = Qs_ch + Qofx;
	Qd = Qd_ch + Qofy;
	Qg = -(Qs+Qd);

	//Sub-circuit initialization
    	I(di,si) <+ dir * Id;
	V(d,di) <+ dir*(Id_elec*Relec+Id_hole*Rhole);
	V(si,s) <+ dir*(Id_elec*Relec+Id_hole*Rhole);

	I(si,g)<+ ddt(Qs);
	I(di,g)<+ ddt(Qd);

end
endmodule