hsm1eval1_0.c

ngspice又一个电子CAD仿真软件代码.功能更全

double  E0 ;
double  E1 , E1_dVbs , E1_dVds , E1_dVgs ;
double  E2 , E2_dVbs , E2_dVds , E2_dVgs ;
double  Etun , Etun_dVbs , Etun_dVds , Etun_dVgs ;
double  Egidl , Egidl_dVbs , Egidl_dVds , Egidl_dVgs ;
double  Igs , Igs_dVbs , Igs_dVds , Igs_dVgs ;
double  Igs_dVbse , Igs_dVdse , Igs_dVgse ;
double  Ilg , Ilg_dVbs , Ilg_dVds , Ilg_dVgs ;
double  Ilg_dVbse , Ilg_dVdse , Ilg_dVgse ;
double  Cox0 ;
double  Lgate ;
double  rp1 , rp1_dVds ;
/* connecting function */
double  FD2 , FD2_dVbs , FD2_dVds , FD2_dVgs ;
double  FMD , FMD_dVds ;
/* Phonon scattering */
double	Wgate ;
double	mueph ;
/* temporary vars.-- */
double  T0 , T1 , T2 , T3 , T4 , T5 , T6 , T7 , T8 , T9 ;
double  TX , TX_dVbs , TX_dVds , TX_dVgs ;
double  TY , TY_dVbs , TY_dVds , TY_dVgs ;
double  T1_dVbs , T1_dVds , T1_dVgs ;
double  T2_dVbs , T2_dVds , T2_dVgs ;
double  T3_dVbs , T3_dVds , T3_dVgs ;
double  T4_dVbs , T4_dVds , T4_dVgs ;
double  T5_dVbs , T5_dVds , T5_dVgs ;
double  T6_dVbs , T6_dVds , T6_dVgs ;
double  T7_dVbs , T7_dVds , T7_dVgs ;
double  T8_dVbs , T8_dVds , T8_dVgs ;
double  T9_dVbs , T9_dVds , T9_dVgs ;
double  T10 , T20 , T21 , T30 , T31 ;

/* intrinsic finging capacitance */
double	Cgdf, Cgdfring;


/*================ Start of executable code.=================*/

    flg_info = sIN.info ;

/*-----------------------------------------------------------*
* Change units into CGS.
* - This section may be moved to an interface routine.
*-----------------*/

/* device instances */
    sIN.xl     *= C_m2cm ;
    sIN.xw     *= C_m2cm ;
    sIN.as     *= C_m2cm_p2 ;
    sIN.ad     *= C_m2cm_p2 ;
    sIN.ps     *= C_m2cm ;
    sIN.pd     *= C_m2cm ;

/* model parameters */
    sIN.tox    *= C_m2cm ;
    sIN.xld    *= C_m2cm ;
    sIN.xwd    *= C_m2cm ;

    sIN.xj     *= C_m2cm ;

    sIN.lp     *= C_m2cm ;
    sIN.xpolyd *= C_m2cm ;
    sIN.tpoly  *= C_m2cm ;

    sIN.rs     *= C_m2cm ;
    sIN.rd     *= C_m2cm ;

    sIN.sc3    *= C_m2cm ;
    sIN.scp3   *= C_m2cm ;
    sIN.parl2  *= C_m2cm ;

    sIN.wfc    *= C_m2cm ;
    sIN.clm2   /= C_m2cm ;

    sIN.rpock1 *= C_m2cm_p1o2 ;

    sIN.qme1   *= C_m2cm ;
    sIN.qme3   *= C_m2cm ;

    sIN.gidl1  /= C_m2cm ;

    sIN.gleak1 /= C_m2cm ;

    sIN.cgso   /= C_m2cm ;
    sIN.cgdo   /= C_m2cm ;
    sIN.cgbo   /= C_m2cm ;

    sIN.js0    /= C_m2cm_p2 ;
    sIN.js0sw  /= C_m2cm ;

    sIN.cj     /= C_m2cm_p2 ;
    sIN.cjsw   /= C_m2cm ;
    sIN.cjswg  /= C_m2cm ;

/*-----------------------------------------------------------*
* Start of the routine. (label)
*-----------------*/
start_of_routine:


 
/*-----------------------------------------------------------*
* Temperature dependent constants. 
*-----------------*/
 
/* Inverse of the thermal voltage */
    beta    = C_QE / ( C_KB * sIN.temp ) ;
    beta2   = beta * beta ;
 
