Implemented the final fitting to perform extrapolation at 10 years, when the model behavior is periodic, by using Fourier basis

src/spicelib/devices/bsim4/b4reliability.c

@@ -8,6 +8,7 @@ Author: 2015 Francesco Lannutti - July 2015
 #include "ngspice/sperror.h"
 
 #include <gsl/gsl_fit.h>
+#include <gsl/gsl_linalg.h>
 
 static int
 BSIM4reliability_internal (BSIM4instance *here, CKTcircuit *ckt, unsigned int mode)
@@ -130,6 +131,9 @@ BSIM4reliability_internal (BSIM4instance *here, CKTcircuit *ckt, unsigned int mo
                 {
                     return (E_INTERN) ;
                 }
+
+                // Update the semiperiod counter
+                here->relStruct->semiPeriods++ ;
             } else {
                 fprintf (stderr, "Reliability Analysis Error\n") ;
             }
@@ -149,6 +153,9 @@ BSIM4reliability_internal (BSIM4instance *here, CKTcircuit *ckt, unsigned int mo
                 {
                     return (E_INTERN) ;
                 }
+
+                // Update the semiperiod counter
+                here->relStruct->semiPeriods++ ;
             } else if (here->relStruct->IsON == 0) {
                 // Until now, the device was OFF - Do NOTHING
                 delta = -1 ;
@@ -173,6 +180,9 @@ BSIM4reliability_internal (BSIM4instance *here, CKTcircuit *ckt, unsigned int mo
             {
                 return (E_INTERN) ;
             }
+
+            // Update the semiperiod counter
+            here->relStruct->semiPeriods++ ;
         } else if (here->relStruct->IsON == 0) {
             // Calculate recovery
             delta = ckt->CKTtime - here->relStruct->time ;
@@ -187,6 +197,9 @@ BSIM4reliability_internal (BSIM4instance *here, CKTcircuit *ckt, unsigned int mo
             {
                 return (E_INTERN) ;
             }
+
+            // Update the semiperiod counter
+            here->relStruct->semiPeriods++ ;
         } else {
             fprintf (stderr, "Reliability Analysis Error\n") ;
         }
@@ -197,7 +210,6 @@ BSIM4reliability_internal (BSIM4instance *here, CKTcircuit *ckt, unsigned int mo
 
 
         /* Calculate fitting */
-        double c0, c1, cov00, cov01, cov11, chisq ;
 
         /* Count how many deltaVth we have */
         unsigned int i ;
@@ -222,22 +234,125 @@ BSIM4reliability_internal (BSIM4instance *here, CKTcircuit *ckt, unsigned int mo
         while (current != NULL)
         {
             timeFit [i] = current->time ;
-            deltaVthFit [i] = exp (current->deltaVth) ;
+            deltaVthFit [i] = current->deltaVth ;
-            printf ("Time: %-.9g\n", current->time) ;
-            printf ("DeltaVth: %-.9g\n", current->deltaVth) ;
-            printf ("DeltaVth (exp): %-.9g\n", exp (current->deltaVth)) ;
             i++ ;
             current = current->next ;
         }
 
-        /* Execute Fitting */
+//        printf ("Number of Semi Periods: %u\n", here->relStruct->semiPeriods) ;
-        gsl_fit_linear (timeFit + 1, 1, deltaVthFit + 1, 1, i - 1, &c0, &c1, &cov00, &cov01, &cov11, &chisq) ;
 
-        /* Perform Extrapolation */
+        if (here->relStruct->semiPeriods > 1)
-        double yFit, yFitErr ;
+        {
-        gsl_fit_linear_est (2e-6, c0, c1, cov00, cov01, cov11, &yFit, &yFitErr) ;
+            /* The model behavior is periodic - Use Fourier basis fitting */
+            /* Generate the fitting matrix */
 
-        printf ("\n\tFitting -> %-.9gs\t%-.9gmV\n\n", 31536000.0, log (yFit) * 1000) ;
+            double f, factor_for_2pi, *fitting_matrix, target ;
+            unsigned int columns, j, number_of_modes, number_of_periods, rows, size ;
+
+            number_of_periods = here->relStruct->semiPeriods / 2 ;
+            number_of_modes = 10 * (number_of_periods + 1) ;
+//            printf ("Number of Periods: %u\n", number_of_periods) ;
+//            printf ("Number of Modes: %u\n", number_of_modes) ;
+
+            rows = i ;
+            columns = 2 * number_of_modes + 1 ;
+            size = rows * columns ;
+            fitting_matrix = TMALLOC (double, size) ;
+//            printf ("Rows: %u\n", rows) ;
+//            printf ("Columns: %u\n", columns) ;
+
+            factor_for_2pi = 2 * 3.14159265359 / (timeFit [rows - 1] - timeFit [0]) ;
+/*            for (i = 0 ; i < rows ; i++)
+            {
+                printf ("Time: %-.9g\n", timeFit [i]) ;
+            }
+            printf ("Max Time: %-.9g\n", timeFit [rows - 1]) ;
+            printf ("Min Time: %-.9g\n", timeFit [0]) ;
+*/
+            for (i = 0 ; i < rows ; i++)
+            {
+                /* The first element of every row is equal to 1 */
+                fitting_matrix [columns * i] = 1 ;
+
+                /* The odd elements of every row are cos(x) */
+                for (j = 0 ; j < number_of_modes ; j++)
+                {
+                    fitting_matrix [columns * i + 2 * j + 1] = cos ((j + 1) * timeFit [i] * factor_for_2pi) ;
+                }
+
+                /* The even elements of every row are sin(x) */
+                for (j = 1 ; j <= number_of_modes ; j++)
+                {
+                    fitting_matrix [columns * i + 2 * j] = sin (j * timeFit [i] * factor_for_2pi) ;
+                }
+            }
+
+            gsl_matrix_view m = gsl_matrix_view_array (fitting_matrix, rows, columns) ;
+            gsl_vector_view b = gsl_vector_view_array (deltaVthFit, rows) ;
+            gsl_vector *tau = gsl_vector_alloc (MIN (rows, columns)) ;
+
+            gsl_vector *x = gsl_vector_alloc (columns) ;
+            gsl_vector *residual = gsl_vector_alloc (rows) ;
+
+            gsl_linalg_QR_decomp (&m.matrix, tau) ;
+            gsl_linalg_QR_lssolve (&m.matrix, tau, &b.vector, x, residual) ;
+
+//            printf ("fitting_matrix = \n") ;
+//            gsl_matrix_fprintf (stdout, &m.matrix, "%g") ;
+//            printf ("x = \n") ;
+//            gsl_vector_fprintf (stdout, x, "%g") ;
+
+            target = 315360000.0 ;
+            f = gsl_vector_get (x, 0) ;
+
+            /* The odd elements of every row are cos(x) */
+            for (j = 0 ; j < number_of_modes ; j++)
+            {
+                f += gsl_vector_get (x, 2 * j + 1) * cos ((j + 1) * target * factor_for_2pi) ;
+            }
+
+            /* The even elements of every row are sin(x) */
+            for (j = 1 ; j <= number_of_modes ; j++)
+            {
+                f += gsl_vector_get (x, 2 * j) * sin (j * target * factor_for_2pi) ;
+            }
+
+            printf ("\n\nExtrapolation at 10 years:\n\t\t\t\tDeltaVth: %-.9gmV\n\n\n\n", f * 1000) ;
+
+            /* Assign the extrapolated DeltaVth to the model */
+            here->relStruct->deltaVth = f ;
+
+            gsl_vector_free (tau) ;
+            gsl_vector_free (x) ;
+            gsl_vector_free (residual) ;
+        } else {
+            /* Execute Fitting */
+            double c0, c1, cov00, cov01, cov11, chisq, *deltaVthFitExp ;
+            unsigned int j ;
+
+            deltaVthFitExp = TMALLOC (double, i) ;
+            for (j = 0 ; j < i ; j++)
+            {
+                deltaVthFitExp [j] = exp (deltaVthFit [j]) ;
+            }
+
+            gsl_fit_linear (timeFit, 1, deltaVthFitExp, 1, i, &c0, &c1, &cov00, &cov01, &cov11, &chisq) ;
+            printf ("\n\n") ;
+            printf ("c0: %-.9g\n", c0) ;
+            printf ("c1: %-.9g\n", c1) ;
+            printf ("cov00: %-.9g\n", cov00) ;
+            printf ("cov01: %-.9g\n", cov01) ;
+            printf ("cov11: %-.9g\n", cov11) ;
+            printf ("chisq: %-.9g\n", chisq) ;
+            printf ("\n\n") ;
+
+            /* Perform Extrapolation */
+            double yFit, yFitErr ;
+            gsl_fit_linear_est (315360000.0, c0, c1, cov00, cov01, cov11, &yFit, &yFitErr) ;
+
+            printf ("\n\tFitting -> %-.9gs\t%-.9gmV\n\n", 315360000.0, log (yFit) * 1000) ;
+            printf ("\n\tFitting -> %-.9gs\t%-.9gmV\n\n", 315360000.0, log (c0 + c1 * 315360000.0) * 1000) ;
+        }
     } else {
         fprintf (stderr, "Reliability Analysis Error\n") ;
     }

src/spicelib/devices/bsim4/b4set.c

@@ -2660,10 +2660,13 @@ do { if((here->ptr = SMPmakeElt(matrix,here->first,here->second))==(double *)NUL
             {
                 here->relStruct = TMALLOC (RELMODELrelStruct, 1) ;
                 here->relStruct->time = 0 ;
+                here->relStruct->offsetTime = 0 ;
                 here->relStruct->deltaVth = 0 ;
+                here->relStruct->deltaVthMax = 0 ;
                 here->relStruct->t_star = 0 ;
                 here->relStruct->IsON = -1 ;
                 here->relStruct->deltaVthList = NULL ;
+                here->relStruct->semiPeriods = 0 ;
             }
 #endif
 

src/spicelib/devices/relmodel/relmodelcalcaging.c

@@ -49,7 +49,7 @@ listInsert (RELMODELrelList **list, double time, double deltaVth)
 int
 RELMODELcalculateAging (GENinstance *inInstance, int modType, double t_aging, unsigned int stress_or_recovery)
 {
-    double A, Nt ;
+    double A, i, Nt, step ;
     BSIM4instance *here ;
     RELMODELmodel *relmodel ;
 
@@ -71,32 +71,44 @@ RELMODELcalculateAging (GENinstance *inInstance, int modType, double t_aging, un
     Nt = pow ((sqrt (relmodel->RELMODELnts)), 3) * 1e-21 ;
     A = (CHARGE / (4 * CONSTepsZero * 1e-9 * relmodel->RELMODELeps_hk)) * pow ((relmodel->RELMODELh_cut / (2 * sqrt (2 * relmodel->RELMODELm_star * relmodel->RELMODELw)) * 1e9), 2) ;
 
-    if (stress_or_recovery)
+//    printf ("\n\nStart Aging...\n") ;
+//    printf ("\tEnd step: %-.9g\n\n", here->relStruct->t_star + t_aging) ;
+    step = 1e-12 ;
+    for (i = 0 ; i < t_aging ; i += step)
     {
-        if (relmodel->RELMODELh_cut / (2 * sqrt (2 * relmodel->RELMODELm_star * relmodel->RELMODELw)) * log (1 + pow (((t_aging + here->relStruct->t_star) / relmodel->RELMODELtau_0), relmodel->RELMODELbeta)) * 1e9 <= 2)
+        if (stress_or_recovery)
         {
-            here->relStruct->deltaVth = Nt * A * pow (log (1 + pow (((t_aging + here->relStruct->t_star) / relmodel->RELMODELtau_0), relmodel->RELMODELbeta)), 2) ;
+            if (relmodel->RELMODELh_cut / (2 * sqrt (2 * relmodel->RELMODELm_star * relmodel->RELMODELw)) * log (1 + pow (((i + here->relStruct->t_star) / relmodel->RELMODELtau_0), relmodel->RELMODELbeta)) * 1e9 <= 2)
+            {
+                here->relStruct->deltaVth = Nt * A * pow (log (1 + pow (((i + here->relStruct->t_star) / relmodel->RELMODELtau_0), relmodel->RELMODELbeta)), 2) ;
+            } else {
+                here->relStruct->deltaVth = pow ((CHARGE / (4 * CONSTepsZero * 1e-9 * relmodel->RELMODELeps_hk)) * Nt * here->BSIM4modPtr->BSIM4toxe * 1e9, 2) ;
+            }
+            here->relStruct->deltaVthMax = here->relStruct->deltaVth ;
         } else {
-            here->relStruct->deltaVth = pow ((CHARGE / (4 * CONSTepsZero * 1e-9 * relmodel->RELMODELeps_hk)) * Nt * here->BSIM4modPtr->BSIM4toxe * 1e9, 2) ;
+//            printf ("\n\nDeltaVth Prior Recovery: %-.9gmV\n\n", here->relStruct->deltaVth * 1000) ;
+            here->relStruct->deltaVth = here->relStruct->deltaVthMax * log (1 + (1.718 / (1 + pow ((i / relmodel->RELMODELtau_e), relmodel->RELMODELbeta1)))) ;
         }
-    } else {
-        here->relStruct->deltaVth = here->relStruct->deltaVth * log (1 + (1.718 / (1 + pow ((t_aging / relmodel->RELMODELtau_e), relmodel->RELMODELbeta1)))) ;
-    }
 
-    /* Insert 'here->relStruct->deltaVth' into the list for the later fitting */
+        /* Insert 'here->relStruct->deltaVth' into the list for the later fitting */
-    listInsert (&(here->relStruct->deltaVthList), here->relStruct->time, here->relStruct->deltaVth) ;
+        listInsert (&(here->relStruct->deltaVthList), here->relStruct->offsetTime + i, here->relStruct->deltaVth) ;
-//    RELMODELrelList *temp ;
+//        printf ("\nStep: %-.9g\tDeltaVth: %-.9gmV\n\n", i, here->relStruct->deltaVth * 1000) ;
-//    temp = TMALLOC (RELMODELrelList, 1) ;
+//        RELMODELrelList *temp ;
-//    temp->time = here->relStruct->time ;
+//        temp = TMALLOC (RELMODELrelList, 1) ;
-//    temp->deltaVth = here->relStruct->deltaVth ;
+//        temp->time = i + here->relStruct->t_star ;
-//    temp->next = here->relStruct->deltaVthList ;
+//        temp->deltaVth = here->relStruct->deltaVth ;
-//    here->relStruct->deltaVthList = temp ;
+//        temp->next = here->relStruct->deltaVthList ;
-//    printf ("QUI\n\n") ;
+//        here->relStruct->deltaVthList = temp ;
+//        printf ("QUI\n\n") ;
+    }
+//    printf ("Stop Aging...\n\n") ;
 
     if (!stress_or_recovery)
     {
         here->relStruct->t_star = pow ((exp (sqrt (here->relStruct->deltaVth / (Nt * A))) - 1), (1 / relmodel->RELMODELbeta)) * relmodel->RELMODELtau_0 ;
     }
 
+    here->relStruct->offsetTime += i ;
+
     return 0 ;
 }

src/spicelib/devices/relmodel/relmodeldefs.h

@@ -19,10 +19,13 @@ typedef struct sRELMODELrelList {
 
 typedef struct sRELMODELrelStruct {
     double time ;
+    double offsetTime ;
     double deltaVth ;
+    double deltaVthMax ;
     double t_star ;
     int IsON ;
     RELMODELrelList *deltaVthList ;
+    unsigned int semiPeriods ;
 } RELMODELrelStruct ;
 
 typedef struct sRELMODELmodel