diff --git a/src/spicelib/devices/dio/dioload.c b/src/spicelib/devices/dio/dioload.c
index 63fd16108..dcc1a2554 100644
--- a/src/spicelib/devices/dio/dioload.c
+++ b/src/spicelib/devices/dio/dioload.c
@@ -2,7 +2,7 @@
 Copyright 1990 Regents of the University of California.  All rights reserved.
 Author: 1985 Thomas L. Quarles
 Modified: 2000 AlansFixes
-Modified by Paolo Nenzi 2003 and Dietmar Warning 2012
+Modified by Paolo Nenzi 2003, Dietmar Warning 2012 and Ste Kulov 2021
 **********/
 
 #include "ngspice/ngspice.h"
@@ -31,7 +31,6 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
     double csat;    /* area-scaled saturation current */
     double czero;
     double czof2;
-
     double delvd;   /* change in diode voltage temporary */
     double evd;
     double evrev;
@@ -50,6 +49,7 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
     int error;
     int SenCond=0;    /* sensitivity condition */
     double diffcharge, deplcharge, diffcap, deplcap;
+
     double tt;
     double vp;
 
@@ -81,6 +81,7 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
             vtebrk = model->DIObrkdEmissionCoeff * vt;
             tt = here->DIOtTransitTime;
             vp = model->DIOsoftRevRecParam;
+
             /*
              *   initialization
              */
@@ -151,10 +152,10 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
                                 ckt->CKTabstol;
                         if (fabs(cdhat- *(ckt->CKTstate0 + here->DIOcurrent))
                                 < tol) {
-                                vd= *(ckt->CKTstate0 + here->DIOvoltage);
-                                cd= *(ckt->CKTstate0 + here->DIOcurrent);
-                                gd= *(ckt->CKTstate0 + here->DIOconduct);
-                                goto load;
+                            vd= *(ckt->CKTstate0 + here->DIOvoltage);
+                            cd= *(ckt->CKTstate0 + here->DIOcurrent);
+                            gd= *(ckt->CKTstate0 + here->DIOconduct);
+                            goto load;
                         }
                     }
                 }
@@ -178,52 +179,34 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
             /*
              *   compute dc current and derivitives
              */