/* Band gap */
    Eg  = C_Eg0 - sIN.temp 
            * ( sIN.bgtmp1 + sIN.temp * sIN.bgtmp2 ) ;

    Eg300  = C_Eg0 - C_T300  
            * ( sIN.bgtmp1 + C_T300 * sIN.bgtmp2 ) ;

/* Intrinsic carrier concentration */
    Nin = C_Nin0 * pow( sIN.temp / C_T300  , 1.5e0 ) 
            * exp( - Eg / 2.0e0 * beta + Eg300 / 2.0e0 * C_b300 ) ;
 
/*-----------------------------------------------------------*
* Fixed part. 
*-----------------*/
 

/* Lgate in [cm] / [m] */
    Lgate = sIN.xl ;
    Wgate = sIN.xw ;
    
/* Phonon Scattering */
    T1 = log( Wgate ) ;
    T2 = sIN.w0 - T1 - sti_dlt ;
    T3 = T2 * T2 ;
    T4 = sqrt( T3 + 4.0 * sti_dlt * sIN.w0 ) ;
    T5 = sIN.w0 - ( T2 - T4 ) / 2 ;
    mueph =  sIN.mueph1 + sIN.mueph2 * T5 ;
 
/* Metallurgical channel geometry */
    Weff    = sIN.xw - 2.0e0 * sIN.xwd ;
    Leff    = sIN.xl - 2.0e0 * sIN.xld ;
    Leff_inv    = 1.0e0 / Leff ;
 
/* Flat band voltage */
    Vfb = sIN.vfbc ;
 
/* Surface impurity profile */

    if( Lgate > sIN.lp ){
      Nsub = ( sIN.nsubc * ( Lgate - sIN.lp ) 
             + sIN.nsubp * sIN.lp ) / Lgate ;
    } else {
      Nsub = sIN.nsubp
           + ( sIN.nsubp - sIN.nsubc ) * ( sIN.lp - Lgate ) / sIN.lp ;
    }

    q_Nsub  = C_QE * Nsub ;
 
/* 2 phi_B */
        /* @temp, with pocket */
    Pb2 =  2.0e0 / beta * log( Nsub / Nin );
        /* @300K, with pocket */
    Pb20 = 2.0e0 / C_b300 * log( Nsub / C_Nin0 ) ;
        /* @300K, w/o pocket */
    Pb2c = 2.0e0 / C_b300 * log( sIN.nsubc / C_Nin0 ) ;

/* Debye length */
    Ldby    = sqrt( C_ESI / beta / q_Nsub ) ;
 
/* Coefficient of the F function for bulk charge */
    cnst0   = q_Nsub * Ldby * C_SQRT_2 ;
 
/* cnst1: n_{p0} / p_{p0} */
    T1  = Nin / Nsub ;
    cnst1   = T1 * T1 ;
/* Cox (clasical) */
    Cox0 = C_EOX / sIN.tox ;



/*-----------------------------------------------------------*
* Exchange bias conditions according to MOS type.
* - Vxse are external biases for HiSIM. ( type=NMOS , Vds >= 0
*   are assumed.) 
*-----------------*/

    /*
    Vbse = sIN.type * sIN.vbs ;
    Vdse = sIN.type * sIN.vds ;
    Vgse = sIN.type * sIN.vgs ;
    Vbde = Vbse - Vdse ;
    */
    /* modified by K. M. for SPICE3f5 */
    Vbse = sIN.vbs;
    Vdse = sIN.vds;
    Vgse = sIN.vgs;
    Vbde = Vbse - Vdse;


/*---------------------------------------------------*
* Cramp too large biases. 
* -note: Quantities are extrapolated in PART-5.
*-----------------*/

    if ( Vbse < Vbs_min ) {
      flg_vbsc = -1 ;
      Vbsc = Vbs_min ;
    } else if ( Vbse > 0.0 ) {
      flg_vbsc = 1 ;
      T1 = Vbse / Vbs_max ;
      T2 = sqrt ( 1.0 + ( T1 * T1 ) ) ;
      Vbsc = Vbse / T2 ;
      Vbsc_dVbse = Vbs_max * Vbs_max 
                 / ( ( Vbs_max * Vbs_max + Vbse * Vbse ) * T2 ) ;
    } else {
      flg_vbsc =  0 ;
      Vbsc = Vbse ;
    }


    if ( Vdse > Vds_max ) {
      flg_vdsc = 1 ;
      Vdsc = Vds_max ;
    } else {
      flg_vdsc =  0 ;
      Vdsc = Vdse ;
    }

    if ( Vgse > Vgs_max ) {
      flg_vgsc = 1 ;
      Vgsc = Vgs_max ;
    } else {
      flg_vgsc =  0 ;
      Vgsc = Vgse ;
    }

    if ( Vbde < Vbd_min ) {
      flg_vbdc = -1 ;
      Vbdc = Vbd_min ;
    } else if ( Vbde > Vbd_max ) {
      Vbdc = Vbd_max ;
      flg_vbdc = 1 ;
    } else {
      Vbdc = Vbde ;
      flg_vbdc =  0 ;
    }

    if ( flg_vbsc == -1 || flg_vdsc != 0 || flg_vgsc != 0 ||
         flg_vbdc != 0 ) {
      flg_vxxc = 1 ;
    }


