1 /***********************************************************************
2 
3  HiSIM (Hiroshima University STARC IGFET Model)
4  Copyright (C) 2014 Hiroshima University & STARC
5 
6  MODEL NAME : HiSIM_HV
7  ( VERSION : 2  SUBVERSION : 2  REVISION : 0 )
8  Model Parameter 'VERSION' : 2.20
9  FILE : hsmhveval_qover.h
10 
11  DATE : 2014.6.11
12 
13  released by
14                 Hiroshima University &
15                 Semiconductor Technology Academic Research Center (STARC)
16 ***********************************************************************/
17 
18 /**********************************************************************
19 
20 The following source code, and all copyrights, trade secrets or other
21 intellectual property rights in and to the source code in its entirety,
22 is owned by the Hiroshima University and the STARC organization.
23 
24 All users need to follow the "HISIM_HV Distribution Statement and
25 Copyright Notice" attached to HiSIM_HV model.
26 
27 -----HISIM_HV Distribution Statement and Copyright Notice--------------
28 
29 Software is distributed as is, completely without warranty or service
30 support. Hiroshima University or STARC and its employees are not liable
31 for the condition or performance of the software.
32 
33 Hiroshima University and STARC own the copyright and grant users a perpetual,
34 irrevocable, worldwide, non-exclusive, royalty-free license with respect
35 to the software as set forth below.
36 
37 Hiroshima University and STARC hereby disclaims all implied warranties.
38 
39 Hiroshima University and STARC grant the users the right to modify, copy,
40 and redistribute the software and documentation, both within the user's
41 organization and externally, subject to the following restrictions
42 
43 1. The users agree not to charge for Hiroshima University and STARC code
44 itself but may charge for additions, extensions, or support.
45 
46 2. In any product based on the software, the users agree to acknowledge
47 Hiroshima University and STARC that developed the software. This
48 acknowledgment shall appear in the product documentation.
49 
50 3. The users agree to reproduce any copyright notice which appears on
51 the software on any copy or modification of such made available
52 to others."
53 
54 Toshimasa Asahara, President, Hiroshima University
55 Mitiko Miura-Mattausch, Professor, Hiroshima University
56 Katsuhiro Shimohigashi, President&CEO, STARC
57 June 2008 (revised October 2011)
58 *************************************************************************/
59 
60 /*  Begin HSMHV2evalQover */
61 
62     /*---------------------------------------------------*
63      * Clamp -Vxbgmt.
64      *-----------------*/
65     T0 = - Vxbgmt;
66     if ( T0 > Vbs_bnd ) {
67       T1 =    T0   - Vbs_bnd;
68       T1_dT =      - Vbs_bnd_dT;
69       T2 =    Vbs_max    - Vbs_bnd;
70       T2_dT = Vbs_max_dT - Vbs_bnd_dT;
71 
72       Fn_SUPoly4m( TY, T1, T2, T11, T0 );
73       TY_dT = T1_dT * T11 + T2_dT * T0;
74 
75       T10 = Vbs_bnd + TY ;
76       T10_dT = Vbs_bnd_dT + TY_dT ;
77     }  else {
78       T10 = T0 ;
79       T11 = 1.0 ;
80       T10_dT = 0.0;
81     }
82     Vxbgmtcl = - T10 - small2 ;
83     Vxbgmtcl_dVxbgmt = T11;
84     Vxbgmtcl_dT = - T10_dT;
85 
86     fac1 = cnst0over_func * Cox0_inv ;
87     fac1_dVbs = 0.0; fac1_dVds = 0.0; fac1_dVgs = 0.0;
88     fac1_dT = cnst0over_func_dT * Cox0_inv ;
89 
90     fac1p2 = fac1 * fac1 ;
91     fac1p2_dT = 2.0 * fac1 * fac1_dT ;
92 
93     VgpLD = - Vgbgmt + pParam->HSMHV2_vfbover;
94     VgpLD_dVgb = - 1.0e0 ;
95 
96     T0 = Nover_func / here->HSMHV2_nin ;
97     Pb2over = 2.0 / beta * log( T0 ) ;
98     T0_dT = - T0 / here->HSMHV2_nin * Nin_dT ;
99     Pb2over_dT = - Pb2over / beta * beta_dT + 2.0 / beta / T0 * T0_dT ;
100 
101     Vgb_fb_LD =  - Vxbgmtcl ;
102 
103     /*-----------------------------------*
104      * QsuLD: total charge = Accumulation | Depletion+inversion
105      *-----------------*/
106     if ( VgpLD  < Vgb_fb_LD ){
107       /*---------------------------*
108        * Accumulation
109        *-----------------*/
110       flg_ovzone = -1 ;
111       T1 = 1.0 / ( beta * cnst0over_func ) ;
112       T1_dT = - T1 * T1 * ( beta_dT * cnst0over_func + beta * cnst0over_func_dT ) ;
113       TY = T1 * Cox0 ;
114       Ac41 = 2.0 + 3.0 * C_SQRT_2 * TY ;
115       Ac4 = 8.0 * Ac41 * Ac41 * Ac41 ;
116       TY_dT = T1_dT  * Cox0 ;
117       Ac41_dT = 3.0 * C_SQRT_2 * TY_dT ;
118       Ac4_dT = 8.0 * 3.0 * Ac41 * Ac41 * Ac41_dT ;
119 
120       Ps0_min = here->HSMHV2_eg - Pb2over ;
121       Ps0_min_dT = Eg_dT - Pb2over_dT ;
122 
123       TX = beta * ( VgpLD + Vxbgmtcl ) ;
124       TX_dVxb = beta * Vxbgmtcl_dVxbgmt ;
125       TX_dVgb = beta * VgpLD_dVgb ;
126       TX_dT = beta_dT * ( VgpLD + Vxbgmtcl ) + beta * Vxbgmtcl_dT;
127 
128       Ac31 = 7.0 * C_SQRT_2 - 9.0 * TY * ( TX - 2.0 ) ;
129       Ac31_dVxb = - 9.0 * TY * TX_dVxb ;
130       Ac31_dVgb = - 9.0 * TY * TX_dVgb ;
131       Ac31_dT = - 9.0 * ( TY_dT * ( TX - 2.0 ) + TY * TX_dT );
132 
133       Ac3 = Ac31 * Ac31 ;
134       T1 = 2.0 * Ac31 ;
135       Ac3_dVxb = T1 * Ac31_dVxb ;
136       Ac3_dVgb = T1 * Ac31_dVgb ;
137       Ac3_dT = T1 * Ac31_dT ;
138 
139       if ( Ac4 < Ac3*1.0e-8 ) {
140         Ac1 = 0.5*Ac4/Ac31 ;
141         Ac1_dVxb =  - 0.5*Ac4/Ac3*Ac31_dVxb ;
142         Ac1_dVgb =  - 0.5*Ac4/Ac3*Ac31_dVxb ;
143         Ac1_dT = 0.5*Ac4_dT/Ac31 - 0.5*Ac4/Ac3*Ac31_dT ;
144       } else {
145         Ac2 = sqrt( Ac4 + Ac3 ) ;
146         T1 = 0.5 / Ac2 ;
147         Ac2_dVxb = T1 *  Ac3_dVxb ;
148         Ac2_dVgb = T1 *  Ac3_dVgb ;
149         Ac2_dT = T1 *  ( Ac4_dT + Ac3_dT );
150 
151         Ac1 = -Ac31 + Ac2 ;
152         Ac1_dVxb = Ac2_dVxb -Ac31_dVxb ;
153         Ac1_dVgb = Ac2_dVgb -Ac31_dVgb ;
154         Ac1_dT   = Ac2_dT -Ac31_dT ;
155       }
156 
157       Acd = pow( Ac1 , C_1o3 ) ;
158       T1 = C_1o3 / ( Acd * Acd ) ;
159       Acd_dVxb = Ac1_dVxb * T1 ;
160       Acd_dVgb = Ac1_dVgb * T1 ;
161       Acd_dT = Ac1_dT * T1 ;
162 
163       Acn = -4.0 * C_SQRT_2 - 12.0 * TY + 2.0 * Acd + C_SQRT_2 * Acd * Acd ;
164       T1 = 2.0 + 2.0 * C_SQRT_2 * Acd ;
165       Acn_dVxb = T1 * Acd_dVxb ;
166       Acn_dVgb = T1 * Acd_dVgb ;
167       Acn_dT = - 12.0 * TY_dT + T1 * Acd_dT ;
168 
169       Chi = Acn / Acd ;
170       T1 = 1.0 / ( Acd * Acd ) ;
171       Chi_dVxb = ( Acn_dVxb * Acd - Acn * Acd_dVxb ) * T1 ;
172       Chi_dVgb = ( Acn_dVgb * Acd - Acn * Acd_dVgb ) * T1 ;
173       Chi_dT = ( Acn_dT * Acd - Acn * Acd_dT ) * T1 ;
174 
175       Psa = Chi * beta_inv - Vxbgmtcl ;
176       Psa_dVxb = Chi_dVxb * beta_inv - Vxbgmtcl_dVxbgmt ;
177       Psa_dVgb = Chi_dVgb * beta_inv ;
178       Psa_dT = Chi_dT * beta_inv + Chi * beta_inv_dT - Vxbgmtcl_dT ;
179 
180       T1 = Psa + Vxbgmtcl ;
181       T1_dT = Psa_dT + Vxbgmtcl_dT ;
182       T2 = T1 / Ps0_min ;
183       T2_dT = ( T1_dT * Ps0_min - T1 * Ps0_min_dT ) / ( Ps0_min * Ps0_min ) ;
184       T3 = sqrt( 1.0 + ( T2 * T2 ) ) ;
185 
186       T9 = T2 / T3 / Ps0_min ;
187       T3_dVd = T9 * ( Psa_dVxb + Vxbgmtcl_dVxbgmt ) ;
188       T3_dVg = T9 * Psa_dVgb ;
189       T3_dT =  T2_dT * T2 / T3;
190 
191       Ps0LD = T1 / T3 - Vxbgmtcl ;
192       T9 = 1.0 / ( T3 * T3 ) ;
193       Ps0LD_dVxb = T9 * ( ( Psa_dVxb + Vxbgmtcl_dVxbgmt ) * T3 - T1 * T3_dVd ) - Vxbgmtcl_dVxbgmt ;
194       Ps0LD_dVgb = T9 * ( Psa_dVgb * T3 - T1 * T3_dVg );
195       Ps0LD_dT = T9 * ( T1_dT * T3 - T1 * T3_dT ) - Vxbgmtcl_dT;
196 
197       T2 = ( VgpLD - Ps0LD ) ;
198       QsuLD = Cox0 * T2 ;
199       QsuLD_dVxb = - Cox0 * Ps0LD_dVxb ;
200       QsuLD_dVgb = Cox0 * ( VgpLD_dVgb - Ps0LD_dVgb ) ;
201       QsuLD_dT = Cox0 * ( - Ps0LD_dT ) ;
202 
203       QbuLD = QsuLD ;
204       QbuLD_dVxb = QsuLD_dVxb ;
205       QbuLD_dVgb = QsuLD_dVgb ;
206       QbuLD_dT = QsuLD_dT ;
207 
208     } else {
209 
210       /*---------------------------*
211        * Depletion and inversion
212        *-----------------*/
213 
214       /* initial value for a few fixpoint iterations
215          to get Ps0_iniA from simplified Poisson equation: */
216        flg_ovzone = 2 ;
217        Chi = znbd3 ;
218        Chi_dVxb = 0.0 ; Chi_dVgb = 0.0 ; Chi_dT = 0.0 ;
219 
220        Ps0_iniA= Chi/beta - Vxbgmtcl ;
221        Ps0_iniA_dVxb = Chi_dVxb/beta - Vxbgmtcl_dVxbgmt ;
222        Ps0_iniA_dVgb = Chi_dVgb/beta ;
223        Ps0_iniA_dT   = Chi_dT/beta - Chi*beta_dT/(beta*beta) - Vxbgmtcl_dT;
224 
225       /* 1 .. 2 relaxation steps should be sufficient */
226       for ( lp_ld = 1; lp_ld <= 2; lp_ld ++ ) {
227         TY = exp(-Chi);
228         TY_dVxb = -Chi_dVxb * TY;
229         TY_dVgb = -Chi_dVgb * TY;
230         TY_dT   = - Chi_dT  * TY;
231         TX = 1.0e0 + 4.0e0
232            * ( beta * ( VgpLD + Vxbgmtcl ) - 1.0e0 + TY ) / ( fac1p2 * beta2 ) ;
233         TX_dVxb = 4.0e0 * ( beta * ( Vxbgmtcl_dVxbgmt ) + TY_dVxb ) / ( fac1p2 * beta2 );
234         TX_dVgb = 4.0e0 * ( beta * ( VgpLD_dVgb       ) + TY_dVgb ) / ( fac1p2 * beta2 );
235         T1 = ( beta * ( VgpLD + Vxbgmtcl ) - 1.0e0 + TY );
236         T1_dT = beta_dT * ( VgpLD + Vxbgmtcl ) + beta * Vxbgmtcl_dT + TY_dT;
237         T3 = fac1p2 * beta2 ;
238         T3_dT = fac1p2_dT * beta2 + fac1p2 * ( 2 * beta * beta_dT ) ;
239         TX_dT = 4 * ( T1_dT * T3 - T1 * T3_dT ) / ( T3 * T3 );
240         if ( TX < epsm10) {
241           TX = epsm10;
242           TX_dVxb = TX_dVgb = TX_dT = 0.0;
243         }
244 
245         Ps0_iniA = VgpLD + fac1p2 * beta / 2.0e0 * ( 1.0e0 - sqrt( TX ) ) ;
246         Ps0_iniA_dVxb =            - fac1p2 * beta / 2.0e0 * TX_dVxb * 0.5 / sqrt( TX );
247         Ps0_iniA_dVgb = VgpLD_dVgb - fac1p2 * beta / 2.0e0 * TX_dVgb * 0.5 / sqrt( TX );
248         T1 = fac1p2 * beta ;
249         T1_dT = fac1p2_dT * beta + fac1p2 * beta_dT ;
250         T2 = 1.0 - sqrt( TX );
251         T2_dT = - 1.0e0 / ( 2.0e0 * sqrt( TX ) ) * TX_dT ;
252         Ps0_iniA_dT = ( T1_dT * T2 + T1 * T2_dT ) / 2.0e0 ;
253 
254         Chi = beta * ( Ps0_iniA + Vxbgmtcl ) ;
255         Chi_dVxb = beta * ( Ps0_iniA_dVxb + Vxbgmtcl_dVxbgmt ) ;
256         Chi_dVgb = beta * ( Ps0_iniA_dVgb ) ;
257         Chi_dT = beta_dT * ( Ps0_iniA + Vxbgmtcl ) + beta * ( Ps0_iniA_dT + Vxbgmtcl_dT );
258       } /* End of iteration */
259 
260       if ( Chi < znbd3 ) {
261 
262         flg_ovzone = 1 ;
263 
264         /*-----------------------------------*
265          * zone-D1
266          * - Ps0_iniA is the analytical solution of QovLD=Qb0 with
267          *   Qb0 being approximated by 3-degree polynomial.
268          *
269          *   new: Inclusion of exp(-Chi) term at right border
270          *-----------------*/
271         Ta =  1.0/(9.0*sqrt(2.0)) - (5.0+7.0*exp(-3.0)) / (54.0*sqrt(2.0+exp(-3.0)));
272         Tb = (1.0+exp(-3.0)) / (2.0*sqrt(2.0+exp(-3.0))) - sqrt(2.0) / 3.0;
273         Tc =  1.0/sqrt(2.0) + 1.0/(beta*fac1);
274         Tc_dT = - (beta_dT*fac1 + beta*fac1_dT)/(beta2*fac1p2);
275         Td = - (VgpLD + Vxbgmtcl) / fac1;
276         Td_dVxb = - Vxbgmtcl_dVxbgmt / fac1;
277         Td_dVgb = - VgpLD_dVgb / fac1;
278         Td_dT   = - (Vxbgmtcl_dT*fac1 - (VgpLD+Vxbgmtcl)*fac1_dT)/fac1p2;
279         Tq = Tb*Tb*Tb / (27.0*Ta*Ta*Ta) - Tb*Tc/(6.0*Ta*Ta) + Td/(2.0*Ta);
280         Tq_dVxb = Td_dVxb/(2.0*Ta);
281         Tq_dVgb = Td_dVgb / (2.0*Ta);
282         Tq_dT   = - Tb/(6.0*Ta*Ta)*Tc_dT + Td_dT/(2.0*Ta);
283         Tp = (3.0*Ta*Tc-Tb*Tb)/(9.0*Ta*Ta);
284         Tp_dT = Tc_dT/(3.0*Ta);
285         T5      = sqrt(Tq*Tq + Tp*Tp*Tp);
286         T5_dVxb = 2.0*Tq*Tq_dVxb / (2.0*T5);
287         T5_dVgb = 2.0*Tq*Tq_dVgb / (2.0*T5);
288         T5_dT   = (2.0*Tq*Tq_dT + 3.0*Tp*Tp*Tp_dT) / (2.0*T5);
289         Tu = pow(-Tq + T5,C_1o3);
290         Tu_dVxb = Tu / (3.0 * (-Tq + T5)) * (-Tq_dVxb + T5_dVxb);
291         Tu_dVgb = Tu / (3.0 * (-Tq + T5)) * (-Tq_dVgb + T5_dVgb);
292         Tu_dT   = Tu / (3.0 * (-Tq + T5)) * (-Tq_dT   + T5_dT);
293         Tv = -pow(Tq + T5,C_1o3);
294         Tv_dVxb = Tv / (3.0 * (-Tq - T5)) * (-Tq_dVxb - T5_dVxb);
295         Tv_dVgb = Tv / (3.0 * (-Tq - T5)) * (-Tq_dVgb - T5_dVgb);
296         Tv_dT   = Tv / (3.0 * (-Tq - T5)) * (-Tq_dT   - T5_dT );
297         TX      = Tu + Tv - Tb/(3.0*Ta);
298         TX_dVxb = Tu_dVxb + Tv_dVxb;
299         TX_dVgb = Tu_dVgb + Tv_dVgb;
300         TX_dT   = Tu_dT   + Tv_dT  ;
301 
302         Ps0_iniA = TX * beta_inv - Vxbgmtcl ;
303         Ps0_iniA_dVxb = TX_dVxb * beta_inv - Vxbgmtcl_dVxbgmt;
304         Ps0_iniA_dVgb = TX_dVgb * beta_inv;
305 	Ps0_iniA_dT = TX_dT * beta_inv + TX * beta_inv_dT - Vxbgmtcl_dT;
306 
307         Chi = beta * ( Ps0_iniA + Vxbgmtcl ) ;
308         Chi_dVxb = beta * ( Ps0_iniA_dVxb + Vxbgmtcl_dVxbgmt ) ;
309         Chi_dVgb = beta * ( Ps0_iniA_dVgb ) ;
310         Chi_dT = beta_dT * ( Ps0_iniA + Vxbgmtcl ) + beta * ( Ps0_iniA_dT + Vxbgmtcl_dT );
311       }
312 
313       if ( model->HSMHV2_coqovsm > 0 ) {
314 	  /*-----------------------------------*
315 	   * - Ps0_iniB : upper bound.
316 	   *-----------------*/
317         flg_ovzone += 2;
318 
319         VgpLD_shift = VgpLD + Vxbgmtcl + 0.1;
320         VgpLD_shift_dVgb = VgpLD_dVgb;
321         VgpLD_shift_dVxb = Vxbgmtcl_dVxbgmt;
322         VgpLD_shift_dT   = Vxbgmtcl_dT;
323         exp_bVbs = exp( beta * - Vxbgmtcl ) + small;
324         exp_bVbs_dVxb = - exp_bVbs * beta * Vxbgmtcl_dVxbgmt;
325         exp_bVbs_dT   = - exp_bVbs * (beta_dT*Vxbgmtcl + beta*Vxbgmtcl_dT);
326         T0 = here->HSMHV2_nin / Nover_func;
327         T0_dT = Nin_dT / Nover_func;
328         cnst1over = T0 * T0;
329         cnst1over_dT = 2.0 * T0 * T0_dT;
330         gamma = cnst1over * exp_bVbs;
331         gamma_dVxb = cnst1over * exp_bVbs_dVxb;
332         gamma_dT   = cnst1over_dT * exp_bVbs + cnst1over * exp_bVbs_dT;
333 
334         T0    = beta2 * fac1p2;
335         T0_dT = 2.0 * beta * fac1 * (beta_dT*fac1+beta*fac1_dT);
336 
337         psi = beta*VgpLD_shift;
338         psi_dVgb = beta*VgpLD_shift_dVgb;
339         psi_dVxb = beta*VgpLD_shift_dVxb;
340         psi_dT   = beta_dT*VgpLD_shift + beta*VgpLD_shift_dT;
341         Chi_1      = log(gamma*T0 + psi*psi) - log(cnst1over*T0) + beta*Vxbgmtcl;
342         Chi_1_dVgb = 2.0*psi*psi_dVgb/ (gamma*T0 + psi*psi);
343         Chi_1_dVxb = (gamma_dVxb*T0+2.0*psi*psi_dVxb)/(gamma*T0+psi*psi)
344                             + beta*Vxbgmtcl_dVxbgmt;
345         Chi_1_dT   = (gamma_dT*T0+gamma*T0_dT+2.0*psi*psi_dT)/(gamma*T0+psi*psi)
346                             - (cnst1over_dT*T0 + cnst1over*T0_dT)/(cnst1over*T0)
347                             + beta_dT*Vxbgmtcl + beta*Vxbgmtcl_dT;
348 
349         Fn_SU2( Chi_1, Chi_1, psi, 1.0, T1, T2 );
350         Chi_1_dVgb = Chi_1_dVgb*T1 + psi_dVgb*T2;
351         Chi_1_dVxb = Chi_1_dVxb*T1 + psi_dVxb*T2;
352         Chi_1_dT   = Chi_1_dT  *T1 + psi_dT  *T2;
353 
354      /* 1 fixpoint step for getting more accurate Chi_B */
355         psi      -= Chi_1 ;
356         psi_dVgb -= Chi_1_dVgb ;
357         psi_dVxb -= Chi_1_dVxb ;
358         psi_dT   -= Chi_1_dT ;
359 
360         psi      += beta*0.1 ;
361         psi_dT   += beta_dT*0.1 ;
362 
363         psi_B = psi;
364         arg_B = psi*psi/(gamma*T0);
365         Chi_B = log(gamma*T0 + psi*psi) - log(cnst1over*T0) + beta*Vxbgmtcl;
366         Chi_B_dVgb = 2.0*psi*psi_dVgb/ (gamma*T0 + psi*psi);
367         Chi_B_dVxb = (gamma_dVxb*T0+2.0*psi*psi_dVxb)/(gamma*T0+psi*psi)
368                             + beta*Vxbgmtcl_dVxbgmt;
369         Chi_B_dT   = (gamma_dT*T0+gamma*T0_dT+2.0*psi*psi_dT)/(gamma*T0+psi*psi)
370                             - (cnst1over_dT*T0 + cnst1over*T0_dT)/(cnst1over*T0)
371                             + beta_dT*Vxbgmtcl + beta*Vxbgmtcl_dT;
372         Ps0_iniB      = Chi_B/beta - Vxbgmtcl ;
373         Ps0_iniB_dVgb = Chi_B_dVgb/beta;
374         Ps0_iniB_dVxb = Chi_B_dVxb/beta- Vxbgmtcl_dVxbgmt;
375         Ps0_iniB_dT   = Chi_B_dT/beta - Chi_B/(beta*beta)*beta_dT - Vxbgmtcl_dT;
376 
377 
378         /* construction of Ps0LD by taking Ps0_iniB as an upper limit of Ps0_iniA
379          *
380          * Limiting is done for Chi rather than for Ps0LD, to avoid shifting
381          * for Fn_SU2 */
382 
383         Chi_A = Chi;
384         Chi_A_dVxb = Chi_dVxb;
385         Chi_A_dVgb = Chi_dVgb;
386         Chi_A_dT   = Chi_dT;
387 
388         Fn_SU2( Chi, Chi_A, Chi_B, c_ps0ini_2*75.00, T1, T2 ); /* org: 50 */
389         Chi_dVgb = Chi_A_dVgb * T1 + Chi_B_dVgb * T2;
390         Chi_dVxb = Chi_A_dVxb * T1 + Chi_B_dVxb * T2;
391         Chi_dT   = Chi_A_dT   * T1 + Chi_B_dT   * T2;
392 
393       }
394 
395         /* updating Ps0LD */
396         Ps0LD= Chi/beta - Vxbgmtcl ;
397         Ps0LD_dVgb = Chi_dVgb/beta;
398         Ps0LD_dVxb = Chi_dVxb/beta- Vxbgmtcl_dVxbgmt;
399         Ps0LD_dT   = Chi_dT/beta - Chi/(beta*beta)*beta_dT - Vxbgmtcl_dT;
400 
401       T1      = Chi - 1.0 + exp(-Chi);
402       T1_dVxb = (1.0 - exp(-Chi)) * Chi_dVxb ;
403       T1_dVgb = (1.0 - exp(-Chi)) * Chi_dVgb ;
404       T1_dT   = (1.0 - exp(-Chi)) * Chi_dT   ;
405       if (T1 < epsm10) {
406          T1 = epsm10 ;
407          T1_dVxb = 0.0 ;
408          T1_dVgb = 0.0 ;
409          T1_dT   = 0.0 ;
410       }
411       T2 = sqrt(T1);
412       QbuLD = cnst0over_func * T2 ;
413       T3 = cnst0over_func * 0.5 / T2 ;
414       QbuLD_dVxb = T3 * T1_dVxb ;
415       QbuLD_dVgb = T3 * T1_dVgb ;
416       QbuLD_dT = T3 * T1_dT + cnst0over_func_dT * T2 ;
417 
418       /*-----------------------------------------------------------*
419        * QsuLD : Qovs or Qovd in unit area.
420        * note: QsuLD = Qdep+Qinv.
421        *-----------------*/
422       QsuLD = Cox0 * ( VgpLD - Ps0LD ) ;
423       QsuLD_dVxb = Cox0 * ( - Ps0LD_dVxb ) ;
424       QsuLD_dVgb = Cox0 * ( VgpLD_dVgb - Ps0LD_dVgb ) ;
425       QsuLD_dT = Cox0 * ( - Ps0LD_dT ) ;
426 
427       if ( model->HSMHV2_coqovsm == 1 ) { /* take initial values from analytical model */
428 
429 
430         /*---------------------------------------------------*
431          * Calculation of Ps0LD. (beginning of Newton loop)
432          * - Fs0 : Fs0 = 0 is the equation to be solved.
433          * - dPs0 : correction value.
434          *-----------------*/
435 
436         /* initial value too close to flat band should not be used */
437 //      if (Ps0LD+Vxbgmtcl < 1.0e-2) Ps0LD = 1.0e-2 - Vxbgmtcl;
438         exp_bVbs = exp( beta * - Vxbgmtcl ) ;
439         T0 = here->HSMHV2_nin / Nover_func;
440         cnst1over = T0 * T0;
441         cnst1over_dT = 2.0 * T0 * ( Nin_dT / Nover_func );
442         cfs1 = cnst1over * exp_bVbs ;
443 
444         flg_conv = 0 ;
445         for ( lp_s0 = 1 ; lp_s0 <= 2*lp_s0_max + 1 ; lp_s0 ++ ) {
446 
447             Chi = beta * ( Ps0LD + Vxbgmtcl ) ;
448 
449             if ( Chi < znbd5 ) {
450               /*-------------------------------------------*
451                * zone-D1/D2. (Ps0LD)
452                *-----------------*/
453               fi = Chi * Chi * Chi
454                 * ( cn_im53 + Chi * ( cn_im54 + Chi * cn_im55 ) ) ;
455               fi_dChi = Chi * Chi
456                 * ( 3 * cn_im53 + Chi * ( 4 * cn_im54 + Chi * 5 * cn_im55 ) ) ;
457 
458               fs01 = cfs1 * fi * fi ;
459               fs01_dPs0 = cfs1 * beta * 2 * fi * fi_dChi ;
460 
461               fb = Chi * ( cn_nc51
462                  + Chi * ( cn_nc52
463                  + Chi * ( cn_nc53
464                  + Chi * ( cn_nc54 + Chi * cn_nc55 ) ) ) ) ;
465               fb_dChi = cn_nc51
466                  + Chi * ( 2 * cn_nc52
467                  + Chi * ( 3 * cn_nc53
468                  + Chi * ( 4 * cn_nc54 + Chi * 5 * cn_nc55 ) ) ) ;
469 
470               fs02 = sqrt( fb * fb + fs01 + small ) ;
471               fs02_dPs0 = ( beta * fb_dChi * 2 * fb + fs01_dPs0 ) / ( fs02 + fs02 ) ;
472 
473             } else {
474              /*-------------------------------------------*
475               * zone-D3. (Ps0LD)
476               *-----------------*/
477              if ( Chi < large_arg ) { /* avoid exp_Chi to become extremely large */
478 	        exp_Chi = exp( Chi ) ;
479 	        fs01 = cfs1 * ( exp_Chi - 1.0e0 ) ;
480 	        fs01_dPs0 = cfs1 * beta * ( exp_Chi ) ;
481              } else {
482                 exp_bPs0 = exp( beta*Ps0LD ) ;
483                 fs01     = cnst1over * ( exp_bPs0 - exp_bVbs ) ;
484                 fs01_dPs0 = cnst1over * beta * exp_bPs0 ;
485              }
486              fs02 = sqrt( Chi - 1.0 + fs01 ) ;
487              fs02_dPs0 = ( beta + fs01_dPs0 ) / fs02 * 0.5 ;
488 
489             } /* end of if ( Chi  ... ) block */
490             /*-----------------------------------------------------------*
491              * Fs0
492              *-----------------*/
493             Fs0 = VgpLD - Ps0LD - fac1 * fs02 ;
494             Fs0_dPs0 = - 1.0e0 - fac1 * fs02_dPs0 ;
495 
496             if ( flg_conv == 1 ) break ;
497 
498             dPs0 = - Fs0 / Fs0_dPs0 ;
499 
500             /*-------------------------------------------*
501              * Update Ps0LD .
502              *-----------------*/
503             dPlim = 0.5*dP_max*(1.0 + Fn_Max(1.e0,fabs(Ps0LD))) ;
504             if ( fabs( dPs0 ) > dPlim ) dPs0 = dPlim * Fn_Sgn( dPs0 ) ;
505 
506             Ps0LD = Ps0LD + dPs0 ;
507 
508             TX = -Vxbgmtcl + ps_conv / 2 ;
509             if ( Ps0LD < TX ) Ps0LD = TX ;
510 
511             /*-------------------------------------------*
512              * Check convergence.
513              *-----------------*/
514             if ( fabs( dPs0 ) <= ps_conv && fabs( Fs0 ) <= gs_conv ) {
515               flg_conv = 1 ;
516             }
517 
518         } /* end of Ps0LD Newton loop */
519 
520         /*-------------------------------------------*
521          * Procedure for diverged case.
522          *-----------------*/
523         if ( flg_conv == 0 ) {
524           fprintf( stderr ,
525                    "*** warning(HiSIM_HV(%s)): Went Over Iteration Maximum (Ps0LD)\n",model->HSMHV2modName ) ;
526           fprintf( stderr , " -Vxbgmtcl = %e   Vgbgmt = %e\n" , -Vxbgmtcl , Vgbgmt ) ;
527         }
528 
529         /*---------------------------------------------------*
530          * Evaluate derivatives of Ps0LD.
531          *-----------------*/
532         Chi_dT = beta_dT * ( Ps0LD + Vxbgmtcl ) + beta * Vxbgmtcl_dT;
533         exp_bVbs_dT = - ( beta_dT * Vxbgmtcl + beta * Vxbgmtcl_dT ) * exp_bVbs;
534         cfs1_dT = exp_bVbs * cnst1over_dT + exp_bVbs_dT * cnst1over;
535 
536         if ( Chi < znbd5 ) {
537           fs01_dVbs = cfs1 * beta * fi * ( - fi + 2 * fi_dChi ) ; /* fs01_dVxbgmtcl */
538           fs01_dT = cfs1 * 2 * fi * fi_dChi * Chi_dT + fi * fi * cfs1_dT ;
539           T2 = 1.0e0 / ( fs02 + fs02 ) ;
540           fs02_dVbs = ( + beta * fb_dChi * 2 * fb + fs01_dVbs ) * T2 ; /* fs02_dVxbgmtcl */
541           fs02_dT = ( 2 * fb * fb_dChi * Chi_dT + fs01_dT ) * T2 ;
542         } else {
543           if ( Chi < large_arg ) {
544             fs01_dVbs = + cfs1 * beta ; /* fs01_dVxbgmtcl */
545             exp_Chi_dT  = exp_Chi * Chi_dT ;
546             fs01_dT     = ( exp_Chi - 1.0e0 ) * cfs1_dT + cfs1 * exp_Chi_dT ;
547           } else {
548             fs01_dVbs   = + cfs1 * beta ;
549             exp_bPs0_dT = exp_bPs0 * Ps0LD * beta_dT ;
550             fs01_dT     = cnst1over_dT*(exp_bPs0-exp_bVbs) + cnst1over*(exp_bPs0_dT-exp_bVbs_dT) ;
551           }
552           T2 = 0.5e0 / fs02 ;
553           fs02_dVbs = ( + beta + fs01_dVbs ) * T2 ; /* fs02_dVxbgmtcl */
554           fs02_dT = T2 * ( Chi_dT + fs01_dT ) ;
555         }
556 
557         T1 = 1.0 / Fs0_dPs0 ;
558         Ps0LD_dVxb = - ( - fac1 * fs02_dVbs ) * T1 ;
559         Ps0LD_dVds = 0.0 ;
560         Ps0LD_dVgb = - ( VgpLD_dVgb - fac1_dVgs * fs02 ) * T1 ;
561         Ps0LD_dT = - ( - ( fac1 * fs02_dT + fac1_dT * fs02 ) ) * T1;
562 
563         Chi_dT = beta_dT * ( Ps0LD + Vxbgmtcl ) + beta * ( Ps0LD_dT + Vxbgmtcl_dT );
564 
565         if ( Chi < znbd5 ) {
566           /*-------------------------------------------*
567            * zone-D1/D2. (Ps0LD)
568            *-----------------*/
569           if ( Chi < znbd3 ) { flg_ovzone = 1; }
570                         else { flg_ovzone = 2; }
571 
572           Xi0 = fb * fb + epsm10 ;
573           T1 = 2 * fb * fb_dChi * beta ;
574           Xi0_dVbs = T1 * ( Ps0LD_dVxb + 1.0 ) ; /* Xi0_dVxbgmtcl */
575           Xi0_dVds = T1 * Ps0LD_dVds ;
576           Xi0_dVgs = T1 * Ps0LD_dVgb ;
577           Xi0_dT = 2 * fb * fb_dChi * Chi_dT ;
578 
579           Xi0p12 = fb + epsm10 ;
580           T1 = fb_dChi * beta ;
581           Xi0p12_dVbs = T1 * ( Ps0LD_dVxb + 1.0 ) ; /* Xi0p12_dVxbgmtcl */
582           Xi0p12_dVds = T1 * Ps0LD_dVds ;
583           Xi0p12_dVgs = T1 * Ps0LD_dVgb ;
584           Xi0p12_dT = fb_dChi * Chi_dT ;
585 
586           Xi0p32 = fb * fb * fb + epsm10 ;
587           T1 = 3 * fb * fb * fb_dChi * beta ;
588           Xi0p32_dVbs = T1 * ( Ps0LD_dVxb + 1.0 ) ; /* Xi0p32_dVxbgmtcl */
589           Xi0p32_dVds = T1 * Ps0LD_dVds ;
590           Xi0p32_dVgs = T1 * Ps0LD_dVgb ;
591           Xi0p32_dT = 3 * fb * fb * fb_dChi * Chi_dT ;
592 
593         } else {
594           /*-------------------------------------------*
595            * zone-D3. (Ps0LD)
596            *-----------------*/
597           flg_ovzone = 3 ;
598 
599           Xi0 = Chi - 1.0e0 ;
600           Xi0_dVbs = beta * ( Ps0LD_dVxb + 1.0e0 ) ; /* Xi0_dVxbgmtcl */
601           Xi0_dVds = beta * Ps0LD_dVds ;
602           Xi0_dVgs = beta * Ps0LD_dVgb ;
603           Xi0_dT = Chi_dT ;
604 
605           Xi0p12 = sqrt( Xi0 ) ;
606           T1 = 0.5e0 / Xi0p12 ;
607           Xi0p12_dVbs = T1 * Xi0_dVbs ;
608           Xi0p12_dVds = T1 * Xi0_dVds ;
609           Xi0p12_dVgs = T1 * Xi0_dVgs ;
610           Xi0p12_dT = T1 * Xi0_dT ;
611 
612           Xi0p32 = Xi0 * Xi0p12 ;
613           T1 = 1.5e0 * Xi0p12 ;
614           Xi0p32_dVbs = T1 * Xi0_dVbs ;
615           Xi0p32_dVds = T1 * Xi0_dVds ;
616           Xi0p32_dVgs = T1 * Xi0_dVgs ;
617           Xi0p32_dT = T1 * Xi0_dT ;
618 
619         } /* end of if ( Chi  ... ) block */
620 
621         /*-----------------------------------------------------------*
622          * - Recalculate the derivatives of fs01 and fs02.
623          *-----------------*/
624         fs01_dVbs = Ps0LD_dVxb * fs01_dPs0 + fs01_dVbs ;
625         fs01_dVds = Ps0LD_dVds * fs01_dPs0 ;
626         fs01_dVgs = Ps0LD_dVgb * fs01_dPs0 ;
627         fs01_dT   = Ps0LD_dT * fs01_dPs0 + fs01_dT;
628         fs02_dVbs = Ps0LD_dVxb * fs02_dPs0 + fs02_dVbs ;
629         fs02_dVxb = Ps0LD_dVds * fs02_dPs0 ;
630         fs02_dVgb = Ps0LD_dVgb * fs02_dPs0 ;
631         fs02_dT   = Ps0LD_dT * fs02_dPs0 + fs02_dT;
632 
633         /*-----------------------------------------------------------*
634          * QbuLD and QiuLD
635          *-----------------*/
636         QbuLD = cnst0over_func * Xi0p12 ;
637         QbuLD_dVxb = cnst0over_func * Xi0p12_dVbs ;
638         QbuLD_dVgb = cnst0over_func * Xi0p12_dVgs ;
639         QbuLD_dT =   cnst0over_func * Xi0p12_dT + cnst0over_func_dT * Xi0p12;
640 
641         T1 = 1.0 / ( fs02 + Xi0p12 ) ;
642         QiuLD = cnst0over_func * fs01 * T1 ;
643         T2 = 1.0 / ( fs01 + epsm10 ) ;
644         QiuLD_dVbs = QiuLD * ( fs01_dVbs * T2 - ( fs02_dVbs + Xi0p12_dVbs ) * T1 ) ;
645         QiuLD_dVgs = QiuLD * ( fs01_dVgs * T2 - ( fs02_dVgb + Xi0p12_dVgs ) * T1 ) ;
646         T1_dT = - T1 * T1 * ( fs02_dT + Xi0p12_dT );
647         QiuLD_dT = cnst0over_func * ( fs01 * T1_dT + T1 * fs01_dT ) + fs01 * T1 * cnst0over_func_dT;
648 
649         /*-----------------------------------------------------------*
650          * Extrapolation: X_dVxbgmt = X_dVxbgmtcl * Vxbgmtcl_dVxbgmt
651          *-----------------*/
652         QbuLD_dVxb *= Vxbgmtcl_dVxbgmt ;
653         QiuLD_dVbs *= Vxbgmtcl_dVxbgmt ;
654 
655         /*-----------------------------------------------------------*
656          * Total overlap charge
657          *-----------------*/
658         QsuLD = QbuLD + QiuLD;
659         QsuLD_dVxb = QbuLD_dVxb + QiuLD_dVbs;
660         QsuLD_dVgb = QbuLD_dVgb + QiuLD_dVgs;
661         QsuLD_dT =   QbuLD_dT   + QiuLD_dT;
662 
663       } /* end of COQOVSM branches */
664 
665     } /* end of Vgbgmt region blocks */
666 
667     /* convert to source ref. */
668     Ps0LD_dVbs = Ps0LD_dVxb * Vxbgmt_dVbs + Ps0LD_dVgb * Vgbgmt_dVbs ;
669     Ps0LD_dVds = Ps0LD_dVxb * Vxbgmt_dVds + Ps0LD_dVgb * Vgbgmt_dVds ;
670     Ps0LD_dVgs = Ps0LD_dVxb * Vxbgmt_dVgs + Ps0LD_dVgb * Vgbgmt_dVgs ;
671 
672     QsuLD_dVbs = QsuLD_dVxb * Vxbgmt_dVbs + QsuLD_dVgb * Vgbgmt_dVbs ;
673     QsuLD_dVds = QsuLD_dVxb * Vxbgmt_dVds + QsuLD_dVgb * Vgbgmt_dVds ;
674     QsuLD_dVgs = QsuLD_dVxb * Vxbgmt_dVgs + QsuLD_dVgb * Vgbgmt_dVgs ;
675 
676     QbuLD_dVbs = QbuLD_dVxb * Vxbgmt_dVbs + QbuLD_dVgb * Vgbgmt_dVbs ;
677     QbuLD_dVds = QbuLD_dVxb * Vxbgmt_dVds + QbuLD_dVgb * Vgbgmt_dVds ;
678     QbuLD_dVgs = QbuLD_dVxb * Vxbgmt_dVgs + QbuLD_dVgb * Vgbgmt_dVgs ;
679 
680     /* inversion charge = total - depletion */
681     QiuLD = QsuLD - QbuLD  ;
682     QiuLD_dVbs = QsuLD_dVbs - QbuLD_dVbs ;
683     QiuLD_dVds = QsuLD_dVds - QbuLD_dVds ;
684     QiuLD_dVgs = QsuLD_dVgs - QbuLD_dVgs ;
685     QiuLD_dT = QsuLD_dT - QbuLD_dT ;
686 
687 /*  End HSMHV2evalQover */
688