-next1:      
-
-            if (vd >= -3*vte) {
+next1:
+            if (vd >= -3*vte) {                 /* forward */
 
                 evd = exp(vd/vte);
-                cd = csat*(evd-1);
-                gd = csat*evd/vte;
+                cd = csat*(evd-1) + ckt->CKTgmin*vd;
+                gd = csat*evd/vte + ckt->CKTgmin;
 
             } else if((!(model->DIObreakdownVoltageGiven)) ||
-                    vd >= -here->DIOtBrkdwnV) {
+                    vd >= -here->DIOtBrkdwnV) { /* reverse */
 
                 arg = 3*vte/(vd*CONSTe);
                 arg = arg * arg * arg;
-                cd = -csat*(1+arg);
-                gd = csat*3*arg/vd;
+                cd = -csat*(1+arg) + ckt->CKTgmin*vd;
+                gd = csat*3*arg/vd + ckt->CKTgmin;
 
-            } else {
+            } else {                            /* breakdown */
 
                 evrev = exp(-(here->DIOtBrkdwnV+vd)/vtebrk);
-                cd = -csat*evrev;
-                gd = csat*evrev/vtebrk;
+                cd = -csat*evrev + ckt->CKTgmin*vd;
+                gd = csat*evrev/vtebrk + ckt->CKTgmin;
 
             }
 
-            /*
-            if (ckt->CKTmode & MODEINITTRAN) {
-                *(ckt->CKTstate1 + here->DIOoldCurr) = cd;
-                *(ckt->CKTstate1 + here->DIOoldCond) = gd;
-                *(ckt->CKTstate2 + here->DIOoldCurr) = cd;
-                *(ckt->CKTstate2 + here->DIOoldCond) = gd;
-            }
-            else {
-                //cd = (*(ckt->CKTstate1 + here->DIOoldCurr) - (ckt->CKTdelta)*cd) / (1 + (ckt->CKTdelta)*vp);
-                //gd = (*(ckt->CKTstate1 + here->DIOoldCond) - (ckt->CKTdelta)*gd) / (1 + (ckt->CKTdelta)*vp);
-                //cd = (1 - vp)*cd + vp*(*(ckt->CKTstate1 + here->DIOoldCurr))*(ckt->CKTdelta);
-                //gd = (1 - vp)*gd + vp*(*(ckt->CKTstate1 + here->DIOoldCond))*(ckt->CKTdelta);
-            }
-            */
-
             if ((ckt->CKTmode & (MODEDCTRANCURVE | MODETRAN | MODEAC | MODEINITSMSIG)) ||
                      ((ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC))) {
               /*
                *   charge storage elements
                */
-
-                /* junction(depletion) charge */
                 czero=here->DIOtJctCap;
                 if (vd < here->DIOtDepCap){
                     arg=1-vd/here->DIOtJctPot;
@@ -237,105 +220,55 @@ next1:
                     deplcap = czof2*(here->DIOtF3+here->DIOtGradingCoeff*vd/here->DIOtJctPot);
                 }
 
-                /* diffusion charge */
-                if (ckt->CKTmode & MODEINITTRAN) {
-                    diffcharge = tt * cd;
-                    diffcap = tt * gd;
-                    *(ckt->CKTstate2 + here->DIOdiffCharge) = diffcharge;
-                    *(ckt->CKTstate2 + here->DIOdiffCap) = diffcap;
-                    *(ckt->CKTstate1 + here->DIOdiffCharge) = diffcharge;
-                    *(ckt->CKTstate1 + here->DIOdiffCap) = diffcap;
-                    *(ckt->CKTstate1 + here->DIOoldCurr) = cd;
-                    *(ckt->CKTstate1 + here->DIOoldCond) = gd;
-                }
-                else {
-                    double dt = ckt->CKTdelta;
-                    /*
-                    //Fixed-point iteration
-                    double diffcharge_old, diffcap_old;
+                if (model->DIOsoftRevRecParamGiven) {
 
-                    //Forward Euler predictor
-                    diffcharge = (*(ckt->CKTstate1 + here->DIOdiffCharge))*(1 - dt / (tt*vp)) + (dt*(*(ckt->CKTstate1 + here->DIOoldCurr)) / vp);
-                    diffcap = (*(ckt->CKTstate1 + here->DIOdiffCap))*(1 - dt / (tt*vp)) + (dt*(*(ckt->CKTstate1 + here->DIOoldCond)) / vp);
-                    do {
-                        diffcharge_old = diffcharge;
-                        diffcharge = (*(ckt->CKTstate1 + here->DIOdiffCharge)) + ((dt / vp)*(cd - (diffcharge_old / tt)));
-                    } while (fabs(diffcharge - diffcharge_old) > 1e-12);
-                    do {
-                        diffcap_old = diffcap;
-                        diffcap = (*(ckt->CKTstate1 + here->DIOdiffCap)) + ((dt / vp)*(gd - (diffcap_old / tt)));
-                    } while (fabs(diffcap - diffcap_old) > 1e-12);
-                    */
+                    if (ckt->CKTmode & MODEINITTRAN) {
+                        diffcharge = tt * cd;
+                        diffcap = tt * gd;
+                        *(ckt->CKTstate2 + here->DIOdiffCharge) = diffcharge;
+                        *(ckt->CKTstate2 + here->DIOdiffCap) = diffcap;
+                        *(ckt->CKTstate1 + here->DIOdiffCharge) = diffcharge;
+                        *(ckt->CKTstate1 + here->DIOdiffCap) = diffcap;
+                        *(ckt->CKTstate1 + here->DIOoldCurr) = cd;
+                        *(ckt->CKTstate1 + here->DIOoldCond) = gd;
+                    }
+                    else {
+                        double dt = ckt->CKTdelta;
 
-                    /*
-                    //Newton's method
-                    double diffcharge_old, diffcap_old;
+                        //Backward Euler
+                        // Qk+1 = tt (VP Qk + h Ik+1) / (tt VP + h)
+                        //diffcharge = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCharge)) + dt*cd) / (tt*vp + dt));
+                        //diffcap = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCap)) + dt*gd) / (tt*vp + dt));
 
-                    //Forward Euler predictor
-                    diffcharge = (*(ckt->CKTstate1 + here->DIOdiffCharge))*(1 - dt / (tt*vp)) + (dt*(*(ckt->CKTstate1 + here->DIOoldCurr)) / vp);
-                    diffcap = (*(ckt->CKTstate1 + here->DIOdiffCap))*(1 - dt / (tt*vp)) + (dt*(*(ckt->CKTstate1 + here->DIOoldCond)) / vp);
-                    do {
-                        diffcharge_old = diffcharge;
-                        diffcharge = diffcharge_old + ((tt*vp*(*(ckt->CKTstate1 + here->DIOdiffCharge)) + tt*dt*cd - diffcharge_old*(dt + tt*vp)) / (dt + tt*vp));
-                    } while (fabs(diffcharge - diffcharge_old) > 1e-12);
+                        //Trap
+                        // Qk+1 = (Qk (2 tt VP - h) + h tt (Ik + Ik+1)) / (2 tt VP + h)
+                        diffcharge = ((*(ckt->CKTstate1 + here->DIOdiffCharge))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCurr)) + cd)) / (2.0*tt*vp + dt);
+                        diffcap = ((*(ckt->CKTstate1 + here->DIOdiffCap))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCond)) + gd)) / (2.0*tt*vp + dt);
 
-                    do {
-                        diffcap_old = diffcap;
-                        diffcap = diffcap_old + ((tt*vp*(*(ckt->CKTstate1 + here->DIOdiffCap)) + tt * dt*gd - diffcap_old * (dt + tt*vp)) / (dt + tt*vp));
-                    } while (fabs(diffcap - diffcap_old) > 1e-12);
-                    */
+                        //Gear
+                        // Qk+1 = tt ((VP (4 Qk - Qk-1) + 2 h Ik+1) / (3 tt VP + 2 h))
+                        //diffcharge = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCharge)) - (*(ckt->CKTstate2 + here->DIOdiffCharge))) + 2.0*dt*cd) / (3.0*tt*vp + 2.0*dt));
+                        //diffcap = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCap)) - (*(ckt->CKTstate2 + here->DIOdiffCap))) + 2.0*dt*gd) / (3.0*tt*vp + 2.0*dt));
 
-                    //Backward Difference Operator
-                    //diffcharge = tt * (cd - vp * (((*(ckt->CKTstate1 + here->DIOdiffCharge)) - (*(ckt->CKTstate2 + here->DIOdiffCharge))) / (ckt->CKTdeltaOld[1])));
-                    //diffcap = tt * (gd - vp * (((*(ckt->CKTstate1 + here->DIOdiffCap)) - (*(ckt->CKTstate2 + here->DIOdiffCap))) / (ckt->CKTdeltaOld[1])));
+                    }
 
-                    //Standard Diffy-Q w/ trap integral
-                    //diffcharge = (dt / (2.0*vt))*(exp(-dt / (tt*vt))*(*(ckt->CKTstate1 + here->DIOoldCurr)) + cd) + exp(-dt / (tt*vt))*(*(ckt->CKTstate1 + here->DIOdiffCharge));
-                    //diffcap = (dt / (2.0*vt))*(exp(-dt / (tt*vt))*(*(ckt->CKTstate1 + here->DIOoldCond)) + gd) + exp(-dt / (tt*vt))*(*(ckt->CKTstate1 + here->DIOdiffCap));
+                    //printf("time: %e vd: %e\n", ckt->CKTtime, vd);
+                    //printf("%e %e %e %e %e %e %e\n", ckt->CKTtime, ckt->CKTdelta, diffcharge, diffcap, vd, cd, gd);
 
-                    //Forward Euler
-                    //diffcharge = ((*(ckt->CKTstate1 + here->DIOdiffCharge))*(tt*vp - dt) + dt*tt*cd) / (tt*vp);
-                    //diffcap = ((*(ckt->CKTstate1 + here->DIOdiffCap))*(tt*vp - dt) + dt*tt*gd) / (tt*vp);
+                } else {
 
-                    //Backward Euler
-                    //diffcharge = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCharge)) + dt*cd) / (tt*vp + dt));
-                    //diffcap = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCap)) + dt*gd) / (tt*vp + dt));
+                    diffcharge = tt*cd;
+                    diffcap = tt*gd;
 
-                    //Trap
-                    diffcharge = ((*(ckt->CKTstate1 + here->DIOdiffCharge))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCurr)) + cd)) / (2.0*tt*vp + dt);
-                    diffcap = ((*(ckt->CKTstate1 + here->DIOdiffCap))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCond)) + gd)) / (2.0*tt*vp + dt);
-
-                    //Gear
-                    //diffcharge = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCharge)) - (*(ckt->CKTstate2 + here->DIOdiffCharge))) + 2.0*dt*cd) / (3.0*tt*vp + 2.0*dt));
-                    //diffcap = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCap)) - (*(ckt->CKTstate2 + here->DIOdiffCap))) + 2.0*dt*gd) / (3.0*tt*vp + 2.0*dt));
-
-                    //Nothing
-                    //diffcharge = tt * cd;
-                    //diffcap = tt * gd;
                 }
 
                 *(ckt->CKTstate0 + here->DIOdiffCharge) = diffcharge;
                 *(ckt->CKTstate0 + here->DIOdiffCap) = diffcap;
