Variational Statement of the Poisson Eq.

Consider the following functional of a scalar field $\phi (\vec r)$:

$$S[\phi] = \int d\vec r [\frac{\epsilon}{8\pi} \vert\vec \nabla \phi (\vec r) \vert^2 - \rho (\vec r) \phi (\vec r) ]$$ (7)

where the integral is over a prescribed region of space. In this expression $\epsilon$ is the dielectric constant (suppressing, for simplicity, any spatial dependence) and $\rho (\vec r)$ is a fixed charged distribution. We shall term this functional an ``Action", based on rough analogies with variational principles of classical mechanics [7] and more direct analogies with field-theoretical derivations of the PB Eq. (cf. Sect. 7 below).

If we consider all possible functions $\phi (\vec r)$satisfying the boundary condition that $\phi$ has prescribed values on the surface, then the $\phi$ field which yields the minimum value for the r.h.s. of Eq. 7 is the one that solves Poisson's Eq. (6) ! This can be shown via a standard application of the calculus of variations. That is, substitute $\phi = \bar \phi + \delta \phi$ (such that $\delta \phi =0$ on the boundary surface) and determine the condition on the function $\bar \phi$ under which the change in S for any small variation $\delta \phi$ around $\bar \phi$ is zero. The condition is precisely that $\bar \phi$ satisfies the Poisson Eq. [7].

Another form of S may be obtained under the condition (relevant for many physical problems) that either $\phi$ or $\vec \nabla \phi$ is zero on the boundary surface. Then we can integrate the first term on the r.h.s. of Eq. 7 by parts and neglect the surface term, so that:

$$S[\phi] = \int d\vec r [-\frac{\epsilon}{8\pi} \phi (\vec r) \nabla^2 \phi (\vec r) - \rho (\vec r) \phi (\vec r) ]$$ (8)

Computing the functional derivative of S with respect to $\phi (\vec r)$ [11] and setting it equal to zero again yields the Poisson equation. These variational statements can be put to practical use. By considering a discretized (lattice) form of Eq. (7) or (8), a simple and effective annealing or relaxation strategy for computing the $\phi$ field is suggested.

In particular, consider the lattice-discretized form of Eq. 8, namely:

$$S = -\frac{\alpha}{2} \sum_{\vec n \vec m} \phi_{\vec n} \De... ... m} \phi_{\vec m} - \sum_{\vec n} \rho_{\vec n} \phi_{\vec n}$$ (9)

where $\alpha = \frac {h \epsilon}{4 \pi}$ (a positive constant with h the grid spacing), $\rho_{\vec n}$is the total charge on lattice site $\vec n$ (the charge density multiplied by h³ ; $\vec n$is a triple of integers that specifies position on the lattice), and $\Delta_{\vec n \vec m}$is h² times the discretized Laplacian operator or matrix. Since this matrix plays an important role in both the theory and practical implementation of Poisson and Poisson-Boltzmann methods, we pause to note its basic properties.

The discretized Laplacian operator is a finite difference (FD) approximation to the second derivative operator, summed over three Cartesian directions. In one dimension, $h^2 \frac{d^2 \phi}{dx^2}_{FD} = (\phi_{i-1} - 2 \phi_i + \phi_{i+1} )$. [Higher order approximations are possible. These are considered in Section 2.4.4.] Thus the second derivative function can be thought of as a vector (indexed by position along the x axis) which can be obtained by linear transformation of the function $\phi$ (also represented by a vector on a 1-d lattice): $h^2 (d^2 \phi / dx^2 )_i = \sum_j \partial_{i,j}^2 \phi_j$.Here $\partial_{i,j}^2$ is the 1-d finite difference second derivative matrix, with non-zero elements of -2 on the diagonal and +1 one off-diagonal. Written out in full, the matrix looks like:

$$\left[ \begin{array} {rrrrrr} -2 & 1 & 0 & . & . & 0\ 1 &-... ... & . & . & . & . & .\ 0 & . & . & . & . & .\end{array}\right]$$ (10)

The element in lower left and upper right hand corners is appropriate for the boundary condition that the field has the value at the edge of the grid. Other boundary conditions are possible, for example, periodic or ``wrap-around" conditions. In this case, we replace with 1 in the lower left and upper right hand corners, but otherwise the matrix is unchanged.

Irrespective of boundary conditions, it is easy to see that this matrix is sparse and highly banded. It connects each (field) point on the lattice to its nearest neighbors, and only to its nearest neighbors.

The 3-D Laplacian matrix is simply the direct product of $\partial_{i,j}^2$ for, say, the x direction and unit matrices in y and z directions, plus analogous terms where $\partial_{i,j}^2$ corresponds to y and z directions. If we label a 3-D Lattice point by $\vec n$, then the non-zero elements of $\Delta_{\vec n \vec m}$ are the diagonal element $\Delta_{\vec n \vec n}= -6$and elements (value +1) that connect $\vec n$ to any of its six nearest neighbors on the lattice $\Delta_{\vec n,\vec n \pm \hat e}$, where $\hat e$ is a unit vector in the x,y, or z directions. Thus each row/column has seven non-zero elements. Just as in the 1-d case, the matrix is highly banded and encompasses self-interaction plus nearest neighbor coupling on the lattice [13].