/* =============================== PNP1 library ===================================These functions let you set up a dynamic workspace adapted to your currentPNP problem (involving a number of ion species and a geometrical grid) andsolve the PNP equations in a given geometry and using given boundaryconditions and PNP parameters. The layout of a pore geometry and thetranslation of far-bulk into domain-boundary conditions are not topics ofthe PNP solver (they are facilitated by other library modules).PNP theory is solved here in one dimension, assuming homogeneity in the others.The metrics of the domain are expressed by an axial coordinate and thearea of the surface of homogeneity that belongs to the axial coordinate.Exact solutions are obtained for both cylindrical and conical domains (with flat or spherically curved surfaces, respectively). Approximate solutions areobtained for domains that are piecewise cylindrical and conical (an exampleis a domain that includes a cylindrical pore proper, tapered atria, andhemispherical subdomains of the bulks). The PNP solver is given the domainmetrics as an array of axial grid coordinates and the corresponding domainsurface areas; it is otherwise unaware of the shape of the domain, and thuscapable of dealing with arbitrary geometries. Furthermore, it works withthe grid defined by the given axial coordinates. This grid generally will benon-uniform (with intervals made proportional to domain surface area).The implemented theory makes provisions to account for two forms of specificshort-range interaction that ions may experience in a PNP domain. (1) binding tosaturable sites (which immobilizes the ion), and (2) non-saturable partitioning(with potential effects on ion mobility as expressed in the diffusion coefficient).For instance, form (1) may apply to protons that protonize structural carboxylgroups, whereas form (2) may apply to alkali cations that interact with an arrayof structural carbonyl groups. Form (2) is suitable also for introducing Bornenergy variations associated with a non-uniform dielectric permittivity (see below).Several PNP parameters are allowed to vary along the axial coordinate and need tobe specified as axial profiles (on the grid of the axial coordinate array):  - area of domain surface (see above)  - relative dielectric permittivity  - diffusion coefficients of ions  - 'partition energies' of ions (standard chemical potentials)  - structural fixed-charge concentrationSaturable binding also is allowed to vary along the pore axis. This is specifiedin a client-provided function that computes 'total' ion concentrations alongthe pore from the given 'free' concentrations (a simple equilibrium computation,since the 'free' concentrations in a stationary flux are constant in time).The PNP solver treats the bulk boundaries as providing flux-independent ionconcentrations and electric potentials. For PNP domains that extend well intothe bulks, the solver can use the far-bulk values. Other boundary treatments(like Donnan or Guy-Chapman theory) are possible; then the PNP solver is giventhe domain-side boundary values produced by the respective boundary treatment.The library comprises the front-end functions:                         - build_PNP_ws: dynamically allocates a workspace memory for the variables that                 transport information between the PNP solver and its clients.                 Several PNP workspaces can coexist and can be used randomly;                 the client is responsible for de-allocation.                   - init_PNP_var: primes internal arrays for the PNP solver that need to be                 computed only once for many uses of the solver (but a change of                 geometry invalidates these internal arrays).                  - solve_PNP: solves the PNP equations for one parameter set  CAVE: This machinery defines its variables dynamically. You cannot expect your       C compiler to protect you from indexing over bounds when you set up PNP       input variables. Furthermore, the functions do not test your inputs: they       execute them verbatim.     */#include "PNP0.C"#include "math.h"#include "stdio.h"/*================================ Semaphore Builder =============================*//* ----------------------------------- front end -----------------------------------   uses a client-supplied function to reserve workspace memory for PNP and sets up   the PNP arrays and the semaphore to access them; info about the allocated memory   (for future de-allocation) is included in the semaphore ('memory', 'nbytes').      The allocator is expected to return a null pointer on failure. Then build_PNP_ws   also returns a null pointer. */PNP_SMP *build_PNP_ws(L nions, L ngrid, void *allocator()){L nb, temp, kion;  B *mem, *free; PNP_SMP *smp;nb = sizeof(PNP_SMP) + ((nions * 5 + 13) * ngrid     + (ngrid-2) * 9  + nions * 4) * sizeof(DF) + nions * 5 * sizeof(DF *);if (!(mem = (B *)allocator(nb))) return(NULL);smp = (PNP_SMP *)mem; free = mem + sizeof(PNP_SMP);smp->memory = mem;smp->nbytes = nb;smp->Nions = nions;smp->Ngrid = ngrid;smp->Ds     = (DF **)free; free += (temp = nions * sizeof(DF *));smp->E0s    = (DF **)free; free += temp;smp->Es     = (DF **)free; free += temp;smp->Cs     = (DF **)free; free += temp;smp->TCs    = (DF **)free; free += temp;smp->X      = (DF *)free; free += (temp = ngrid * sizeof(DF));smp->Area   = (DF *)free; free += temp;smp->Radius = (DF *)free; free += temp;smp->F      = (DF *)free; free += temp;smp->EPS    = (DF *)free; free += temp;smp->V      = (DF *)free; free += temp;smp->WA     = (DF *)free; free += temp;smp->WB     = (DF *)free; free += temp;smp->WC     = (DF *)free; free += temp;smp->WR     = (DF *)free; free += temp;smp->WU     = (DF *)free; free += temp;smp->WT     = (DF *)free; free += temp;smp->H      = (DF *)free; free += temp;for (kion=0; kion<nions; kion++)    { (smp->Ds)[kion]   = (DF *)free; free += temp;      (smp->E0s)[kion]  = (DF *)free; free += temp;      (smp->Es)[kion]   = (DF *)free; free += temp;      (smp->Cs)[kion]   = (DF *)free; free += temp;      (smp->TCs)[kion]  = (DF *)free; free += temp;    }    smp->Z     = (DF *)free; free += (temp = nions * sizeof(DF));smp->CL    = (DF *)free; free += temp;smp->CR    = (DF *)free; free += temp;smp->FLUX  = (DF *)free; free += temp;smp->al1   = (DF *)free; free += (temp = (ngrid-2) * sizeof(DF));smp->al2   = (DF *)free; free += temp;smp->al3   = (DF *)free; free += temp;smp->bt1   = (DF *)free; free += temp;smp->bt2   = (DF *)free; free += temp;smp->bt3   = (DF *)free; free += temp;smp->owa   = (DF *)free; free += temp;smp->owb   = (DF *)free; free += temp;smp->owc   = (DF *)free; free += temp;return(smp);}/*======================= Initializer of PNP variables ========================   Initializes difference coefficients, voltage profile, and some PNP control   variables. Use after building semaphore and before solving PNP.      Note: this defaults the solver to the 'verbose' mode, in which it prints   the maximal voltage change for each iteration as an index of convergence.   The usefulness of this info hardly can be overstated.*/void init_PNP_var(PNP_SMP *smp){L k; DF temp, *hp, *h1p, *al1p, *al2p, *al3p, *bt1p, *bt2p, *bt3p;for (k=1; k<smp->Ngrid; k++) (smp->H)[k] = (smp->X)[k] - (smp->X)[k-1];for (hp=smp->H+1, h1p=hp+1,     al1p=smp->al1, al2p=smp->al2, al3p=smp->al3,     bt1p=smp->bt1, bt2p=smp->bt2, bt3p=smp->bt3,     k=1; k<smp->Ngrid-1;     k++, hp++, h1p++    )    { *al1p = 0.5 / *h1p; *al3p = -0.5 / *hp; *(al2p++) = - ( *(al1p++) + *(al3p++) );      temp = 2.0 / ( *h1p + *hp );      *bt1p = temp / *h1p; *bt3p = temp / *hp; *(bt2p++) = - (*(bt1p++) + *(bt3p++) );    }for (k=0; k<smp->Ngrid; k++) (smp->V)[k] = 0.0;smp->MaxIter = 200;smp->Tolerance = 1e-8;smp->verbose = TRUE;}/*=============================== PNP solver =====================================*//*----------------------------- solver library -----------------------------------*//*--- find absolute maximum of an array */static DF maximum(DF max, DF *arr, L n){ L k; DF temp;for (k=0;k<n;k++) { if ((temp = fabs(*(arr++))) > max) max = temp; }return(max);}/*--- integrate array in place and return value of integral (results are *2) */static DF integrate(DF *h, DF *y, L n){DF temp, sum; L k;temp = *(y); *(y++) = sum = 0.0;for (k=1; k<n; k++) {  sum += (temp + *y) * *(++h); temp = *y; *(y++) = sum; }return(sum);}/*--- compute energy and concentration profiles, and raw fluxes */static void comp_EC_prof(PNP_SMP *smp){L kion, k;DF *Ep, *Vp, *E0p, *WRp, *Dp, *Areap, *Cp, *cp, *wrp, *ep, z, temp, temp1, intL;for (kion=0; kion<smp->Nions; kion++)    { z = (smp->Z)[kion];      for (Ep = (smp->Es)[kion], Vp = smp->V, E0p = (smp->E0s)[kion],           WRp = smp->WR, Dp=(smp->Ds)[kion], Areap=smp->Area,           temp = smp->kT_e, k=0; k<smp->Ngrid; k++          )          { *Ep = *(Vp++) * z + *(E0p++) / temp;            *(WRp++) = exp(*(Ep++)) / (*(Dp++) * *(Areap++));          }      intL = integrate(smp->H, smp->WR, smp->Ngrid);      Cp = (smp->Cs)[kion]; Ep = (smp->Es)[kion];      temp = exp(Ep[0]) * Cp[0];      temp1 = exp(Ep[smp->Ngrid-1]) * Cp[smp->Ngrid-1];      for (cp=Cp+1, wrp=smp->WR+1, ep=Ep+1,           k=1; k<smp->Ngrid-1; k++, ep++, wrp++          )          { *(cp++) = ( (intL - *wrp) * temp + *wrp * temp1                      ) / ( exp(*ep) * intL );          }      (smp->FLUX)[kion] = ( temp - temp1) * 2.0 / intL;    }(*smp->comp_TC)(smp);}/*--- build discretisized Poisson equations */static void build_Poisson(PNP_SMP *smp){L kion, k;DF *wap, *wbp, *wcp, *owap, *owbp, *owcp, *wrp, *Fp, temp, temp2,   *Vp, z, z2, *vp, *tcp;Vp = smp->V; temp2 = 1.0 / smp->Cscale;for ( owap=smp->owa, owbp=smp->owb, owcp=smp->owc,      wap=smp->WA+1, wbp=smp->WB+1, wcp=smp->WC+1,      wrp=smp->WR+1, Fp=smp->F+1,      k=1; k<smp->Ngrid-1; k++    )    { *(wap++) = *(owap++); *(wbp++) = *(owbp++); *(wcp++) = *(owcp++);      *(wrp++) = - *(Fp++) * temp2;    }for (kion=0; kion<smp->Nions; kion++)    { z = (smp->Z)[kion]; z2 = z * z;       for (wbp=smp->WB+1, wrp=smp->WR+1, vp=smp->V+1, tcp = (smp->TCs)[kion]+1,           k=1; k<smp->Ngrid-1; k++          )          { temp = *tcp * z2;            *(wbp++) -= temp; *(wrp++) -= temp * *(vp++) + *(tcp++) * z;          }    }smp->WR[smp->Ngrid-2] -= smp->V[smp->Ngrid-1] * smp->WC[smp->Ngrid-2];smp->WR[1] -= smp->V[0] * smp->WA[1];}/*--- solve tridiagonal linear eqn system   - wa holds the subdiagonal matrix elements (the first element is ignored)   - wb holds the diagonal elements   - wc holds the supradiagonal elements (the last element is ignored)   - wr holds the right-hand sides   - wu receives the solution   - wg is an array for temporaries   - returned bool indicates whether system was solvable*/static BOOL solvetridiag(PNP_SMP *smp){DF *wap, *wbp, *wcp, *wrp, *wup, *wtp, bet; L j, n;wap = smp->WA+1; wbp = smp->WB+1; wcp = smp->WC+1;wrp = smp->WR+1; wup = smp->WU+1; wtp = smp->WT+1;n = smp->Ngrid - 2;if (wbp[0] == 0.0) return(FALSE);wup[0] = wrp[0] / (bet = wbp[0]);for (j=1; j<n; j++, wap++, wbp++, wcp++, wrp++, wup++, wtp++)	{ wtp[1] = wcp[0] / bet;	  bet = wbp[1] - wap[1] * wtp[1];	  if (bet == 0.0) return(FALSE);	  wup[1] = (wrp[1] - wap[1] * wup[0]) / bet;	}for (j = (n-2), wup=smp->WU+1+j, wtp=smp->WT+1+j;     j >=0;     j--, wup--, wtp--    ) wup[0] -= wtp[1] * wup[1];return(TRUE);}/*--- compute error and new V */static DF new_V_and_err(PNP_SMP *smp){L k;DF *wap, *wup, *vp;for ( wap=smp->WA+1, wup=smp->WU+1, vp=smp->V+1,      k=1; k<smp->Ngrid-1; k++    )    { *(wap++) = *wup - *vp; *(vp++) = *(wup++); }return(maximum(0.0,smp->WA+1,smp->Ngrid-2));}/* ------------------------------ front end ---------------------------------------- - imports/exports data via semaphore structure (see PNP0.C) - the returned boolean signals convergence within the specified voltage tolerance - expects pore geometry to be defined in X, Area (you can use the geometry builder   for setting these up) - expects boundary concentrations to be presented in CL, CR (it is your    responsibility to compute these from bulk concentrations if, e.g., surface charge   effects are to be considered) - uses V profile from preceding instance as first approximation (adding a linear   update for a change in VM); therefore, V should be nulled prior to the first   instance (as is done automatically by the initializer function)   */BOOL solve_PNP(PNP_SMP *smp){ L k, kiter;DF PoissonF, Fluxscale, temp, temp1;DF *al1p, *al2p, *al3p, *bt1p, *bt2p, *bt3p, *areap, *epsp,   *owap, *owbp, *owcp;DF *Vp, *C;/*-- design scaling factors */temp = maximum(0.0,smp->CL,smp->Nions);smp->Cscale = maximum(temp,smp->CR,smp->Nions);smp->kT = Boltzmann * (Celsius0 + smp->Celsius);smp->kT_e = smp->kT / e0;PoissonF = epsilon0 * smp->kT_e * Liter /          ( e0 * Avogadro * smp->Cscale * Angstrom * Angstrom );Fluxscale = cm2 * Avogadro * smp->Cscale * Angstrom * Angstrom / (Liter * Angstrom);/*-- set up originals of left-hand arrays for Poisson */for ( al1p=smp->al1, al2p=smp->al2, al3p=smp->al3,      bt1p=smp->bt1, bt2p=smp->bt2, bt3p=smp->bt3,      epsp=smp->EPS+1, areap=smp->Area+1,       owap=smp->owa, owbp=smp->owb, owcp=smp->owc,      k=1; k<smp->Ngrid-1;      k++, al1p++,al2p++,al3p++,bt1p++,bt2p++,bt3p++,epsp++,areap++,      owap++,owbp++,owcp++    )    {      temp = ( al1p[0] * areap[1] + al2p[0] * areap[0] + al3p[0] * areap[-1] )             * epsp[0] / areap[0]              + al1p[0] * epsp[1] + al2p[0] * epsp[0] + al3p[0] * epsp[-1];      owcp[0] = ( epsp[0] * bt1p[0] + temp * al1p[0] ) * PoissonF;      owbp[0] = ( epsp[0] * bt2p[0] + temp * al2p[0] ) * PoissonF;        owap[0] = ( epsp[0] * bt3p[0] + temp * al3p[0] ) * PoissonF;    }/* printf("\nbt = %e, %e, %e",(smp->bt1)[500],(smp->bt2)[500],(smp->bt3[500])); *//*-- insert scaled boundary concentrations; update the V array of the previous     instance for a change in VM (adding a linear correction to the profile)*/for (k=0; k<smp->Nions; k++)    { C = (smp->Cs)[k]; C[0] = smp->CL[k] / smp->Cscale;      C[smp->Ngrid-1] = smp->CR[k] / smp->Cscale;    }if((temp = smp->VM / smp->kT_e - (smp->V)[0]) != 0.0)    { for (k=0, Vp = smp->V, temp1 = smp->Ngrid - 1; k<smp->Ngrid; k++)          { *(Vp++) += temp * (temp1 - k) / temp1;          }    }     /*-- solve PNP by iteration *//* printf("\n%e",(smp->V)[smp->Ngrid/2]); */kiter = smp->MaxIter;while ((kiter--) > 0)	{ 	  comp_EC_prof(smp);/*	  printf("\n%e %e",((smp->TCs)[0])[smp->Ngrid/2],((smp->TCs)[1])[smp->Ngrid/2]);*/ 	  build_Poisson(smp);/*	  printf("\n%e %e %e %e",smp->WA[500],smp->WB[500],smp->WC[500],smp->WR[500]); */	  if (!solvetridiag(smp)) return(FALSE); 	  temp = new_V_and_err(smp); 	  if (smp->verbose) printf("\ndVmax = %e",temp);	  if (temp < smp->Tolerance) break;	}/*-- compute final concentration profiles, calibrated fluxes, and net current */comp_EC_prof(smp);smp->NetCurrent = 0.0;for (k=0; k<smp->Nions; k++)    { temp = ((smp->FLUX)[k] *= Fluxscale);      smp->NetCurrent += temp * (smp->Z)[k] * e0 / pA;    }return(TRUE);}