+                *(ckt->CKTstate0 + here->DIOcapCharge) = deplcharge + diffcharge;
                 *(ckt->CKTstate0 + here->DIOoldCurr) = cd;
                 *(ckt->CKTstate0 + here->DIOoldCond) = gd;
 
-                //printf("Time: %e\n", ckt->CKTtime);
-                //printf("Charge: %e\n", diffcharge);
-                //printf("Cap: %e\n", diffcap);
-                //printf("Curr: %e\n", cd);
-                //printf("Cond: %e\n", gd);
-                //printf("DelCharge: %e\n", diffcharge-(*(ckt->CKTstate1 + here->DIOdiffCharge)));
-                //printf("DelCap: %e\n", diffcap-(*(ckt->CKTstate1 + here->DIOdiffCap)));
-                //printf("DelCurr: %e\n", cd-(*(ckt->CKTstate1 + here->DIOoldCurr)));
-                //printf("DelCond: %e\n\n", gd-(*(ckt->CKTstate1 + here->DIOoldCond)));
-
-                //diffcharge = here->DIOtTransitTime*cd;
-                *(ckt->CKTstate0 + here->DIOcapCharge) =
-                        diffcharge + deplcharge;
-
-                //diffcap = here->DIOtTransitTime*gd;
-
-                capd = diffcap + deplcap + here->DIOcmetal + here->DIOcpoly;
+                capd = deplcap + diffcap;
 
                 here->DIOcap = capd;
 
