////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//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 virtual-source (VS) based self-consistent transport/capacitance model for Si MOSFET
// transport model:  A. Khakifirooz, et al, p. 1674, T-ED 2009.
// charge model: L. Wei et al, p. 1263, T-ED 2012.
// Implemented on July 15, 2013 by S. Rakheja
// Modified on Sep. 19, 2013 by S. Rakheja
// Modified on July 17, 2015 by S. Rakheja 
// List of changes: 
// 1. Energy transfer parameter zeta added
// 2. Qinv definition corrected to use FF instead of FFpre
// 3. DIBL correction normalized
// 4. Ballistic charge calculated for linear potential profile (default) 
// 5. Inner fringing charge corrected 
// 6. Unused variables removed

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

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

// Original VS parameters
parameter real 		version 	= 1.10;                                         		//      MVS model version = 1.0.1
parameter integer 	type 		= 1 			from [-1 : 1] exclude 0;  		//	type of transistor. nFET type=1; pFET type=-1
parameter real 		W          	= 1e-4 			from (0:inf);	        		//	Transistor width [cm]
parameter real 		Lgdr 		= 80e-7 		from (0:inf);                   	//	Physical gate length [cm]. //	This is the designed gate length for litho printing.
parameter real 		dLg 		= 10.5e-7 		from (0:inf);                   	//	Overlap length including both source and drain sides [cm]  
parameter real 		Cg          	= 2.2e-6 		from (0:inf);                   	//	Gate-to-channel areal capacitance at the virtual source [F/cm^2]
parameter real		etov 		= 1.3e-3 		from (0:inf);                   	//	Equivalent thickness of dielectric at S/D-G overlap [cm]
parameter real 		delta 		= 0.10 			from [0:inf);                   	//	Drain-induced-barrier-lowering (DIBL) [V/V]
parameter real 		n0         	= 1.5 			from [0:inf);                   	// 	Subthreshold swing factor [unit-less] {typically between 1.0 and 2.0}
parameter real 		Rs0 		= 100 			from (0:inf);                   	// 	Access resistance on s-terminal [Ohms-micron]
parameter real 		Rd0 		= 100 			from (0:inf);                   	// 	Access resistance on d-terminal [Ohms-micron] 
                                                                                			//	Generally, Rs0 = Rd0 for symmetric source and drain
parameter real 		Cif 		= 1e-12 		from [0:inf);                   	//	Inner fringing S or D capacitance [F/cm] 
parameter real 		Cof 		= 2e-13 		from [0:inf);                   	//	Outer fringing S or D capacitance [F/cm] 
parameter real 		vxo 		= 0.765e7 		from (0:inf);                   	//	Virtual source injection velocity [cm/s]
parameter real 		mu         	= 200 			from (0:inf);                   	//	Low-field mobility [cm^2/V.s]
parameter real 		beta 		= 1.7 			from (0:inf);                   	//	Saturation factor. Typ. nFET=1.8, pFET=1.6
parameter real      	Tjun       	= 298           	from [173:inf);                 	//  	Junction temperature [K]
parameter real 		phib 		= 1.2;                                          		//	~abs(2*phif)>0 [V]
parameter real 		gamma 		= 0.0 			from [0:inf);                   	//	Body factor  [sqrt(V)]
parameter real 		Vt0 		= 0.486;                                        		//	Strong inversion threshold voltage [V] 
parameter real 		alpha 		= 3.5;                                          		//	Empirical parameter for threshold voltage shift between strong and weak inversion.
parameter real 		mc          	= 0.2 			from [0.01 : 10];               	// 	Choose an appropriate value between 0.01 to 10 
                                                                                			//	For, values outside of this range,convergence or accuracy of results is not guaranteed
parameter integer	CTM_select 	= 1             	from [1 : inf);                		// 	If CTM_select = 1, then classic DD-NVSAT model is used
                                                                                			//	For CTM_select other than 1,blended DD-NVSAT and ballistic charge transport model is used 
parameter real 		CC 		= 0                 	from [0:inf);                   	//	Fitting parameter to adjust Vg-dependent inner fringe capacitances(Not used in this version)
parameter real 		nd 		= 0                 	from [0:inf);                   	// 	Punch-through factor [1/V]
parameter real		zeta		= 1.0			from (0:1];				// 	Energy-transfer factor

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

