This module helps to solve the 1-D and radial Poisson-Boltzmann equation.

Examples

Simple solve

Solves a PB equation for a planar geometry, a surface density of 1, using a fixed derivative as a left boundary condition. Then outputs the resulting potential profile to a file.

import polypbrenpkg / pbsolv

let
  ds=0.1
  zv = [1.0, -1.0]
  conc = [0.5, 0.5]
var eq = initPBEquation(smin = 0.0, smax = 10.0, sbin = 0.1,
                        kind=PlanePBE, zv = zv, concv=conc,
                        bvL=BCder, bcR=BCDer, left = 1.0, right = 0.0,
                        charge = 1.0, lambdab = 1.0)
# initialize the potential profile to 1.
eq.initialize(1.0)
# solves the equation.
let number_iterations = eq.solve(niter = 5000, tol = 0.00001)
var fout = open("phi.dat",fmWrite)
for i,x in eq.mesh:
  fout.writeLine(x, ' ',eq.phi[i])
fout.close

With a condition

Solves a PB equation for radial geometry with a condition that should converge to 0 after enough iterations. Here it is a simplified pH condition.

import polypbrenpkg / pbsolv

let
  pKa = 7.7
  siteDens = 5.5  # density of sites
  dqSite = -1.0    # charge difference of a site
  pHref = 9.0     # imposed pH

proc pHcond(e: PBEquation): float =
  let
    alpha = abs(sigma(e))/siteDens   # sigma() computes the charge density (defined in this module)
    pH = pKa + (ln(alpha/(1.0-alpha)) + dqSite*e.phi[0])/ln(10.0)
  result = pH - pHref

var eq = initPBEquation(smin = 8.0, smax = 30.0, sbin = 0.1,
                        kind=RadialPBE, zv = [1.0,-1.0], concv=[6.022e-3,6.022e-3],
                        lambdab = 0.7105)
# initialize the potential profile to -1.
eq.initialize(-1.0)
# solves the equation.
let number_iterations = eq.solveByCondition(condition = pHcond, minVal = -10.0, maxVal = 0.0,
                                              nitermax = 50, niterin = 5000, tolIn = 1.0e-4)

Types

EqModelKind = enum PlanePBE, RadialPBE
BoundaryCondition = enum BCval, BCder
PBEquation = object pars: PBEquationPars grid: int h, h2: float phi*: seq[float] phi2: seq[float]

Procs

proc defMesh(pbe: var PBEquation; smin, smax, sbin: float; reDefSbin = true; gridStep = 2) {...}{. raises: [ValueError], tags: [].}

Defines (or redefines) the space grid for equation pbe. smin is the lower bound, smax the higher, sbin the space interval. If reDefSbin is set to true, then sbin is redefined to match the [sbin, smax] interval and gridStep. Dimensions of potential tables : N = int((smin-smax)/pas) raises ValueError if the grid is too small (number of points < 3)

proc charge=(e: var PBEquation; charge: float) {...}{.raises: [], tags: [].}

sets the surface charge density of the left (inner) surface

proc rho=(e: var PBEquation; rho: float) {...}{.raises: [], tags: [].}

sets the density of an uniform background charge

proc lambdab=(e: var PBEquation; lambdab: float) {...}{.raises: [], tags: [].}

changes the bjerrum length

proc defBoundCond(e: var PBEquation; bcL, bcR: BoundaryCondition; left = 0.0; right = 0.0) {...}{. raises: [], tags: [].}

Defines the boundary conditions : left (or inner radius) (type is bcL, value is left), and right (or outer radius) (type is bcL, value is left) conditions.

proc initPBEquation(smin, smax, sbin: float; zv, concv: openArray[float]; kind = PlanePBE; bcL = BCval; bcR = BCder; left = 0.0; right = 0.0; lambdab = 0.715; charge = 0.0; rho = 0.0; gridStep = 2): PBEquation {...}{. raises: [ValueError], tags: [].}

Returns a new PBEquation object.

Parameters:

smin: minimum of radius/x coordinate
smax: maximum of radius/x coordinate
sbin: requested separation between points
zv: charges of ion species
conv: bulk concentration of ion species
kind: defines a planar geometry or a spherical geometry
bcL: type of boundary condition on the "left" (inner radius)
bcR: type of boundary condition on the "right" (inner radius)
left: value of the left boundary condition
right: value of the right boundary condition
charge: surface charge density of the left (inner) surface
rho: charge density of the background charge

proc initialize(e: var PBEquation; f: proc (x: float): float {...}{.closure.}) {...}{.raises: [], tags: [].}

Initialize the potential with a function.

proc initialize(e: var PBEquation; x: float) {...}{.raises: [], tags: [].}

Initialize the potential with an uniform value.

proc mesh(e: var PBEquation): seq[float] {...}{.raises: [], tags: [].}

gives the grid used by e as a sequence of x coordinates (or radial coordinates).

proc sigma(e: PBEquation): float {...}{.inline, raises: [], tags: [].}

Computes the inner (ion+background) surface charge density of the system

proc sigmaCum(e: PBEquation): seq[float] {...}{.noSideEffect, raises: [], tags: [].}

Computes the cumulative (ion+background) surface charge density of the system

proc sigmaCumSansBack(e: PBEquation): seq[float] {...}{.noSideEffect, raises: [], tags: [].}

Computes the cumulative (ions only) surface charge density of the system

proc residual(e: PBEquation): float {...}{.raises: [], tags: [].}

proc solve(eq: var PBEquation; niter = 10000; tol = 1e-05): int {...}{.noSideEffect, raises: [], tags: [].}

Solves the Equation for the boundary conditions previously set, in a maximum of niter iterations, for a tolerance of tol on the residue. Returns the number of iterations effectively computed.

proc solveByNeutralization(eq: var PBEquation; tolF = 0.001; tolIn = 1e-05; nitermax = 500; niterin = 10000; minVal = 0.0; maxVal = 10.0): int {...}{. raises: [IOError, ValueError], tags: [WriteIOEffect].}

Solves the Equation by matching the input eq.charge (placed on left or inner surface) with the calculated ionic + outer static charge iteratively, with a maximum of nitermax iterations, with niterin steps for each iteration. Returns the number of iterations effectively computed.

proc solveByCondition(eq: var PBEquation; condition: proc (e: PBEquation): float {...}{.closure.}; tolF = 0.001; tolIn = 1e-05; nitermax = 500; niterin = 10000; minVal = 0.0; maxVal = 10.0): int {...}{.raises: [IOError, ValueError], tags: [WriteIOEffect].}

Solves the Equation by nullifying iteratively the conditional function condition with a maximum of nitermax iterations, with niterin steps for each iteration. Returns the number of iterations actually computed.

proc residuals(e: PBEquation): seq[float] {...}{.raises: [], tags: [].}