@@ -394,7 +327,7 @@ next1:
              *   check convergence
              */
             if ( (!(ckt->CKTmode & MODEINITFIX)) || (!(here->DIOoff))  ) {
-                if (Check == 1) {
+                if (Check == 1)  {
                     ckt->CKTnoncon++;
                     ckt->CKTtroubleElt = (GENinstance *) here;
                 }
diff --git a/src/spicelib/devices/dio/dioload_separate.c b/src/spicelib/devices/dio/dioload_separate.c
new file mode 100644
index 000000000..723afb089
--- /dev/null
+++ b/src/spicelib/devices/dio/dioload_separate.c
@@ -0,0 +1,372 @@
+/**********
+Copyright 1990 Regents of the University of California.  All rights reserved.
+Author: 1985 Thomas L. Quarles
+Modified: 2000 AlansFixes
+Modified by Paolo Nenzi 2003, Dietmar Warning 2012 and Ste Kulov 2021
+**********/
+
+#include "ngspice/ngspice.h"
+#include "ngspice/devdefs.h"
+#include "ngspice/cktdefs.h"
+#include "diodefs.h"
+#include "ngspice/const.h"
+#include "ngspice/trandefs.h"
+#include "ngspice/sperror.h"
+#include "ngspice/suffix.h"
+
+int
+DIOload(GENmodel *inModel, CKTcircuit *ckt)
+        /* actually load the current resistance value into the
+         * sparse matrix previously provided
+         */
+{
+    DIOmodel *model = (DIOmodel*)inModel;
+    DIOinstance *here;
+    double arg;
+    double capd;
+    double cd;
+    double cdeq;
+    double cdhat;
+    double ceq;
+    double csat;    /* area-scaled saturation current */
+    double czero;
+    double czof2;
+    double delvd;   /* change in diode voltage temporary */
+    double evd;
+    double evrev;
+    double gd;
+    double geq;
+    double gspr;    /* area-scaled conductance */
+    double sarg;
+#ifndef NOBYPASS
+    double tol;     /* temporary for tolerence calculations */
+#endif
+    double vd;      /* current diode voltage */
+    double vdtemp;
+    double vt;      /* K t / Q */
+    double vte, vtebrk;
+    int Check=0;
+    int error;
+    int SenCond=0;    /* sensitivity condition */
+    double diffcharge, deplcharge, diffcap, deplcap;
+
+    double tt;
+    double vp;
+
+    /*  loop through all the diode models */
+    for( ; model != NULL; model = DIOnextModel(model)) {
+
+        /* loop through all the instances of the model */
+        for (here = DIOinstances(model); here != NULL ;
+                here=DIOnextInstance(here)) {
+
+            /*
+             *     this routine loads diodes for dc and transient analyses.
+             */
+
+            if(ckt->CKTsenInfo){
+                if((ckt->CKTsenInfo->SENstatus == PERTURBATION)
+                        && (here->DIOsenPertFlag == OFF))continue;
+                SenCond = here->DIOsenPertFlag;
+
+#ifdef SENSDEBUG
+                printf("DIOload \n");
+#endif /* SENSDEBUG */
+
+            }
+            csat = here->DIOtSatCur;
+            gspr = here->DIOtConductance;
+            vt = CONSTKoverQ * here->DIOtemp;
+            vte = model->DIOemissionCoeff * vt;
+            vtebrk = model->DIObrkdEmissionCoeff * vt;
+            tt = here->DIOtTransitTime;
+            vp = model->DIOsoftRevRecParam;
+
+            /*
+             *   initialization
+             */
+
+            if(SenCond){
+
+#ifdef SENSDEBUG
+                printf("DIOsenPertFlag = ON \n");
+#endif /* SENSDEBUG */
+
+                if((ckt->CKTsenInfo->SENmode == TRANSEN)&&
+                        (ckt->CKTmode & MODEINITTRAN)) {
+                    vd = *(ckt->CKTstate1 + here->DIOvoltage);
+                } else{
+                    vd = *(ckt->CKTstate0 + here->DIOvoltage);
+                }
+
+#ifdef SENSDEBUG
+                printf("vd = %.7e \n",vd);
+#endif /* SENSDEBUG */
+                goto next1;
+            }
+
+            Check=1;
+            if(ckt->CKTmode & MODEINITSMSIG) {
+                vd= *(ckt->CKTstate0 + here->DIOvoltage);
+            } else if (ckt->CKTmode & MODEINITTRAN) {
+                vd= *(ckt->CKTstate1 + here->DIOvoltage);
+            } else if ( (ckt->CKTmode & MODEINITJCT) &&
+                    (ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC) ) {
+                vd=here->DIOinitCond;
+            } else if ( (ckt->CKTmode & MODEINITJCT) && here->DIOoff) {
+                vd=0;
+            } else if ( ckt->CKTmode & MODEINITJCT) {
+                vd=here->DIOtVcrit;
+            } else if ( ckt->CKTmode & MODEINITFIX && here->DIOoff) {
+                vd=0;
+            } else {
+#ifndef PREDICTOR
+                if (ckt->CKTmode & MODEINITPRED) {
+                    *(ckt->CKTstate0 + here->DIOvoltage) =
+                            *(ckt->CKTstate1 + here->DIOvoltage);
+                    vd = DEVpred(ckt,here->DIOvoltage);
+                    *(ckt->CKTstate0 + here->DIOcurrent) =
+                            *(ckt->CKTstate1 + here->DIOcurrent);
+                    *(ckt->CKTstate0 + here->DIOconduct) =
+                            *(ckt->CKTstate1 + here->DIOconduct);
+                } else {
+#endif /* PREDICTOR */
+                    vd = *(ckt->CKTrhsOld+here->DIOposPrimeNode)-
+                            *(ckt->CKTrhsOld + here->DIOnegNode);
+#ifndef PREDICTOR
+                }
+#endif /* PREDICTOR */
+                delvd=vd- *(ckt->CKTstate0 + here->DIOvoltage);
+                cdhat= *(ckt->CKTstate0 + here->DIOcurrent) +
+                        *(ckt->CKTstate0 + here->DIOconduct) * delvd;
+                /*
+                 *   bypass if solution has not changed
+                 */
+#ifndef NOBYPASS
+                if ((!(ckt->CKTmode & MODEINITPRED)) && (ckt->CKTbypass)) {
+                    tol=ckt->CKTvoltTol + ckt->CKTreltol*
+                        MAX(fabs(vd),fabs(*(ckt->CKTstate0 +here->DIOvoltage)));
+                    if (fabs(delvd) < tol){
+                        tol=ckt->CKTreltol* MAX(fabs(cdhat),
+                                fabs(*(ckt->CKTstate0 + here->DIOcurrent)))+
+                                ckt->CKTabstol;
+                        if (fabs(cdhat- *(ckt->CKTstate0 + here->DIOcurrent))
+                                < tol) {
+                            vd= *(ckt->CKTstate0 + here->DIOvoltage);
+                            cd= *(ckt->CKTstate0 + here->DIOcurrent);
+                            gd= *(ckt->CKTstate0 + here->DIOconduct);
+                            goto load;
+                        }
+                    }
+                }
+#endif /* NOBYPASS */
+                /*
+                 *   limit new junction voltage
+                 */
+                if ( (model->DIObreakdownVoltageGiven) &&
+                        (vd < MIN(0,-here->DIOtBrkdwnV+10*vtebrk))) {
+                    vdtemp = -(vd+here->DIOtBrkdwnV);
+                    vdtemp = DEVpnjlim(vdtemp,
+                            -(*(ckt->CKTstate0 + here->DIOvoltage) +
+                            here->DIOtBrkdwnV),vtebrk,
+                            here->DIOtVcrit,&Check);
+                    vd = -(vdtemp+here->DIOtBrkdwnV);
+                } else {
+                    vd = DEVpnjlim(vd,*(ckt->CKTstate0 + here->DIOvoltage),
+                            vte,here->DIOtVcrit,&Check);
+                }
+            }
+            /*
+             *   compute dc current and derivitives
+             */
+next1:
+            if (vd >= -3*vte) {                 /* forward */
+
+                evd = exp(vd/vte);
+                cd = csat*(evd-1) + ckt->CKTgmin*vd;
+                gd = csat*evd/vte + ckt->CKTgmin;
+
+            } else if((!(model->DIObreakdownVoltageGiven)) ||
+                    vd >= -here->DIOtBrkdwnV) { /* reverse */
+
+                arg = 3*vte/(vd*CONSTe);
+                arg = arg * arg * arg;
+                cd = -csat*(1+arg) + ckt->CKTgmin*vd;
+                gd = csat*3*arg/vd + ckt->CKTgmin;
+
+            } else {                            /* breakdown */
+
+                evrev = exp(-(here->DIOtBrkdwnV+vd)/vtebrk);
+                cd = -csat*evrev + ckt->CKTgmin*vd;
+                gd = csat*evrev/vtebrk + ckt->CKTgmin;
+
+            }
+
+            if ((ckt->CKTmode & (MODEDCTRANCURVE | MODETRAN | MODEAC | MODEINITSMSIG)) ||
+                     ((ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC))) {
+              /*
+               *   charge storage elements
+               */
+                czero=here->DIOtJctCap;
+                if (vd < here->DIOtDepCap){
+                    arg=1-vd/here->DIOtJctPot;
+                    sarg=exp(-here->DIOtGradingCoeff*log(arg));
+                    deplcharge = here->DIOtJctPot*czero*(1-arg*sarg)/(1-here->DIOtGradingCoeff);
+                    deplcap = czero*sarg;
+                } else {
+                    czof2=czero/here->DIOtF2;
+                    deplcharge = czero*here->DIOtF1+czof2*(here->DIOtF3*(vd-here->DIOtDepCap)+
+                                 (here->DIOtGradingCoeff/(here->DIOtJctPot+here->DIOtJctPot))*(vd*vd-here->DIOtDepCap*here->DIOtDepCap));
+                    deplcap = czof2*(here->DIOtF3+here->DIOtGradingCoeff*vd/here->DIOtJctPot);
+                }
+
+                *(ckt->CKTstate0 + here->DIOcapCharge) = deplcharge;
+
+                if (model->DIOsoftRevRecParamGiven) {
+
+                    if (ckt->CKTmode & MODEINITTRAN) {
+                        diffcharge = tt * cd;
+                        diffcap = tt * gd;
+                        *(ckt->CKTstate2 + here->DIOdiffCharge) = diffcharge;
+                        *(ckt->CKTstate2 + here->DIOdiffCap) = diffcap;
+                        *(ckt->CKTstate1 + here->DIOdiffCharge) = diffcharge;
+                        *(ckt->CKTstate1 + here->DIOdiffCap) = diffcap;
+                        *(ckt->CKTstate1 + here->DIOoldCurr) = cd;
+                        *(ckt->CKTstate1 + here->DIOoldCond) = gd;
+                    }
+                    else {
+                        double dt = ckt->CKTdelta;
+
+                        //Backward Euler
+                        // Qk+1 = tt (VP Qk + h Ik+1) / (tt VP + h)
+                        //diffcharge = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCharge)) + dt*cd) / (tt*vp + dt));
+                        //diffcap = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCap)) + dt*gd) / (tt*vp + dt));
+
+                        //Trap
+                        // Qk+1 = (Qk (2 tt VP - h) + h tt (Ik + Ik+1)) / (2 tt VP + h)
+                        diffcharge = ((*(ckt->CKTstate1 + here->DIOdiffCharge))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCurr)) + cd)) / (2.0*tt*vp + dt);
+                        diffcap = ((*(ckt->CKTstate1 + here->DIOdiffCap))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCond)) + gd)) / (2.0*tt*vp + dt);
+
+                        //Gear
+                        // Qk+1 = tt ((VP (4 Qk - Qk-1) + 2 h Ik+1) / (3 tt VP + 2 h))
+                        //diffcharge = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCharge)) - (*(ckt->CKTstate2 + here->DIOdiffCharge))) + 2.0*dt*cd) / (3.0*tt*vp + 2.0*dt));
+                        //diffcap = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCap)) - (*(ckt->CKTstate2 + here->DIOdiffCap))) + 2.0*dt*gd) / (3.0*tt*vp + 2.0*dt));
+
+                    }
+
+                    //printf("time: %e vd: %e\n", ckt->CKTtime, vd);
+                    //printf("%e %e %e %e %e %e\n", ckt->CKTtime, ckt->CKTdelta, diffcharge, diffcap, cd, gd);
+
+                } else {
+
+                    diffcharge = tt*cd;
+                    diffcap = tt*gd;
+
+                }
+
+                *(ckt->CKTstate0 + here->DIOdiffCharge) = diffcharge;
+                *(ckt->CKTstate0 + here->DIOdiffCap) = diffcap;
+                //*(ckt->CKTstate0 + here->DIOoldCurr) = cd;
+                //*(ckt->CKTstate0 + here->DIOoldCond) = gd;
+
+                capd = deplcap + diffcap;
+
+                here->DIOcap = capd;
+
+                /*
+                 *   store small-signal parameters
+                 */
+                if( (!(ckt->CKTmode & MODETRANOP)) ||
+                        (!(ckt->CKTmode & MODEUIC)) ) {
+                    if (ckt->CKTmode & MODEINITSMSIG){
+                        *(ckt->CKTstate0 + here->DIOcapCurrent) = capd;
+
+                        if(SenCond){
+                            *(ckt->CKTstate0 + here->DIOcurrent) = cd;
+                            *(ckt->CKTstate0 + here->DIOconduct) = gd;
+#ifdef SENSDEBUG
+                            printf("storing small signal parameters\n");
+                            printf("cd = %.7e,vd = %.7e\n",cd,vd);
+                            printf("capd = %.7e ,gd = %.7e \n",capd,gd);
+#endif /* SENSDEBUG */
+                        }
+                        continue;
+                    }
+
+                    /*
+                     *   transient analysis
+                     */
+                    if(SenCond && (ckt->CKTsenInfo->SENmode == TRANSEN)){
+                        *(ckt->CKTstate0 + here->DIOcurrent) = cd;
+#ifdef SENSDEBUG
+                        printf("storing parameters for transient sensitivity\n"
+                                );
+                        printf("qd = %.7e, capd = %.7e,cd = %.7e\n",
+                                *(ckt->CKTstate0 + here->DIOcapCharge),capd,cd);
+#endif /* SENSDEBUG */
+                        continue;
+                    }
+
+                    if (ckt->CKTmode & MODEINITTRAN) {
+                        *(ckt->CKTstate1 + here->DIOcapCharge) =
+                                *(ckt->CKTstate0 + here->DIOcapCharge);
+                        *(ckt->CKTstate1 + here->DIOdiffCharge) =
+                                *(ckt->CKTstate0 + here->DIOdiffCharge);
+                    }
+                    error = NIintegrate(ckt,&geq,&ceq,deplcap,here->DIOcapCharge);
+                    if(error) return(error);
+                    gd=gd+geq;
+                    cd=cd+*(ckt->CKTstate0 + here->DIOcapCurrent);
+                    error = NIintegrate(ckt,&geq,&ceq,diffcap,here->DIOdiffCharge);
+                    if(error) return(error);
+                    gd=gd+geq;
+                    cd=cd+*(ckt->CKTstate0 + here->DIOdiffCurrent);
+                    if (ckt->CKTmode & MODEINITTRAN) {
+                        *(ckt->CKTstate1 + here->DIOcapCurrent) =
+                                *(ckt->CKTstate0 + here->DIOcapCurrent);
+                        *(ckt->CKTstate1 + here->DIOdiffCurrent) =
+                                *(ckt->CKTstate0 + here->DIOdiffCurrent);
+                    }
+                }
+            }
+
+            if(SenCond) goto next2;
+
+            /*
+             *   check convergence
+             */
+            if ( (!(ckt->CKTmode & MODEINITFIX)) || (!(here->DIOoff))  ) {
+                if (Check == 1)  {
+                    ckt->CKTnoncon++;
+                    ckt->CKTtroubleElt = (GENinstance *) here;
+                }
+            }
+next2:      *(ckt->CKTstate0 + here->DIOvoltage) = vd;
+            *(ckt->CKTstate0 + here->DIOcurrent) = cd;
+            *(ckt->CKTstate0 + here->DIOconduct) = gd;
+
+            if(SenCond)  continue;
+
+#ifndef NOBYPASS
+            load:
+#endif
+            /*
+             *   load current vector
+             */
+            cdeq=cd-gd*vd;
+            *(ckt->CKTrhs + here->DIOnegNode) += cdeq;
+            *(ckt->CKTrhs + here->DIOposPrimeNode) -= cdeq;
+            /*
+             *   load matrix
+             */
+            *(here->DIOposPosPtr) += gspr;
+            *(here->DIOnegNegPtr) += gd;
+            *(here->DIOposPrimePosPrimePtr) += (gd + gspr);
+            *(here->DIOposPosPrimePtr) -= gspr;
+            *(here->DIOnegPosPrimePtr) -= gd;
+            *(here->DIOposPrimePosPtr) -= gspr;
+            *(here->DIOposPrimeNegPtr) -= gd;
+        }
+    }
+    return(OK);
+}
diff --git a/src/spicelib/devices/dio/dioload_sk.c b/src/spicelib/devices/dio/dioload_sk.c
index 6692d682b..99fdc5add 100644
--- a/src/spicelib/devices/dio/dioload_sk.c
+++ b/src/spicelib/devices/dio/dioload_sk.c
@@ -2,7 +2,7 @@
 Copyright 1990 Regents of the University of California.  All rights reserved.
 Author: 1985 Thomas L. Quarles
 Modified: 2000 AlansFixes
