1/* $Id: d_mos5.model,v 26.92 2008/08/23 05:40:00 al Exp $ -*- C++ -*-
2 * Copyright (C) 2001 Albert Davis
3 * Author: Albert Davis <aldavis@gnu.org>
4 *
5 * This file is part of "Gnucap", the Gnu Circuit Analysis Package
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 3, or (at your option)
10 * any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
20 * 02110-1301, USA.
21 *------------------------------------------------------------------
22 * Spice BSIM2 model
23 * derived from Spice3f4,Copyright 1990 Regents of the University of California
24 * 1988 Min-Chie Jeng, Hong J. Park, Thomas L. Quarles
25 * Recoded for Gnucap model compiler, Al Davis, 2000
26 */
27h_headers {
28#include "d_mos_base.h"
29}
30cc_headers {
31}
32/*--------------------------------------------------------------------------*/
33model BUILT_IN_MOS5 {
34  level 5;
35  public_keys {
36    nmos5 polarity=pN;
37    pmos5 polarity=pP;
38  }
39  dev_type BUILT_IN_MOS;
40  inherit BUILT_IN_MOS_BASE;
41  independent {
42    override {
43      double mjsw "" final_default=.33;
44      double pb "" final_default=0.1 quiet_min=0.1;
45      double pbsw "" final_default=pb quiet_min=0.1;
46      double cjo "" default=0.0;
47      int cmodel "CMODEL" print_test="cmodel!=1"
48	calculate="((!cmodel)?1:cmodel)";
49      int mos_level "back-annotate for diode" name=DIODElevel
50	print_test="mos_level != LEVEL" default=LEVEL;
51    }
52    raw_parameters {
53      double dl_u "Channel length reduction"
54	name=DL default=0.;
55      double dw_u "Channel width reduction"
56	name=DW default=0.;
57      double tox_u "Gate oxide thickness"
58	name=TOX default=0. quiet_min=1e-20;
59      double vdd "Max Vds"
60	name=VDD default=0.;
61      double vgg "Max Vgs"
62	name=VGG default=0.;
63      double vbb "Max Vbs"
64	name=VBB default=0.;
65      double wdf "Default width of source drain diffusion (ignored)"
66	name=WDF default=0.;
67      double dell "Length reduction of source drain diff (ignored)"
68	name=DELL default=0.;
69      double temp_c "temperature"
70	name=TEMP default=27.;
71      double xpart "Flag for channel charge partitioning"
72	name=XPART default=0.;
73    }
74    calculated_parameters {
75      double dl "" calculate="dl_u*MICRON2METER";
76      double dw "" calculate="dw_u*MICRON2METER";
77      double tox "" calculate="tox_u*MICRON2METER";
78      double cox "" calculate="3.453e-11 /*E_OX*/ / tox";
79      double vdd2 "" calculate="2 * vdd";
80      double vgg2 "" calculate="2 * vgg";
81      double vbb2 "" calculate="2 * vbb";
82      double Vtm "" calculate="8.625e-5 /*K/Q*/ * (temp_c + P_CELSIUS0 -.15)";
83    }
84  }
85  size_dependent {
86    raw_parameters {
87      double phi "Strong inversion surface potential"
88	name=PHI default=0.;
89      double vfb "flat band voltage at given L and W"
90	name=VFB default=0.;
91      double k1  "bulk effect coefficient 1"
92	name=K1 default=0.;
93      double k2  "bulk effect coefficient 2"
94	name=K2 default=0.;
95      double eta0 "drain induced barrier lowering"
96	name=ETA0 default=0.;
97      double etaB "Vbs dependence of Eta"
98	name=ETAB default=0.;
99
100      double mob0 "" name=MU0 default=0.;
101      double mob0B "" name=MU0B default=0.;
102      double mobs0 "" name=MUS0 default=0.;
103      double mobsB "" name=MUSB default=0.;
104      double mob20 "" name=MU20 default=0.;
105      double mob2B "" name=MU2B default=0.;
106      double mob2G "" name=MU2G default=0.;
107      double mob30 "" name=MU30 default=0.;
108      double mob3B "" name=MU3B default=0.;
109      double mob3G "" name=MU3G default=0.;
110      double mob40 "" name=MU40 default=0.;
111      double mob4B "" name=MU4B default=0.;
112      double mob4G "" name=MU4G default=0.;
113
114      double ua0 "Linear Vgs dependence of Mobility"
115	name=UA0 default=0.;
116      double uaB "Vbs dependence of Ua"
117	name=UAB default=0.;
118      double ub0 "Quadratic Vgs dependence of Mobility"
119	name=UB0 default=0.;
120      double ubB "Vbs dependence of Ub"
121	name=UBB default=0.;
122      double u10 "Drift Velocity Saturation due to Vds"
123	name=U10 default=0.;
124      double u1B "Vbs dependence of U1"
125	name=U1B default=0.;
126      double u1D "Vds dependence of U1"
127	name=U1D default=0.;
128      double n0  "Subthreshold slope at Vds=0, Vbs=0"
129	name=N0 default=0. positive;
130      double nB  "Vbs dependence of n"
131	name=NB default=0.;
132      double nD  "Vds dependence of n"
133	name=ND default=0.;
134      double vof0 "Vth offset at Vds=0, Vbs=0"
135	name=VOF0 default=0.;
136      double vofB "Vbs dependence of Vof"
137	name=VOFB default=0.;
138      double vofD "Vds dependence of Vof"
139	name=VOFD default=0.;
140      double ai0 "Pre-factor in hot-electron effects"
141	name=AI0 default=0.;
142      double aiB "Vbs dependence of Ai"
143	name=AIB default=0.;
144      double bi0 "Exp-factor in hot-electron effects"
145	name=BI0 default=0.;
146      double biB "Vbs dependence of Bi"
147	name=BIB default=0.;
148      double vghigh "Upper bound of cubic spline function"
149	name=VGHIGH default=0.;
150      double vglow  "Lower bound of cubic spline function"
151	name=VGLOW default=0.;
152    }
153    calculated_parameters {
154      double beta0  "Beta at Vds = 0 and Vgs = Vth"
155	calculate="mob0 * CoxWoverL";
156      double beta0B "Vbs dependence of Beta0"
157	calculate="mob0B * CoxWoverL";
158      double betas0  "Beta at Vds=Vdd and Vgs=Vth"
159	calculate="mobs0 * CoxWoverL" quiet_min="1.01*beta0";
160      double betasB "Vbs dependence of Betas"
161	calculate="mobsB * CoxWoverL";
162      double beta20 "Vds dependence of Beta in tanh term"
163	calculate="mob20";
164      double beta2B "Vbs dependence of Beta2"
165	calculate="mob2B";
166      double beta2G "Vgs dependence of Beta2"
167	calculate="mob2G";
168      double beta30 "Vds dependence of Beta in linear term"
169	calculate="mob30 * CoxWoverL";
170      double beta3B "Vbs dependence of Beta3"
171	calculate="mob3B * CoxWoverL";
172      double beta3G "Vgs dependence of Beta3"
173	calculate="mob3G * CoxWoverL";
174      double beta40 "Vds dependence of Beta in quadra term"
175	calculate="mob40 * CoxWoverL";
176      double beta4B "Vbs dependence of Beta4"
177	calculate="mob4B * CoxWoverL";
178      double beta4G "Vgs dependence of Beta4"
179	calculate="mob4G * CoxWoverL";
180      double Phis3 "" calculate="sqrt(phi) * phi";
181      double One_Third_CoxWL "" calculate="cgate / 3.0";
182      double Two_Third_CoxWL "" calculate="2.0 * One_Third_CoxWL";
183      double Arg;
184    }
185    code_pre {
186      l_eff -= m->dl;
187      w_eff -= m->dw;
188      cgate = m->cox * w_eff * l_eff;
189      double L = l_eff/MICRON2METER;
190      double W = w_eff/MICRON2METER;
191      double CoxWoverL = 1e-4 * m->cox * w_eff / l_eff;
192    }
193    override {
194      double cgate "" calculate="m->cox * w_eff * l_eff";
195    }
196    code_post {
197      double tmp = betas0 - beta0 - beta0B * m->vbb;
198      if ((-betasB * m->vbb) > tmp) {
199	untested();
200	betasB = -tmp / m->vbb;
201      }
202      Arg = betasB - beta0B - m->vdd * (beta3B - m->vdd * beta4B);
203    }
204  }
205  /*-----------------------------------------------------------------------*/
206  tr_eval {
207    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
208    trace3("", d->vds, d->vgs, d->vbs);
209    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
210    d->reverse_if_needed();
211    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
212    trace4("", c->lo, m->dl, c->wo, m->dw);
213
214    double Vbs = std::max(m->vbb2, d->vbs);
215    double Vgs = std::min(m->vgg2, d->vgs);
216    double Vds = std::min(m->vdd2, d->vds);
217    trace3("", Vbs, Vgs, Vds);
218
219    /* Threshold Voltage. */
220    double Phisb, dPhisb_dVb, T1s, dT1s_dVb;
221    if (Vbs <= 0.0) {
222      d->sbfwd = false;
223      Phisb = s->phi - Vbs;
224      dPhisb_dVb = -1.0;
225      T1s = sqrt(Phisb);
226      dT1s_dVb = -0.5 / T1s;
227      trace4("-", Phisb, dPhisb_dVb, T1s, dT1s_dVb);
228    }else{
229      d->sbfwd = true;
230      double tmp = s->phi / (s->phi + Vbs);
231      Phisb = s->phi * tmp;
232      dPhisb_dVb = -tmp * tmp;
233      T1s = s->Phis3 / (s->phi + 0.5 * Vbs);
234      dT1s_dVb = -0.5 * T1s * T1s / s->Phis3;
235      trace4("+", Phisb, dPhisb_dVb, T1s, dT1s_dVb);
236    }
237
238    double Eta = s->eta0 + s->etaB * Vbs;
239    double Ua = s->ua0 + s->uaB * Vbs;
240    double Ub = s->ub0 + s->ubB * Vbs;
241    double U1s = s->u10 + s->u1B * Vbs;
242    trace4("", Eta, Ua, Ub, U1s);
243
244    d->von = s->vfb + s->phi + s->k1*T1s - s->k2*Phisb - Eta*Vds;
245    double dVth_dVd = -Eta;
246    double dVth_dVb = s->k1 * dT1s_dVb + s->k2 - s->etaB * Vds;
247    d->vgst = Vgs - d->von;
248    trace4("", d->von, dVth_dVd, dVth_dVb, d->vgst);
249
250    double Aa, dAa_dVb, Inv_Aa;
251    {
252      double tmp = 1.0 / (1.744 + 0.8364 * Phisb);
253      double Gg = 1.0 - tmp;
254      double dGg_dVb = 0.8364 * tmp * tmp * dPhisb_dVb;
255      double T0 = Gg / T1s;
256      trace4("", tmp, Gg, dGg_dVb, T0);
257      double tmp1 = 0.5 * T0 * s->k1;
258      Aa = 1.0 + tmp1;
259      dAa_dVb = (Aa - 1.0) * (dGg_dVb / Gg - dT1s_dVb / T1s);
260      Inv_Aa = 1.0 / Aa;
261      trace4("", tmp1, Aa, dAa_dVb, Inv_Aa);
262    }
263
264    trace3("", s->vghigh, s->vglow, d->vgst);
265
266    double Exp0, Exp1, n, Vgeff, dVgeff_dVg, dVgeff_dVd, dVgeff_dVb;
267    if ((d->vgst >= s->vghigh) || (s->n0 == 0.0)) {
268      Exp0 = NOT_VALID;
269      Exp1 = NOT_VALID;
270      n = NOT_VALID;
271      Vgeff = d->vgst;
272      dVgeff_dVg = 1.0;
273      dVgeff_dVd = -dVth_dVd;
274      dVgeff_dVb = -dVth_dVb;
275      trace0("vgst>vghigh");
276    }else{
277      double Vof = s->vof0 + s->vofB * Vbs + s->vofD * Vds;
278      n = s->n0 + s->nB / T1s + s->nD * Vds;
279      double tmp = 0.5 / (n * m->Vtm);
280      trace3("", Vof, n, tmp);
281
282      double ExpArg1 = -Vds / m->Vtm;
283      ExpArg1 = std::max(ExpArg1, -30.0);
284      Exp1 = exp(ExpArg1);
285      double tmp1 = 1.0 - Exp1;
286      tmp1 = std::max(tmp1, 1.0e-18);
287      double tmp2 = 2.0 * Aa * tmp1;
288      trace4("", ExpArg1, Exp1, tmp1, tmp2);
289
290      // exports Exp0, Vgeff, dVgeff_dVg, dVgeff_dVd, dVgeff_dVb
291      if (d->vgst <= s->vglow) {
292	double ExpArg = d->vgst * tmp;
293	ExpArg = std::max(ExpArg, -30.0);
294	Exp0 = exp(0.5 * Vof + ExpArg);
295	Vgeff = sqrt(tmp2) * m->Vtm * Exp0;
296	double T0 = n * m->Vtm;
297	dVgeff_dVg = Vgeff * tmp;
298	dVgeff_dVd = dVgeff_dVg * (n / tmp1 * Exp1 - dVth_dVd
299				   - d->vgst * s->nD / n + T0 * s->vofD);
300	dVgeff_dVb = dVgeff_dVg * (s->vofB * T0 - dVth_dVb
301				 + s->nB * d->vgst / (n * T1s * T1s) * dT1s_dVb
302				 + T0 * Inv_Aa * dAa_dVb);
303	trace0("vgst<vglow");
304      }else{
305	double ExpArg = s->vglow * tmp;
306	ExpArg = std::max(ExpArg, -30.0);
307	Exp0 = exp(0.5 * Vof + ExpArg);
308	Vgeff = sqrt(2.0 * Aa * (1.0 - Exp1)) * m->Vtm * Exp0;
309	double Con1 = s->vghigh;
310	double Con3 = Vgeff;
311	double Con4 = Con3 * tmp;
312	double SqrVghigh = s->vghigh * s->vghigh;
313	double SqrVglow = s->vglow * s->vglow;
314	double CubVghigh = s->vghigh * SqrVghigh;
315	double CubVglow = s->vglow * SqrVglow;
316	double T0 = 2.0 * s->vghigh;
317	double T1 = 2.0 * s->vglow;
318	double T2 = 3.0 * SqrVghigh;
319	double T3 = 3.0 * SqrVglow;
320	double T4 = s->vghigh - s->vglow;
321	double T5 = SqrVghigh - SqrVglow;
322	double T6 = CubVghigh - CubVglow;
323	double T7 = Con1 - Con3;
324	double delta = (T1-T0) * T6 + (T2-T3) * T5 + (T0*T3 - T1*T2) * T4;
325	delta = 1.0 / delta;
326	double Coeffb = (T1 - Con4 * T0) * T6 + (Con4 * T2 - T3) * T5
327	  + (T0 * T3 - T1 * T2) * T7;
328	double Coeffc = (Con4-1.0) * T6 + (T2-T3) * T7 + (T3 - Con4*T2) * T4;
329	double Coeffd = (T1-T0) * T7 + (1.0-Con4) * T5 + (Con4*T0 - T1) * T4;
330	double Coeffa = SqrVghigh * (Coeffc + Coeffd * T0);
331	Vgeff = (Coeffa + d->vgst * (Coeffb + d->vgst*(Coeffc+d->vgst*Coeffd)))
332	  *delta;
333	dVgeff_dVg = (Coeffb + d->vgst*(2.0*Coeffc+3.0*d->vgst*Coeffd)) *delta;
334	T7 = Con3 * tmp;
335	double T8 = dT1s_dVb * s->nB / (T1s * T1s * n);
336	double T9 = n * m->Vtm;
337	double dCon3_dVd = T7*(n*Exp1/tmp1 -s->vglow*s->nD/n + T9*s->vofD);
338	double dCon3_dVb = T7*(T9*Inv_Aa*dAa_dVb + s->vglow*T8 + T9*s->vofB);
339	double dCon4_dVd = tmp * dCon3_dVd - T7 * s->nD / n;
340	double dCon4_dVb = tmp * dCon3_dVb + T7 * T8;
341
342	double dCoeffb_dVd = dCon4_dVd*(T2*T5-T0*T6) + dCon3_dVd*(T1*T2-T0*T3);
343	double dCoeffc_dVd = dCon4_dVd * (T6 - T2*T4) + dCon3_dVd * (T3 - T2);
344	double dCoeffd_dVd = dCon4_dVd * (T0*T4 - T5) + dCon3_dVd * (T0 - T1);
345	double dCoeffa_dVd = SqrVghigh * (dCoeffc_dVd + dCoeffd_dVd * T0);
346
347	dVgeff_dVd = -dVgeff_dVg * dVth_dVd + (dCoeffa_dVd
348		+ d->vgst * (dCoeffb_dVd + d->vgst
349			    * (dCoeffc_dVd + d->vgst * dCoeffd_dVd))) * delta;
350
351	double dCoeffb_dVb = dCon4_dVb*(T2*T5-T0*T6) + dCon3_dVb*(T1*T2-T0*T3);
352	double dCoeffc_dVb = dCon4_dVb * (T6 - T2*T4) + dCon3_dVb * (T3 - T2);
353	double dCoeffd_dVb = dCon4_dVb * (T0*T4 - T5) + dCon3_dVb * (T0 - T1);
354	double dCoeffa_dVb = SqrVghigh * (dCoeffc_dVb + dCoeffd_dVb * T0);
355
356	dVgeff_dVb = -dVgeff_dVg * dVth_dVb + (dCoeffa_dVb
357		+ d->vgst * (dCoeffb_dVb + d->vgst
358			     * (dCoeffc_dVb + d->vgst * dCoeffd_dVb))) * delta;
359	trace0("else");
360      }
361    }
362    trace3("", Exp0, Exp1, n);
363    trace4("", Vgeff, dVgeff_dVg, dVgeff_dVd, dVgeff_dVb);
364
365    double dVdsat_dVd, dVdsat_dVg, dVdsat_dVb;
366    if (Vgeff > 0.0) {	// normal operation
367      d->cutoff = false;
368      double Uvert = 1.0 + Vgeff * (Ua + Vgeff * Ub);
369      Uvert = std::max(Uvert, 0.2);
370      double Inv_Uvert = 1.0 / Uvert;
371      double dUvert_dVg, dUvert_dVd, dUvert_dVb;
372      {
373	double T8 = Ua + 2.0 * Ub * Vgeff;
374	dUvert_dVg = T8 * dVgeff_dVg;
375	dUvert_dVd = T8 * dVgeff_dVd;
376	dUvert_dVb = T8 * dVgeff_dVb + Vgeff * (s->uaB + Vgeff * s->ubB);
377	trace2("", T8, Uvert);
378	trace3("", dUvert_dVg, dUvert_dVd, dUvert_dVb);
379      }
380
381      double Vc, dVc_dVg, dVc_dVd, dVc_dVb;
382      {
383	double T8 = U1s * Inv_Aa * Inv_Uvert;
384	Vc = T8 * Vgeff;
385	double T9 = Vc * Inv_Uvert;
386	dVc_dVg = T8 * dVgeff_dVg - T9 * dUvert_dVg;
387	dVc_dVd = T8 * dVgeff_dVd - T9 * dUvert_dVd;
388	dVc_dVb = T8 * dVgeff_dVb
389	  + s->u1B * Vgeff * Inv_Aa * Inv_Uvert
390	  - Vc * Inv_Aa * dAa_dVb
391	  - T9 * dUvert_dVb;
392	trace3("", T8, T9, Vc);
393	trace3("", dVc_dVg, dVc_dVd, dVc_dVb);
394      }
395
396      double Kk, dKk_dVc;
397      {
398	double tmp2 = sqrt(1.0 + 2.0 * Vc);
399	Kk = 0.5 * (1.0 + Vc + tmp2);
400	dKk_dVc = 0.5  + 0.5 / tmp2;
401	trace3("", tmp2, Kk, dKk_dVc);
402      }
403
404      {
405	double T8 = Inv_Aa / sqrt(Kk);
406	d->vdsat = std::max(Vgeff * T8, 1.0e-18);
407	double T9 = 0.5 * d->vdsat * dKk_dVc / Kk;
408	dVdsat_dVd = T8 * dVgeff_dVd - T9 * dVc_dVd;
409	dVdsat_dVg = T8 * dVgeff_dVg - T9 * dVc_dVg;
410	dVdsat_dVb = T8 * dVgeff_dVb - T9 * dVc_dVb - d->vdsat*Inv_Aa*dAa_dVb;
411	trace3("", T8, T9, d->vdsat);
412	trace3("", dVdsat_dVd, dVdsat_dVg, dVdsat_dVb);
413      }
414
415      double Beta, dBeta_dVd, dBeta_dVg, dBeta_dVb;
416      {
417	double Beta0 = s->beta0  + s->beta0B * Vbs;
418	double Betas = s->betas0 + s->betasB * Vbs;
419	double Beta2 = s->beta20 + s->beta2B * Vbs + s->beta2G * Vgs;
420	double Beta3 = s->beta30 + s->beta3B * Vbs + s->beta3G * Vgs;
421	double Beta4 = s->beta40 + s->beta4B * Vbs + s->beta4G * Vgs;
422	double Beta1 = Betas - (Beta0 + m->vdd * (Beta3 - m->vdd * Beta4));
423	trace4("", Beta0, s->beta0, s->beta0B, Vbs);
424	trace4("", Betas, s->betas0, s->betasB, Vgs);
425	trace4("", Beta2, s->beta20, s->beta2B, s->beta2G);
426	trace4("", Beta3, s->beta30, s->beta3B, s->beta3G);
427	trace4("", Beta4, s->beta40, s->beta4B, s->beta4G);
428	trace2("", Beta1, m->vdd);
429
430	double T0 = Vds * Beta2 / d->vdsat;
431	T0 = std::min(T0, 30.0);
432	double T1 = exp(T0);
433	double T2 = T1 * T1;
434	double T3 = T2 + 1.0;
435	trace4("", T0, T1, T2, T3);
436	double tanh = (T2 - 1.0) / T3;
437	double Sqrsech = 4.0 * T2 / (T3 * T3);
438	trace2("", tanh, Sqrsech);
439
440	Beta = Beta0 + Beta1 * tanh + Vds * (Beta3 - Beta4 * Vds);
441	double T4 = Beta1 * Sqrsech / d->vdsat;
442	double T5 = m->vdd * tanh;
443	dBeta_dVd = Beta3 - 2.0*Beta4*Vds + T4*(Beta2-T0*dVdsat_dVd);
444	dBeta_dVg = T4 * (s->beta2G * Vds - T0 * dVdsat_dVg)
445	  + s->beta3G * (Vds - T5)
446	  - s->beta4G * (Vds * Vds - m->vdd * T5);
447	double dBeta1_dVb = s->Arg;
448	dBeta_dVb = s->beta0B
449	  + dBeta1_dVb * tanh
450	  + Vds * (s->beta3B - Vds * s->beta4B)
451	  + T4 * (s->beta2B * Vds - T0 * dVdsat_dVb);
452	trace3("", T4, T5, dBeta1_dVb);
453	trace4("", Beta, dBeta_dVd, dBeta_dVg, dBeta_dVb);
454      }
455
456      if (d->vgst > s->vglow) {	// not subthreshold
457	d->subthreshold = false;
458	if (Vds <= d->vdsat) {		// triode region
459	  d->saturated = false;
460	  double T3 = Vds / d->vdsat;
461	  double T4 = T3 - 1.0;
462	  double T2 =  1.0 - s->u1D * T4 * T4;
463	  double U1 =  U1s * T2;
464	  double Utot = Uvert + U1 * Vds;
465	  Utot = std::max(Utot, 0.5);
466	  double Inv_Utot = 1.0 / Utot;
467	  double T5 = 2.0 * U1s * s->u1D * T4 / d->vdsat;
468	  double dU1_dVd = T5 * (T3 * dVdsat_dVd - 1.0);
469	  double dU1_dVg = T5 * T3 * dVdsat_dVg;
470	  double dU1_dVb = T5 * T3 * dVdsat_dVb + s->u1B * T2;
471	  double dUtot_dVd = dUvert_dVd + U1 + Vds * dU1_dVd;
472	  double dUtot_dVg = dUvert_dVg + Vds * dU1_dVg;
473	  double dUtot_dVb = dUvert_dVb + Vds * dU1_dVb;
474
475	  double tmp1 = (Vgeff - 0.5 * Aa * Vds);
476	  double tmp3 = tmp1 * Vds;
477	  double Betaeff = Beta * Inv_Utot;
478	  d->ids = Betaeff * tmp3;
479	  double T6 = d->ids / Betaeff * Inv_Utot;
480	  d->gds = T6 * (dBeta_dVd - Betaeff * dUtot_dVd)
481	    + Betaeff * (tmp1 + (dVgeff_dVd - 0.5 * Aa) * Vds);
482	  d->gmf = T6 * (dBeta_dVg - Betaeff * dUtot_dVg)
483	    + Betaeff * Vds * dVgeff_dVg;
484	  d->gmbf = T6 * (dBeta_dVb - Betaeff * dUtot_dVb)
485	    + Betaeff * Vds * (dVgeff_dVb - 0.5 * Vds * dAa_dVb);
486	}else{	// Saturation
487	  d->saturated = true;
488	  double Inv_Kk = 1.0 / Kk;
489	  double tmp1 = Vgeff * Inv_Aa * Inv_Kk;
490	  double tmp3 = 0.5 * Vgeff * tmp1;
491	  double Betaeff = Beta * Inv_Uvert;
492	  d->ids = Betaeff * tmp3;
493	  double T0 = d->ids / Betaeff * Inv_Uvert;
494	  double T1 = Betaeff * Vgeff * Inv_Aa * Inv_Kk;
495	  double T2 = d->ids * Inv_Kk * dKk_dVc;
496
497	  if (s->ai0 != 0.0) {
498	    double Ai = s->ai0 + s->aiB * Vbs;
499	    double Bi = s->bi0 + s->biB * Vbs;
500	    double T5 = Bi / (Vds - d->vdsat);
501	    T5 = std::min(T5, 30.0);
502	    double T6 = exp(-T5);
503	    double FR = 1.0 + Ai * T6;
504	    double T7 = T5 / (Vds - d->vdsat);
505	    double T8 = (1.0 - FR) * T7;
506	    double dFR_dVd = T8 * (dVdsat_dVd - 1.0);
507	    double dFR_dVg = T8 * dVdsat_dVg;
508	    double dFR_dVb = T8 * dVdsat_dVb
509	      + T6 * (s->aiB - Ai * s->biB / (Vds - d->vdsat));
510
511	    d->gds = (T0 * (dBeta_dVd - Betaeff * dUvert_dVd)
512		     + T1 * dVgeff_dVd - T2 * dVc_dVd) * FR + d->ids * dFR_dVd;
513	    d->gmf = (T0 * (dBeta_dVg - Betaeff * dUvert_dVg)
514		     + T1 * dVgeff_dVg - T2 * dVc_dVg) * FR + d->ids * dFR_dVg;
515	    d->gmbf = (T0 * (dBeta_dVb - Betaeff * dUvert_dVb)
516		      + T1 * dVgeff_dVb - T2 * dVc_dVb
517		      - d->ids * Inv_Aa * dAa_dVb) * FR + d->ids * dFR_dVb;
518	    d->ids *= FR;
519	  }else{
520	    d->gds = T0 * (dBeta_dVd - Betaeff * dUvert_dVd)
521	      + T1 * dVgeff_dVd - T2 * dVc_dVd;
522	    d->gmf = T0 * (dBeta_dVg - Betaeff * dUvert_dVg)
523	      + T1 * dVgeff_dVg - T2 * dVc_dVg;
524	    d->gmbf = T0 * (dBeta_dVb - Betaeff * dUvert_dVb)
525	      + T1 * dVgeff_dVb - T2 * dVc_dVb - d->ids * Inv_Aa * dAa_dVb;
526	  }
527	}
528      }else{	// subthreshold
529	d->subthreshold = true;
530	assert(Exp0 != NOT_VALID);
531	double T0 = Exp0 * Exp0;
532	assert(Exp1 != NOT_VALID);
533	double T1 = Exp1;
534	trace4("sub", Exp0, Exp1, T0, T1);
535	trace2("", n, m->Vtm);
536	d->ids = Beta * m->Vtm * m->Vtm * T0 * (1.0 - T1);
537	double T2 = d->ids / Beta;
538	double T4 = n * m->Vtm;
539	double T3 = d->ids / T4;
540	trace4("", d->ids, T2, T4, T3);
541	double FR, dFR_dVd, dFR_dVg, dFR_dVb;
542	if ((Vds > d->vdsat) && s->ai0 != 0.0) {
543	  d->saturated = true;
544	  double Ai = s->ai0 + s->aiB * Vbs;
545	  double Bi = s->bi0 + s->biB * Vbs;
546	  double T5 = Bi / (Vds - d->vdsat);
547	  trace3("", Ai, Bi, T5);
548	  T5 = std::min(T5, 30.0);
549	  double T6 = exp(-T5);
550	  FR = 1.0 + Ai * T6;
551	  double T7 = T5 / (Vds - d->vdsat);
552	  double T8 = (1.0 - FR) * T7;
553	  trace4("", T5, T6, T7, T8);
554	  dFR_dVd = T8 * (dVdsat_dVd - 1.0);
555	  dFR_dVg = T8 * dVdsat_dVg;
556	  dFR_dVb = T8 * dVdsat_dVb + T6 * (s->aiB-Ai*s->biB/(Vds-d->vdsat));
557	  trace4("ai0!=0", FR, dFR_dVd, dFR_dVg, dFR_dVb);
558	}else{
559	  d->saturated = false;
560	  FR = 1.0;
561	  dFR_dVd = 0.0;
562	  dFR_dVg = 0.0;
563	  dFR_dVb = 0.0;
564	  trace4("ai0==0", FR, dFR_dVd, dFR_dVg, dFR_dVb);
565	}
566
567	d->gds = (T2 * dBeta_dVd
568		  + T3 * (s->vofD * T4 - dVth_dVd - s->nD * d->vgst / n)
569		  + Beta * m->Vtm * T0 * T1) * FR + d->ids * dFR_dVd;
570	d->gmf = (T2 * dBeta_dVg + T3) * FR + d->ids * dFR_dVg;
571	d->gmbf = (T2 * dBeta_dVb + T3 * (s->vofB * T4 - dVth_dVb + s->nB
572		* d->vgst / (n * T1s * T1s) * dT1s_dVb)) * FR + d->ids*dFR_dVb;
573	d->ids *= FR;
574      }
575    }else{	// cutoff???
576      d->cutoff = true;
577      // reachable only if vghigh and vgst both negative
578      d->vdsat = 0.0;
579      dVdsat_dVd = dVdsat_dVg = dVdsat_dVb = 0.0;
580      d->ids = 0.0;
581      d->gmf  = 0.0;
582      d->gds = 0.0;
583      d->gmbf = 0.0;
584    }
585    trace4("", d->vdsat, dVdsat_dVd, dVdsat_dVg, dVdsat_dVb);
586    trace4("", d->ids, d->gmf, d->gds, d->gmbf);
587
588    /* Some Limiting of DC Parameters */
589    d->gds = std::max(d->gds, 1.0e-20);
590    d->ids = std::max(d->ids, 1.0e-50);
591    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
592    double Vbseff, dVbseff_dVb;
593    if (Vbs < 0.0) {
594      Vbseff = Vbs;
595      dVbseff_dVb = 1.0;
596    }else{
597      Vbseff = s->phi - Phisb;
598      dVbseff_dVb = -dPhisb_dVb;
599    }
600    trace3("", Vbs, Vbseff, dVbseff_dVb);
601
602    double Arg1 = Vgs - Vbseff - s->vfb;
603    double Arg2 = Arg1 - d->vgst;
604    trace2("", Arg1, Arg2);
605    double Qbulk = s->One_Third_CoxWL * Arg2;
606    double dQbulk_dVb = s->One_Third_CoxWL * (dVth_dVb - dVbseff_dVb);
607    double dQbulk_dVd = s->One_Third_CoxWL * dVth_dVd;
608    trace3("", Qbulk, dQbulk_dVb, dQbulk_dVd);
609    if (Arg1 <= 0.0) { // accumulation region
610      d->qgate = s->cgate * Arg1;
611      d->qbulk = -(d->qgate);
612      d->qdrn = 0.0;
613
614      d->cggb = s->cgate;
615      d->cgdb = 0.0;
616      d->cgsb = -d->cggb * (1.0 - dVbseff_dVb);
617
618      d->cdgb = 0.0;
619      d->cddb = 0.0;
620      d->cdsb = 0.0;
621
622      d->cbgb = -s->cgate;
623      d->cbdb = 0.0;
624      d->cbsb = -d->cgsb;
625      trace4("", d->qgate, d->cggb, d->cgdb, d->cgsb);
626      trace4("", d->qbulk, d->cbgb, d->cbdb, d->cbsb);
627      trace4("acc", d->qdrn, d->cdgb, d->cddb, d->cdsb);
628    }else if (d->vgst <= 0.0) { // subthreshold
629      double T2 = Arg1 / Arg2;
630      double T3 = T2 * T2 * (s->cgate - s->Two_Third_CoxWL * T2);
631
632      d->qgate = s->cgate * Arg1 * (1.0 - T2 * (1.0 - T2 / 3.0));
633      d->qbulk = -(d->qgate);
634      d->qdrn = 0.0;
635
636      d->cggb = s->cgate * (1.0 - T2 * (2.0 - T2));
637      d->cgdb = T3 * dVth_dVd;
638      double tmp = T3 * dVth_dVb - (d->cggb + T3) * dVbseff_dVb;
639      d->cgsb = -(d->cggb + d->cgdb + tmp);
640
641      d->cbgb = -d->cggb;
642      d->cbdb = -d->cgdb;
643      d->cbsb = -d->cgsb;
644
645      d->cdgb = 0.0;
646      d->cddb = 0.0;
647      d->cdsb = 0.0;
648      trace3("", T2, T3, tmp);
649      trace4("", d->qgate, d->cggb, d->cgdb, d->cgsb);
650      trace4("", d->qbulk, d->cbgb, d->cbdb, d->cbsb);
651      trace4("sub", d->qdrn, d->cdgb, d->cddb, d->cdsb);
652    }else{
653      double Vdsat; // changes dVdsat_dVd, dVdsat_dVg, dVdsat_dVb;
654      if (d->vgst < s->vghigh) {
655	double Uvert = 1.0 + d->vgst * (Ua + d->vgst * Ub);
656	Uvert = std::max(Uvert, 0.2);
657	double Inv_Uvert = 1.0 / Uvert;
658	double dUvert_dVg = Ua + 2.0 * Ub * d->vgst;
659	double dUvert_dVd = -dUvert_dVg * dVth_dVd;
660	double dUvert_dVb = -dUvert_dVg * dVth_dVb
661	  + d->vgst * (s->uaB + d->vgst * s->ubB);
662	trace2("", Uvert, Inv_Uvert);
663	trace3("", dUvert_dVg, dUvert_dVd, dUvert_dVb);
664
665	double T8 = U1s * Inv_Aa * Inv_Uvert;
666	double Vc = T8 * d->vgst;
667	double T9 = Vc * Inv_Uvert;
668	double dVc_dVg = T8 - T9 * dUvert_dVg;
669	double dVc_dVd = -T8 * dVth_dVd - T9 * dUvert_dVd;
670	double dVc_dVb = -T8 * dVth_dVb
671	  + s->u1B * d->vgst * Inv_Aa * Inv_Uvert
672	  - Vc * Inv_Aa * dAa_dVb - T9 * dUvert_dVb;
673	trace3("", T8, T9, Vc);
674	trace3("", dVc_dVg, dVc_dVd, dVc_dVb);
675
676	double tmp2 = sqrt(1.0 + 2.0 * Vc);
677	double Kk = 0.5 * (1.0 + Vc + tmp2);
678	double dKk_dVc = 0.5  + 0.5 / tmp2;
679	trace3("", tmp2, Kk, dKk_dVc);
680
681	T8 = Inv_Aa / sqrt(Kk);
682	Vdsat = d->vgst * T8;
683	T9 = 0.5 * Vdsat * dKk_dVc / Kk;
684	trace2("", T8, T9);
685	dVdsat_dVd = -T8 * dVth_dVd - T9 * dVc_dVd;
686	dVdsat_dVg = T8 - T9 * dVc_dVg;
687	dVdsat_dVb = -T8*dVth_dVb - T9*dVc_dVb - Vdsat*Inv_Aa*dAa_dVb;
688	trace2("new", d->vdsat, Vdsat);
689	trace3("", dVdsat_dVd, dVdsat_dVg, dVdsat_dVb);
690      }else{
691	Vdsat = d->vdsat;
692	trace2("keep", d->vdsat, Vdsat);
693	trace3("", dVdsat_dVd, dVdsat_dVg, dVdsat_dVb);
694	// keep dVdsat_dVd, dVdsat_dVg, dVdsat_dVb;
695      }
696      if (Vds >= Vdsat) {       /* saturation region */
697	d->cggb = s->Two_Third_CoxWL;
698	d->cgdb = -d->cggb * dVth_dVd + dQbulk_dVd;
699	double tmp = -d->cggb * dVth_dVb + dQbulk_dVb;
700	d->cgsb = -(d->cggb + d->cgdb + tmp);
701	trace1("", tmp);
702
703	d->cbgb = 0.0;
704	d->cbdb = -dQbulk_dVd;
705	d->cbsb = dQbulk_dVd + dQbulk_dVb;
706
707	d->cdgb = -0.4 * d->cggb;
708	d->cddb = -d->cdgb * dVth_dVd;
709	tmp = -d->cdgb * dVth_dVb;
710	d->cdsb = -(d->cdgb + d->cddb + tmp);
711	trace1("", tmp);
712
713	d->qbulk = -Qbulk;
714	d->qgate = s->Two_Third_CoxWL * d->vgst + Qbulk;
715	d->qdrn = d->cdgb * d->vgst;
716	trace4("", d->qgate, d->cggb, d->cgdb, d->cgsb);
717	trace4("", d->qbulk, d->cbgb, d->cbdb, d->cbsb);
718	trace4("sat", d->qdrn, d->cdgb, d->cddb, d->cdsb);
719      }else{       /* linear region  */
720	double T7 = Vds / Vdsat;
721	double T8 = d->vgst / Vdsat;
722	double T6 = T7 * T8;
723	double T9 = 1.0 - T7;
724	double Vgdt = d->vgst * T9;
725	double T0 = d->vgst / (d->vgst + Vgdt);
726	double T1 = Vgdt / (d->vgst + Vgdt);
727	double T5 = T0 * T1;
728	double T2 = 1.0 - T1 + T5;
729	double T3 = 1.0 - T0 + T5;
730	trace4("", T7, T8, T6, T9);
731	trace2("", Vgdt, T0);
732	trace4("", T1, T5, T2, T3);
733
734	double dVgdt_dVg = T9 + T6 * dVdsat_dVg;
735	double dVgdt_dVd = T6 * dVdsat_dVd - T8 -T9 * dVth_dVd;
736	double dVgdt_dVb = T6 * dVdsat_dVb -T9 * dVth_dVb;
737	trace3("", dVgdt_dVg, dVgdt_dVd, dVgdt_dVb);
738
739	d->qgate = s->Two_Third_CoxWL * (d->vgst + Vgdt
740				   - Vgdt * T0) + Qbulk;
741	d->qbulk = -Qbulk;
742	d->qdrn = -s->One_Third_CoxWL * (0.2 * Vgdt
743				    + 0.8 * d->vgst + Vgdt * T1
744				    + 0.2 * T5 * (Vgdt - d->vgst));
745
746	d->cggb = s->Two_Third_CoxWL * (T2 + T3 * dVgdt_dVg);
747	d->cgdb = s->Two_Third_CoxWL * (T3*dVgdt_dVd-T2*dVth_dVd) + dQbulk_dVd;
748	double tmp = dQbulk_dVb +s->Two_Third_CoxWL*(T3*dVgdt_dVb-T2*dVth_dVb);
749	d->cgsb = -(d->cggb + d->cgdb + tmp);
750	trace1("", tmp);
751
752	T2 = 0.8 - 0.4 * T1 * (2.0 * T1 + T0 + T0 * (T1 - T0));
753	T3 = 0.2 + T1 + T0 * (1.0 - 0.4 * T0 * (T1 + 3.0 * T0));
754	d->cdgb = -s->One_Third_CoxWL * (T2 + T3 * dVgdt_dVg);
755	d->cddb = s->One_Third_CoxWL * (T2 * dVth_dVd - T3 * dVgdt_dVd);
756	tmp = s->One_Third_CoxWL * (T2 * dVth_dVb - T3 * dVgdt_dVb);
757	d->cdsb = -(d->cdgb + tmp + d->cddb);
758	trace3("", T2, T3, tmp);
759
760	d->cbgb = 0.0;
761	d->cbdb = -dQbulk_dVd;
762	d->cbsb = dQbulk_dVd + dQbulk_dVb;
763	trace4("", d->qgate, d->cggb, d->cgdb, d->cgsb);
764	trace4("", d->qbulk, d->cbgb, d->cbdb, d->cbsb);
765	trace4("lin", d->qdrn, d->cdgb, d->cddb, d->cdsb);
766      }
767    }
768    if (d->reversed) {
769      d->ids *= -1;
770      d->gmr = d->gmf;
771      d->gmbr = d->gmbf;
772      d->gmf = d->gmbf = 0;
773    }else{
774      d->gmr = d->gmbr = 0.;
775    }
776  }
777}
778/*--------------------------------------------------------------------------*/
779/*--------------------------------------------------------------------------*/
780