real Rs, Rd, Vds, Vgs, Vgsraw, Vgdraw, Vbs, Vdsi, Vgsi, Vbsi, dir;
real Leff, me, phit;
real n, nphit, aphit, Vtpcorr, eVgpre, FFpre, ab, Vcorr, Vgscorr, Vbscorr, Vt0bs, Vt0bs0, Vtp, Vtp0;
real eVg, FF, eVg0, FF0, Qref, eta, eta0;
real Qinv, Qinv_corr, vx0, Vdsats, Vdsat,Vdratio, Vdbeta, Vdbetabeta, Fsat, Id ;
real Vgt, psis, Vgta, Vdsatq, Fsatq, x, den;
real qsc, qdc, qi, kq, kq2, kq4, tol, qsb, qdb, qs, qd, Qs, Qd;
real Qb, etai, Qinvi, dQinv, dibl_corr; 
real Qinvs, Qinvd, Qsov, Qdov, Vt0x, Vt0y, Fs_arg, Fs, Fd_arg, Fd, FFx, FFy, Qsif, Qdif, Qg, a, Cofs, Cofd;
real A, B;

analog begin

	//Voltage definitions
	Vgsraw			=	type * ( V(g) - V(si) );
	Vgdraw			=	type * ( V(g) - V(di) );
	if (Vgsraw >= Vgdraw) begin
		Vds		=	type * ( V(d) - V(s) ); 
		Vgs		=	type * ( V(g) - V(s) );
		Vbs		=	type * ( V(b) - V(s) );
		Vdsi		=	type * ( V(di) - V(si) );
		Vgsi		=	Vgsraw;
		Vbsi		=	type * ( V(b) - V(si) );
		dir		=	1;
	end
	else begin
		Vds		=	type * ( V(s) - V(d) );
		Vgs		=	type * ( V(g) - V(d) );
		Vbs		=	type * ( V(b) - V(d) );
		Vdsi		=	type * ( V(si) - V(di) );
		Vgsi		=	Vgdraw;
		Vbsi		=	type * ( V(b) - V(di) );
		dir		=	-1;
	end


    //Parasitic element definition
	Rs			=	1e-4/ W * Rs0;                                          	// 	s-terminal resistance [ohms]
	Rd			=	Rs;                                                     	// 	d-terminal resistance [ohms] For symmetric source and drain Rd = Rs. 
      //Rd			=	1e-4/ W * Rd0;                                          	// 	d-terminal resistance [ohms] {Uncomment for asymmetric source and drain resistance.}
    	Cofs			= 	( 0.345e-12/ etov ) * dLg/ 2.0 + Cof; 				// 	s-terminal outer fringing cap [F/cm]
   	Cofd 			= 	( 0.345e-12/ etov ) * dLg/ 2.0 + Cof; 				// 	d-terminal outer fringing cap [F/cm]
	Leff 			= 	Lgdr - dLg;                                         		// 	Effective channel length [cm]. After subtracting overlap lengths on s and d side 
	
    	phit        		=   	$vt(Tjun);                                             		// 	Thermal voltage, kT/q [V]                
	me			=	(9.1e-31) * mc;                                        		// 	Carrier mass [Kg]
	n			=	n0 + nd * Vdsi;                                          	// 	Total subthreshold swing factor taking punchthrough into account [unit-less]
	nphit			= 	n * phit;                                           		// 	Product of n and phit [used as one variable]
	aphit			=	alpha * phit;                                       		// 	Product of alpha and phit [used as one variable]

	
	//Correct Vgsi and Vbsi
	//Vcorr is computed using external Vbs and Vgs but internal Vdsi, Qinv and Qinv_corr are computed with uncorrected Vgs, Vbs and corrected Vgs, Vbs respectively.	
	Vtpcorr			=	Vt0 + gamma * (sqrt(abs(phib - Vbs))- sqrt(phib))- Vdsi * delta;// 	Calculated from extrinsic Vbs
	eVgpre 			= 	exp(( Vgs - Vtpcorr )/ ( aphit * 1.5 )); 			// 	Calculated from extrinsic Vgs
    	FFpre 			= 	1.0/ ( 1.0 + eVgpre );						// 	Only used to compute the correction factor
   	ab			=	2 * ( 1 - 0.99 * FFpre ) * phit;  
    	Vcorr			=	( 1.0 + 2.0 * delta ) * ( ab/ 2.0 ) * ( exp( -Vdsi/ ab )); 	// 	Correction to intrinsic Vgs
	Vgscorr			=	Vgsi + Vcorr; 							// 	Intrinsic Vgs corrected (to be used for charge and current computation)
	Vbscorr			=	Vbsi + Vcorr; 							// 	Intrinsic Vgs corrected (to be used for charge and current computation)
	Vt0bs			=	Vt0 + gamma * (sqrt( abs( phib - Vbscorr)) - sqrt( phib )); 	// 	Computed from corrected intrinsic Vbs
	Vt0bs0			=	Vt0 + gamma * (sqrt( abs( phib - Vbsi)) - sqrt( phib )); 	// 	Computed from uncorrected intrinsic Vbs
	Vtp			=	Vt0bs - Vdsi * delta - 0.5 * aphit; 				// 	Computed from corrected intrinsic Vbs and intrinsic Vds
   	Vtp0			=	Vt0bs0 - Vdsi * delta - 0.5 * aphit; 				// 	Computed from uncorrected intrinsic Vbs and intrinsic Vds
	eVg			=	exp(( Vgscorr - Vtp )/ ( aphit )); 				// 	Compute eVg factor from corrected intrinsic Vgs
	FF			=	1.0/ ( 1.0 + eVg );
	eVg0			=	exp(( Vgsi - Vtp0 )/ ( aphit )); 				// 	Compute eVg factor from uncorrected intrinsic Vgs
	FF0			=	1.0/ ( 1.0 + eVg0 );
	Qref			=	Cg * nphit;	
	eta			=	( Vgscorr - ( Vt0bs - Vdsi * delta - FF * aphit ))/ ( nphit ); 	// 	Compute eta factor from corrected intrinsic Vgs and intrinsic Vds
	eta0			=	( Vgsi - ( Vt0bs0 - Vdsi * delta - FF * aphit ))/ ( nphit );	//	Compute eta0 factor from uncorrected intrinsic Vgs and internal Vds. 
													//	Using FF instead of FF0 in eta0 gives smoother capacitances.

	//Charge at VS in saturation (Qinv)
    	if (eta  <= `LARGE_VALUE) begin
		Qinv_corr  	= 	Qref * ln( 1.0 + exp(eta) );
	end
	else begin
		Qinv_corr	=	Qref * eta;
	end	
	if (eta0 <= `LARGE_VALUE) begin 
		Qinv		=	Qref * ln( 1.0 + exp(eta0) );					// 	Compute charge w/ uncorrected intrinsic Vgs for use later on in charge partitioning
	end
	else	begin
		Qinv		=	Qref * eta0;
	end


	//Transport equations
	vx0			=	vxo;	
	Vdsats			=	vx0 * Leff/ mu;							
	Vdsat			=	Vdsats * ( 1.0 - FF ) + phit * FF; 				// 	Saturation drain voltage for current
	Vdratio			=	abs( Vdsi/ Vdsat);
	Vdbeta			=	pow( Vdratio, beta);
	Vdbetabeta 		= 	pow( 1.0 + Vdbeta, 1.0/ beta);
	Fsat			=	Vdratio / Vdbetabeta;						// 	Transition function from linear to saturation. 
													//	Fsat = 1 when Vds>>Vdsat; Fsat= Vds when Vds<<Vdsat

	//Total drain current 										
	Id			=	Qinv_corr * vx0 * Fsat * W;		


	//Calculation of intrinsic charge partitioning factors (qs and qd)
	Vgt			=	Qinv/ Cg; 							// 	Use charge computed from uncorrected intrinsic Vgs

	// Approximate solution for psis is weak inversion
	if (gamma == 0) begin
		a 		= 	1.0;
		if (eta0 <= `LARGE_VALUE) begin
            		psis    = 	phib + phit * ( 1.0 + ln( ln( 1.0 + `SMALL_VALUE + exp( eta0 ))));
		end
		else begin
            		psis    =	phib + phit * ( 1.0 + ln( eta0 ));
        	end
	end
	else begin
        	if (eta0 <= `LARGE_VALUE) begin
           		psis	=	phib + ( 1.0 - gamma )/ ( 1.0 + gamma ) * phit * ( 1.0 + ln( ln( 1.0 + `SMALL_VALUE + exp( eta0 ))));
        	end
       		else begin	
            		psis	=	phib + ( 1.0 - gamma )/ ( 1.0 + gamma ) * phit * ( 1.0 + ln( eta0 ));
        	end
        	a		=	1.0 + gamma/ ( 2.0 * sqrt( abs( psis - ( Vbsi ))));
    	end
	Vgta			=	Vgt/ a; 							// 	Vdsat in strong inversion
	Vdsatq			=	sqrt( FF0 * aphit * aphit + Vgta * Vgta);			// 	Vdsat approx. to extend to weak inversion; 
													// 	The multiplier of phit has strong effect on Cgd discontinuity at Vd=0.

	// Modified Fsat for calculation of charge partitioning
	//DD-NVSAT charge
	Fsatq			=	abs( Vdsi/ Vdsatq )/ ( pow( 1.0 + pow( abs( Vdsi/ Vdsatq ), beta ), 1.0/ beta ));
	x			=	1.0 - Fsatq;
	

	 // MVS 1.0.1 release uses Tsividis model for implementing DD-NVSAT charges. 
        //  In MVS 1.1.0, we use the formulation from L. Wei's 2012 TED paper referenced in the documentation. 
        // Both models are identical. 

        // *** From Tsividis (DD-NVSAT)
	// den			=	15 * ( 1 + x ) * ( 1 + x );
	// qsc 			=	Qinv *(6 + 12 * x + 8 * x * x + 4 * x * x * x)/ den;
	// qdc			=	Qinv *(4 + 8 * x + 12 * x * x + 6 * x * x * x)/ den;
	// qi			=	qsc + qdc; 							// 	Charge in the channel 

        // *** From L. Wei (DD-NVSAT)
        A = pow((1-x),2)/(12.0*(1-(1-x)/2.0));
	B = (5.0-2.0*(1-x))/(10.0-(1-x)/2.0);
	qsc = Qinv*(0.5-(1-x)/6.0+A*(1-B));
	qdc = Qinv*(0.5-(1-x)/3.0+A*B);
	qi = qsc+qdc;

		

	//QB charge	
	//kq 			= 	0.0;
	tol			= 	( `SMALL_VALUE * vxo/ 100.0 ) * ( `SMALL_VALUE * vxo/ 100.0 ) * me/ ( 2 * `P_Q );
	if (Vdsi <= tol) begin
		kq2		=	( 2.0 * `P_Q/ me * Vdsi * zeta)/ ( vx0 * vx0 ) * 10000.0;
		kq4		=	kq2 * kq2;
		
		// Uncomment following lines for linear potential profile
		qsb		= 	Qinv * (0.5 - kq2/12.0 +kq4/32.0);
		qdb		=	Qinv * (0.5 - kq2/6.0 +3.0*kq4/32.0);
		
		// Uncomment following lines for parabolic potential profile
		//qsb		=	Qinv * ( 0.5 - kq2/ 24.0 + kq4/ 80.0 );
		//qdb		=	Qinv * ( 0.5 - 0.125 * kq2 + kq4/ 16.0 );
	end
	else begin	
		kq		=	sqrt( 2.0 * `P_Q/ me * Vdsi * zeta )/ vx0 * 100.0;
		kq2		=	kq * kq;
		
		// Uncomment following lines for linear potential profile
		qsb		=	Qinv * (4.0*(kq2+1.0)*sqrt(kq2+1.0)-(6.0*kq2+4.0))/(3.0*kq2*kq2); 
		qdb		=	Qinv * (2.0*(kq2-2)*sqrt(kq2+1.0)+4.0)/(3.0*kq2*kq2);

		// Uncomment following lines for parabolic potential profile
		// qsb		=	Qinv * ( asinh( kq )/ kq - ( sqrt( kq2 + 1.0 ) - 1.0 )/ kq2);
		// qdb		=	Qinv * (( sqrt( kq2 + 1.0 )- 1.0 )/ kq2);
	end
	

	// Flag for classic or ballistic charge partitioning:
	if (CTM_select == 1) begin					 				// 	Ballistic blended with classic DD-NVSAT
		qs		=	qsc;								// 	Calculation of "ballistic" channel charge partitioning factors, qsb and qdb.
		qd		=	qdc;								// 	Here it is assumed that the potential increases parabolically from the
	end												// 	virtual source point, where Qinv_corr is known to Vds-dvd at the drain.
	else begin 											// 	Hence carrier velocity increases linearly by kq (below) depending on the
		qs		=	qsc * ( 1 - Fsatq * Fsatq ) + qsb * Fsatq * Fsatq;		// 	efecive ballistic mass of the carriers.
		qd		=	qdc * ( 1 - Fsatq * Fsatq ) + qdb * Fsatq * Fsatq;				
	end												
													
								
	//Body charge based on approximate surface potential (psis) calculation with delta=0 using psis=phib in Qb gives continuous Cgs, Cgd, Cdd in SI, while Cdd is smooth anyway.
	Qb 			=	-type * W * Leff * ( Cg * gamma * sqrt( abs( psis - Vbsi )) + ( a - 1.0 )/ ( 1.0 * a ) * Qinv * ( 1.0 - qi ));


	//DIBL effect on drain charge calculation.
	//Calculate dQinv at virtual source due to DIBL only. Then:Correct the qd factor to reflect this channel charge change due to Vd
	//Vt0bs0 and FF=FF0 causes least discontinuity in Cgs and Cgd but produces a spike in Cdd at Vds=0 (in weak inversion.  But bad in strong inversion)
	etai 			=	( Vgsi - ( Vt0bs0 - FF * aphit ))/ ( nphit );
	if (etai <= `LARGE_VALUE) begin
		Qinvi		=	Qref * ln( 1.0 + exp( etai ));
	end
	else begin
		Qinvi		=	Qref * etai;
	end
	dQinv			=	Qinv - Qinvi;
	dibl_corr		=	( 1.0 - FF0 ) * ( 1.0 - Fsatq ) * qi * dQinv/Qinv;
	qd			=	qd - dibl_corr;

		 
	//Inversion charge partitioning to terminals s and d
	Qinvs			=	type * Leff * (( 1 + dir ) * qs + ( 1 - dir ) * qd)/ 2.0;
	Qinvd			=	type * Leff * (( 1 - dir ) * qs + ( 1 + dir ) * qd)/ 2.0;


	//Outer fringing capacitance
	Qsov			=	Cofs * ( V(g) - V(si) );
	Qdov			=	Cofd * ( V(g) - V(di) );


	//Inner fringing capacitance
	Vt0x			=	Vt0 + gamma * ( sqrt( abs( phib - type * ( V(b) - V(si) ))) - sqrt(phib));
	Vt0y			=	Vt0 + gamma * ( sqrt( abs( phib - type * ( V(b) - V(di) ))) - sqrt(phib));
	Fs_arg			=	( Vgsraw - ( Vt0x - Vdsi * delta * Fsat ) + aphit * 0.5 )/ ( 1.1 * nphit );
	if (Fs_arg <= `LARGE_VALUE) begin
		Fs		=	1.0 + exp( Fs_arg );
		FFx		=	Vgsraw - nphit * ln( Fs );
	end
	else begin
		Fs		=	0.0; 								//	Not used
		FFx		=	Vgsraw - nphit * Fs_arg;
	end
	Fd_arg			=	( Vgdraw - ( Vt0y - Vdsi * delta * Fsat ) + aphit * 0.5 )/ ( 1.1 * nphit );
	if (Fd_arg <= `LARGE_VALUE) begin
		Fd		=	1.0 + exp( Fd_arg );
		FFy		=	Vgdraw - nphit * ln( Fd );
	end
	else begin
		Fd		=	0.0; 								//	Not used
		FFy		=	Vgdraw - nphit * Fd_arg;
	end
	Qsif			=	( type * Cif + CC * Vgsraw ) * FFx;
	Qdif			=	( type * Cif + CC * Vgdraw ) * FFy;
		
	
	//Partitioned charge
	Qs			=	-W * ( Qinvs + Qsov + Qsif ); 					// 	s-terminal charge
	Qd			=	-W * ( Qinvd + Qdov + Qdif ); 					// 	d-terminal charge
	Qg			=	-( Qs + Qd + Qb );						// 	g-terminal charge


	//Sub-circuit initialization
    	I(di,si) 		<+ 	type * dir * Id;
    	I(d,di)			<+ 	( V(d) - V(di) )/ Rd;
    	I(si,s) 		<+ 	( V(si) - V(s) )/ Rs;

	I(si,b) 		<+ 	ddt( Qs );							// 	charge term: node si to node b
	I(di,b) 		<+ 	ddt( Qd );							// 	charge term: node di to node b
	I(g,b) 			<+ 	ddt( Qg );							// 	charge term: node g to node b


end
endmodule