-Modified by Paolo Nenzi 2003, Dietmar Warning 2012 and Ste Kulov 2021
+Modified by Paolo Nenzi 2003 and Dietmar Warning 2012
 **********/
 
 #include "ngspice/ngspice.h"
@@ -23,6 +23,7 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
     DIOmodel *model = (DIOmodel*)inModel;
     DIOinstance *here;
     double arg;
+    double capd;
     double cd;
     double cdeq;
     double cdhat;
@@ -30,6 +31,7 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
     double csat;    /* area-scaled saturation current */
     double czero;
     double czof2;
+
     double delvd;   /* change in diode voltage temporary */
     double evd;
     double evrev;
@@ -48,15 +50,12 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
     int error;
     int SenCond=0;    /* sensitivity condition */
     double diffcharge, deplcharge, diffcap, deplcap;
-
     double tt;
     double vp;
 
     /*  loop through all the diode models */
     for( ; model != NULL; model = DIOnextModel(model)) {
 
-        vp = model->DIOsoftRevRecParam;
-
         /* loop through all the instances of the model */
         for (here = DIOinstances(model); here != NULL ;
                 here=DIOnextInstance(here)) {
@@ -81,6 +80,7 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
             vte = model->DIOemissionCoeff * vt;
             vtebrk = model->DIObrkdEmissionCoeff * vt;
             tt = here->DIOtTransitTime;
+            vp = model->DIOsoftRevRecParam;
             /*
              *   initialization
              */
@@ -151,10 +151,10 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
                                 ckt->CKTabstol;
                         if (fabs(cdhat- *(ckt->CKTstate0 + here->DIOcurrent))
                                 < tol) {
-                            vd= *(ckt->CKTstate0 + here->DIOvoltage);
-                            cd= *(ckt->CKTstate0 + here->DIOcurrent);
-                            gd= *(ckt->CKTstate0 + here->DIOconduct);
-                            goto load;
+                                vd= *(ckt->CKTstate0 + here->DIOvoltage);
+                                cd= *(ckt->CKTstate0 + here->DIOcurrent);
+                                gd= *(ckt->CKTstate0 + here->DIOconduct);
+                                goto load;
                         }
                     }
                 }