/*-------------------------------------------------------------------*
* Set flags. 
*-----------------*/

    flg_rsrd = 0 ;
    flg_iprv = 0 ;
    flg_pprv = 0 ;

    Rs  = sIN.rs / Weff ;
    Rd  = sIN.rd / Weff ;

    if ( Rs + Rd >= epsm10 && sIN.corsrd >= 1 ) {
        flg_rsrd  = 1 ;
    }

    if ( sIN.has_prv == 1 ) {

      Vbsc_dif = Vbsc - sIN.vbsc_prv ;
      Vdsc_dif = Vdsc - sIN.vdsc_prv ;
      Vgsc_dif = Vgsc - sIN.vgsc_prv ;

      sum_vdif  = fabs( Vbsc_dif ) + fabs( Vdsc_dif ) 
                + fabs( Vgsc_dif ) ;

      if ( sIN.coiprv >= 1 && sum_vdif <= vtol_iprv ) { flg_iprv = 1 ;}
      if ( sIN.copprv >= 1 && sum_vdif <= vtol_pprv ) { flg_pprv = 1 ;}
    }

    if ( flg_rsrd == 0  ) {
        lp_bs_max = 1 ;
        flg_iprv  = 0 ;
    }



/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*
* Bias loop: iteration to solve the system of equations of 
*            the small circuit taking account Rs and Rd.
* - Vxs are internal (or effective) biases.
* - Equations:
*     Vbs = Vbsc - Rs * Ids
*     Vds = Vdsc - ( Rs + Rd ) * Ids
*     Vgs = Vgsc - Rs * Ids
*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

/*-----------------------------------------------------------*
* Initial guesses for biases.
*-----------------*/

    if ( flg_iprv == 1 ) {

        sIN.ids_dvbs_prv = Fn_Max( 0.0 , sIN.ids_dvbs_prv ) ;
        sIN.ids_dvds_prv = Fn_Max( 0.0 , sIN.ids_dvds_prv ) ;
        sIN.ids_dvgs_prv = Fn_Max( 0.0 , sIN.ids_dvgs_prv ) ;

        dVbs = Vbsc_dif * ( 1.0 - 
               1.0 / ( 1.0 + Rs * sIN.ids_dvbs_prv ) ) ;
        dVds = Vdsc_dif * ( 1.0 - 
               1.0 / ( 1.0 + ( Rs + Rd ) * sIN.ids_dvds_prv ) ) ;
        dVgs = Vgsc_dif * ( 1.0 - 
               1.0 / ( 1.0 + Rs * sIN.ids_dvgs_prv ) ) ;

	/*
	Ids = sIN.type * sIN.ids_prv 
            + sIN.ids_dvbs_prv * dVbs 
            + sIN.ids_dvds_prv * dVds 
            + sIN.ids_dvgs_prv * dVgs ;
	*/
        Ids = sIN.ids_prv 
            + sIN.ids_dvbs_prv * dVbs 
            + sIN.ids_dvds_prv * dVds 
            + sIN.ids_dvgs_prv * dVgs ;

        T1  = ( Ids - sIN.ids_prv ) ;
        T2  = fabs( T1 ) ;
        if ( Ids_maxvar * sIN.ids_prv < T2 ) {
            Ids = sIN.ids_prv * ( 1.0 + Fn_Sgn( T1 ) * Ids_maxvar ) ;
        }

        if ( Ids < 0 ) {
          Ids = 0.0 ;
        }

    } else {
        Ids = 0.0 ;

        if ( flg_pprv == 1 ) {
            dVbs = Vbsc_dif ;
            dVds = Vdsc_dif ;
            dVgs = Vgsc_dif ;
        }
    }

    Vbs = Vbsc - Ids * Rs ;

    Vds = Vdsc - Ids * ( Rs + Rd ) ;
    if ( Vds * Vdsc <= 0.0 ) { Vds = 0.0 ; } 

    Vgs = Vgsc - Ids * Rs ;

    if ( flg_pprv == 1 ) {

        Ps0 = sIN.ps0_prv ;

        Ps0_dVbs = sIN.ps0_dvbs_prv ;
        Ps0_dVds = sIN.ps0_dvds_prv ;
        Ps0_dVgs = sIN.ps0_dvgs_prv ;

        Pds = sIN.pds_prv ;

        Pds_dVbs = sIN.pds_dvbs_prv ;
        Pds_dVds = sIN.pds_dvds_prv ;
        Pds_dVgs = sIN.pds_dvgs_prv ;
    }
 