@@ -178,22 +178,23 @@ DIOload(GENmodel *inModel, CKTcircuit *ckt)
             /*
              *   compute dc current and derivitives
              */
-next1:
-            if (vd >= -3*vte) {                 /* forward */
+next1:      
+
+            if (vd >= -3*vte) {
 
                 evd = exp(vd/vte);
                 cd = csat*(evd-1);
                 gd = csat*evd/vte;
 
             } else if((!(model->DIObreakdownVoltageGiven)) ||
-                    vd >= -here->DIOtBrkdwnV) { /* reverse */
+                    vd >= -here->DIOtBrkdwnV) {
 
                 arg = 3*vte/(vd*CONSTe);
                 arg = arg * arg * arg;
                 cd = -csat*(1+arg);
                 gd = csat*3*arg/vd;
 
-            } else {                            /* breakdown */
+            } else {
 
                 evrev = exp(-(here->DIOtBrkdwnV+vd)/vtebrk);
                 cd = -csat*evrev;
@@ -201,14 +202,28 @@ next1:
 
             }
 
-            gd = gd + ckt->CKTgmin;
-            cd = cd + ckt->CKTgmin*vd;
+            /*
+            if (ckt->CKTmode & MODEINITTRAN) {
+                *(ckt->CKTstate1 + here->DIOoldCurr) = cd;
+                *(ckt->CKTstate1 + here->DIOoldCond) = gd;
+                *(ckt->CKTstate2 + here->DIOoldCurr) = cd;
+                *(ckt->CKTstate2 + here->DIOoldCond) = gd;
+            }
+            else {
+                //cd = (*(ckt->CKTstate1 + here->DIOoldCurr) - (ckt->CKTdelta)*cd) / (1 + (ckt->CKTdelta)*vp);
+                //gd = (*(ckt->CKTstate1 + here->DIOoldCond) - (ckt->CKTdelta)*gd) / (1 + (ckt->CKTdelta)*vp);
+                //cd = (1 - vp)*cd + vp*(*(ckt->CKTstate1 + here->DIOoldCurr))*(ckt->CKTdelta);
+                //gd = (1 - vp)*gd + vp*(*(ckt->CKTstate1 + here->DIOoldCond))*(ckt->CKTdelta);
+            }
+            */
 
             if ((ckt->CKTmode & (MODEDCTRANCURVE | MODETRAN | MODEAC | MODEINITSMSIG)) ||
                      ((ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC))) {
               /*
                *   charge storage elements
                */
+
+                /* junction(depletion) charge */
                 czero=here->DIOtJctCap;
                 if (vd < here->DIOtDepCap){
                     arg=1-vd/here->DIOtJctPot;
@@ -222,56 +237,107 @@ next1:
                     deplcap = czof2*(here->DIOtF3+here->DIOtGradingCoeff*vd/here->DIOtJctPot);
                 }
 
-                if (model->DIOsoftRevRecParamGiven) {
+                /* diffusion charge */
+                if (ckt->CKTmode & MODEINITTRAN) {
+                    diffcharge = tt * cd;
+                    diffcap = tt * gd;
+                    *(ckt->CKTstate2 + here->DIOdiffCharge) = diffcharge;
+                    *(ckt->CKTstate2 + here->DIOdiffCap) = diffcap;
+                    *(ckt->CKTstate1 + here->DIOdiffCharge) = diffcharge;
+                    *(ckt->CKTstate1 + here->DIOdiffCap) = diffcap;
+                    *(ckt->CKTstate1 + here->DIOoldCurr) = cd;
+                    *(ckt->CKTstate1 + here->DIOoldCond) = gd;
+                }
+                else {
+                    double dt = ckt->CKTdelta;
+                    /*
+                    //Fixed-point iteration
+                    double diffcharge_old, diffcap_old;
 
-                    if (ckt->CKTmode & MODEINITTRAN) {
-                        diffcharge = tt * cd;
-                        diffcap = tt * gd;
-                        *(ckt->CKTstate2 + here->DIOdiffCharge) = diffcharge;
-                        *(ckt->CKTstate2 + here->DIOdiffCap) = diffcap;
-                        *(ckt->CKTstate1 + here->DIOdiffCharge) = diffcharge;
-                        *(ckt->CKTstate1 + here->DIOdiffCap) = diffcap;
-                        *(ckt->CKTstate1 + here->DIOoldCurr) = cd;
-                        *(ckt->CKTstate1 + here->DIOoldCond) = gd;
-                    }
-                    else {
-                        double dt = 1*ckt->CKTdelta;
+                    //Forward Euler predictor
+                    diffcharge = (*(ckt->CKTstate1 + here->DIOdiffCharge))*(1 - dt / (tt*vp)) + (dt*(*(ckt->CKTstate1 + here->DIOoldCurr)) / vp);
+                    diffcap = (*(ckt->CKTstate1 + here->DIOdiffCap))*(1 - dt / (tt*vp)) + (dt*(*(ckt->CKTstate1 + here->DIOoldCond)) / vp);
+                    do {
+                        diffcharge_old = diffcharge;
+                        diffcharge = (*(ckt->CKTstate1 + here->DIOdiffCharge)) + ((dt / vp)*(cd - (diffcharge_old / tt)));
+                    } while (fabs(diffcharge - diffcharge_old) > 1e-12);
+                    do {
+                        diffcap_old = diffcap;
+                        diffcap = (*(ckt->CKTstate1 + here->DIOdiffCap)) + ((dt / vp)*(gd - (diffcap_old / tt)));
+                    } while (fabs(diffcap - diffcap_old) > 1e-12);
+                    */
 
-                        //Backward Euler
-                        // Qk+1 = tt (VP Qk + h Ik+1) / (tt VP + h)
-                        diffcharge = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCharge)) + dt*cd) / (tt*vp + dt));
-                        diffcap = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCap)) + dt*gd) / (tt*vp + dt));
+                    /*
+                    //Newton's method
+                    double diffcharge_old, diffcap_old;
 
-                        //Trap
-                        // Qk+1 = (Qk (2 tt VP - h) + h tt (Ik + Ik+1)) / (2 tt VP + h)
-                        //diffcharge = ((*(ckt->CKTstate1 + here->DIOdiffCharge))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCurr)) + cd)) / (2.0*tt*vp + dt);
-                        //diffcap = ((*(ckt->CKTstate1 + here->DIOdiffCap))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCond)) + gd)) / (2.0*tt*vp + dt);
+                    //Forward Euler predictor
+                    diffcharge = (*(ckt->CKTstate1 + here->DIOdiffCharge))*(1 - dt / (tt*vp)) + (dt*(*(ckt->CKTstate1 + here->DIOoldCurr)) / vp);
+                    diffcap = (*(ckt->CKTstate1 + here->DIOdiffCap))*(1 - dt / (tt*vp)) + (dt*(*(ckt->CKTstate1 + here->DIOoldCond)) / vp);
+                    do {
+                        diffcharge_old = diffcharge;
+                        diffcharge = diffcharge_old + ((tt*vp*(*(ckt->CKTstate1 + here->DIOdiffCharge)) + tt*dt*cd - diffcharge_old*(dt + tt*vp)) / (dt + tt*vp));
+                    } while (fabs(diffcharge - diffcharge_old) > 1e-12);
 
-                        //Gear
-                        // Qk+1 = tt ((VP (4 Qk - Qk-1) + 2 h Ik+1) / (3 tt VP + 2 h))
-                        //diffcharge = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCharge)) - (*(ckt->CKTstate2 + here->DIOdiffCharge))) + 2.0*dt*cd) / (3.0*tt*vp + 2.0*dt));
-                        //diffcap = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCap)) - (*(ckt->CKTstate2 + here->DIOdiffCap))) + 2.0*dt*gd) / (3.0*tt*vp + 2.0*dt));
+                    do {
+                        diffcap_old = diffcap;
+                        diffcap = diffcap_old + ((tt*vp*(*(ckt->CKTstate1 + here->DIOdiffCap)) + tt * dt*gd - diffcap_old * (dt + tt*vp)) / (dt + tt*vp));
+                    } while (fabs(diffcap - diffcap_old) > 1e-12);
+                    */
 
-                    }
+                    //Backward Difference Operator
+                    //diffcharge = tt * (cd - vp * (((*(ckt->CKTstate1 + here->DIOdiffCharge)) - (*(ckt->CKTstate2 + here->DIOdiffCharge))) / (ckt->CKTdeltaOld[1])));
+                    //diffcap = tt * (gd - vp * (((*(ckt->CKTstate1 + here->DIOdiffCap)) - (*(ckt->CKTstate2 + here->DIOdiffCap))) / (ckt->CKTdeltaOld[1])));
 
-                    *(ckt->CKTstate0 + here->DIOdiffCap) = diffcap;
-                    *(ckt->CKTstate0 + here->DIOoldCurr) = cd;
-                    *(ckt->CKTstate0 + here->DIOoldCond) = gd;
+                    //Standard Diffy-Q w/ trap integral
+                    //diffcharge = (dt / (2.0*vt))*(exp(-dt / (tt*vt))*(*(ckt->CKTstate1 + here->DIOoldCurr)) + cd) + exp(-dt / (tt*vt))*(*(ckt->CKTstate1 + here->DIOdiffCharge));
+                    //diffcap = (dt / (2.0*vt))*(exp(-dt / (tt*vt))*(*(ckt->CKTstate1 + here->DIOoldCond)) + gd) + exp(-dt / (tt*vt))*(*(ckt->CKTstate1 + here->DIOdiffCap));
 
-                    //printf("time: %e vd: %e\n", ckt->CKTtime, vd);
-                    //printf("%e %e %e %e %e %e\n", ckt->CKTtime, ckt->CKTdelta, diffcharge, diffcap, cd, gd);
+                    //Forward Euler
+                    //diffcharge = ((*(ckt->CKTstate1 + here->DIOdiffCharge))*(tt*vp - dt) + dt*tt*cd) / (tt*vp);
+                    //diffcap = ((*(ckt->CKTstate1 + here->DIOdiffCap))*(tt*vp - dt) + dt*tt*gd) / (tt*vp);
 