/*-----------------------------------------------------------*
* start of the loop.
*-----------------*/

  for ( lp_bs = 1 ; lp_bs <= lp_bs_max ; lp_bs ++ ) {


    Ids_last = Ids ;


/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
* PART-1: Basic device characteristics. 
*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
 
/*-----------------------------------------------------------*
* Initialization. 
*-----------------*/
    /* Initialization of counters is needed for restart. */
    lp_s0   = 0 ;
    lp_sl   = 0 ;
 
/*-----------------------------------------------------------*
* Vxsz: Modified bias introduced to realize symmetry at Vds=0.
*-----------------*/

    T1 = exp( - Vbsc_dVbse * Vds / ( 2.0 * sIN.vzadd0 ) ) ;
    Vzadd = sIN.vzadd0 * T1 ;
    Vzadd_dVds = - 0.5 * Vbsc_dVbse * T1 ;


    if ( Vzadd < ps_conv ) {
        Vzadd = 0.0 ;
        Vzadd_dVds = 0.0 ;
    }

    Vbsz = Vbs + Vzadd ;
    Vbsz_dVbs = 1.0 ;
    Vbsz_dVds = Vzadd_dVds ;

    Vdsz = Vds + 2 * Vzadd ;
    Vdsz_dVds = 1.0 + 2 * Vzadd_dVds ;

    Vgsz = Vgs + Vzadd ;
    Vgsz_dVgs = 1.0 ;
    Vgsz_dVds = Vzadd_dVds ;

 
/*-----------------------------------------------------------*
* Quantum effect
*-----------------*/

    T1 = 2.0 * q_Nsub * C_ESI ;
    T2 = sqrt( T1 * ( Pb20 - Vbsz ) ) ;

    Vthq = Pb20 + Vfb + sIN.tox / C_EOX * T2 + sIN.qme2 ;

    T3 = - 0.5 * sIN.tox / C_EOX * T1 / T2 ;
    Vthq_dVbs = T3 * Vbsz_dVbs ;
    Vthq_dVds = T3 * Vbsz_dVds ;


    T1 = - Vfb - sIN.qme2 ;
    T2 = Vgsz - Vthq ;
    T3 = sIN.qme1 * T1 * T1 + sIN.qme3 ;
    T4 = sIN.qme1 * T2 * T2 + sIN.qme3 ;

    T5 = T4 - T3 - qme_dlt ;

    T6 = sqrt( T5 * T5 + 4.0 * qme_dlt * T4 ) ;

    dTox = T4 - 0.5 * ( T5 + T6 ) ;

    /* dTox_dT4 */
    T7 = 1.0 - 0.5 * ( 1.0 + ( T4 - T3 + qme_dlt ) / T6 ) ;

    T8 = 2.0 * sIN.qme1 * T2 * T7 ;

    dTox_dVbs = T8 * ( - Vthq_dVbs ) ;
    dTox_dVds = T8 * ( Vgsz_dVds - Vthq_dVds ) ;
    dTox_dVgs = T8 * ( Vgsz_dVgs ) ;

    Tox = sIN.tox + dTox ;
    Tox_dVbs = dTox_dVbs ;
    Tox_dVds = dTox_dVds ;
    Tox_dVgs = dTox_dVgs ;

    Cox = C_EOX / Tox ;
    T1  = - C_EOX / ( Tox * Tox ) ;
    Cox_dVbs = T1 * Tox_dVbs ;  
    Cox_dVds = T1 * Tox_dVds ;  
    Cox_dVgs = T1 * Tox_dVgs ;  

    Cox_inv  = Tox / C_EOX ; 
    T1  = 1.0 / C_EOX ;
    Cox_inv_dVbs = T1 * Tox_dVbs ;  
    Cox_inv_dVds = T1 * Tox_dVds ;  
    Cox_inv_dVgs = T1 * Tox_dVgs ;  


/*-----------------------------------------------------------*
* Threshold voltage. 
*-----------------*/

    Delta   = 0.1 ;

    Vbs1    = 2.0 - 0.25 * Vbsz ;
    Vbs2    = - Vbsz ;