-                } else {
+                    //Backward Euler
+                    diffcharge = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCharge)) + dt*cd) / (tt*vp + dt));
+                    diffcap = tt*((vp*(*(ckt->CKTstate1 + here->DIOdiffCap)) + dt*gd) / (tt*vp + dt));
 
-                    diffcharge = tt*cd;
-                    diffcap = tt*gd;
+                    //Trap
+                    //diffcharge = ((*(ckt->CKTstate1 + here->DIOdiffCharge))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCurr)) + cd)) / (2.0*tt*vp + dt);
+                    //diffcap = ((*(ckt->CKTstate1 + here->DIOdiffCap))*(2.0*tt*vp - dt) + dt*tt*((*(ckt->CKTstate1 + here->DIOoldCond)) + gd)) / (2.0*tt*vp + dt);
 
+                    //Gear
+                    //diffcharge = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCharge)) - (*(ckt->CKTstate2 + here->DIOdiffCharge))) + 2.0*dt*cd) / (3.0*tt*vp + 2.0*dt));
+                    //diffcap = tt*((vp*(4.0*(*(ckt->CKTstate1 + here->DIOdiffCap)) - (*(ckt->CKTstate2 + here->DIOdiffCap))) + 2.0*dt*gd) / (3.0*tt*vp + 2.0*dt));
+
+                    //Nothing
+                    //diffcharge = tt * cd;
+                    //diffcap = tt * gd;
                 }
 
-                *(ckt->CKTstate0 + here->DIOcapCharge) = deplcharge;
                 *(ckt->CKTstate0 + here->DIOdiffCharge) = diffcharge;
+                *(ckt->CKTstate0 + here->DIOdiffCap) = diffcap;
+                //*(ckt->CKTstate0 + here->DIOoldCurr) = cd;
+                //*(ckt->CKTstate0 + here->DIOoldCond) = gd;
 
-                here->DIOcap = deplcap + diffcap;
+                //printf("Time: %e\n", ckt->CKTtime);
+                //printf("Charge: %e\n", diffcharge);
+                //printf("Cap: %e\n", diffcap);
+                //printf("Curr: %e\n", cd);
+                //printf("Cond: %e\n", gd);
+                //printf("DelCharge: %e\n", diffcharge-(*(ckt->CKTstate1 + here->DIOdiffCharge)));
+                //printf("DelCap: %e\n", diffcap-(*(ckt->CKTstate1 + here->DIOdiffCap)));
+                //printf("DelCurr: %e\n", cd-(*(ckt->CKTstate1 + here->DIOoldCurr)));
+                //printf("DelCond: %e\n\n", gd-(*(ckt->CKTstate1 + here->DIOoldCond)));
+
+                //diffcharge = here->DIOtTransitTime*cd;
+                *(ckt->CKTstate0 + here->DIOcapCharge) =
+                        diffcharge + deplcharge;
+
+                //diffcap = here->DIOtTransitTime*gd;
+
+                capd = diffcap + deplcap + here->DIOcmetal + here->DIOcpoly;
+
+                here->DIOcap = capd;
 
                 /*
                  *   store small-signal parameters
@@ -279,7 +345,7 @@ next1:
                 if( (!(ckt->CKTmode & MODETRANOP)) ||
                         (!(ckt->CKTmode & MODEUIC)) ) {
                     if (ckt->CKTmode & MODEINITSMSIG){
-                        *(ckt->CKTstate0 + here->DIOcapCurrent) = deplcap + diffcap;
+                        *(ckt->CKTstate0 + here->DIOcapCurrent) = capd;
 
                         if(SenCond){
                             *(ckt->CKTstate0 + here->DIOcurrent) = cd;
@@ -287,7 +353,7 @@ next1:
 #ifdef SENSDEBUG
                             printf("storing small signal parameters\n");
                             printf("cd = %.7e,vd = %.7e\n",cd,vd);
-                            printf("capd = %.7e ,gd = %.7e \n",deplcap+diffcap,gd);
+                            printf("capd = %.7e ,gd = %.7e \n",capd,gd);
 #endif /* SENSDEBUG */
                         }
                         continue;
@@ -302,7 +368,7 @@ next1:
                         printf("storing parameters for transient sensitivity\n"
                                 );
                         printf("qd = %.7e, capd = %.7e,cd = %.7e\n",
-                                *(ckt->CKTstate0 + here->DIOcapCharge),deplcap+diffcap,cd);
+                                *(ckt->CKTstate0 + here->DIOcapCharge),capd,cd);
 #endif /* SENSDEBUG */
                         continue;
                     }
@@ -310,22 +376,14 @@ next1:
                     if (ckt->CKTmode & MODEINITTRAN) {
                         *(ckt->CKTstate1 + here->DIOcapCharge) =
                                 *(ckt->CKTstate0 + here->DIOcapCharge);
-                        *(ckt->CKTstate1 + here->DIOdiffCharge) =
-                                *(ckt->CKTstate0 + here->DIOdiffCharge);
                     }
-                    error = NIintegrate(ckt,&geq,&ceq,deplcap,here->DIOcapCharge);
+                    error = NIintegrate(ckt,&geq,&ceq,capd,here->DIOcapCharge);
                     if(error) return(error);
                     gd=gd+geq;
                     cd=cd+*(ckt->CKTstate0 + here->DIOcapCurrent);
-                    error = NIintegrate(ckt,&geq,&ceq,diffcap,here->DIOdiffCharge);
-                    if(error) return(error);
-                    gd=gd+geq;
-                    cd=cd+*(ckt->CKTstate0 + here->DIOdiffCurrent);
                     if (ckt->CKTmode & MODEINITTRAN) {
                         *(ckt->CKTstate1 + here->DIOcapCurrent) =
                                 *(ckt->CKTstate0 + here->DIOcapCurrent);
-                        *(ckt->CKTstate1 + here->DIOdiffCurrent) =
-                                *(ckt->CKTstate0 + here->DIOdiffCurrent);
                     }
                 }
             }
@@ -336,7 +394,7 @@ next1:
              *   check convergence
              */
             if ( (!(ckt->CKTmode & MODEINITFIX)) || (!(here->DIOoff))  ) {
-                if (Check == 1)  {
+                if (Check == 1) {
                     ckt->CKTnoncon++;
                     ckt->CKTtroubleElt = (GENinstance *) here;
                 }