%% LyX 1.5.6 created this file.  For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[english,preprint]{revtex4}
\usepackage[T1]{fontenc}
\usepackage[latin9]{inputenc}
\usepackage{graphicx}
\usepackage{babel}

\begin{document}

\title{MATLAB version of the APBS}


\date{\today}

\begin{abstract}
This is the first version. This solver uses the biconjugate gradient
stabilized method and the inexact LU decomposition to numerically
solve the linearized PB equation on the finest (target) 3D-grid. This
version requires the shifted dielectric and the ion accessibility
coefficient (kappa function) maps as generated by the APBS code as
well as the corresponding pqr file generated by the pdb2pqr code.
It uses standard three-linear splines (spl0) to spread the charge
density along the nearest grid points if needed. The resulting electrostatic
potential and charge maps are saved in dx format. For visualization
purpose, this code also generates two files (.fig and .tiff) corresponding
to the graphical representation of the electrostatic potential surface.
\end{abstract}
\maketitle

\subsection*{Description}

This code is based on Michel Holst's thesis and Nathan Baker's APBS
approach. The box-method is used to discretize the following (linearized)
PB equation 

\begin{equation}
-\nabla.\left(\epsilon\left(\mathbf{r}\right)\nabla u\left(\mathbf{r}\right)\right)+\bar{\kappa}\left(\mathbf{r}\right)u\left(\mathbf{r}\right)=magic\sum_{i=1}^{N}z_{i}\delta\left(\mathbf{r}-\mathbf{r}_{i}\right)\label{eq:one}\end{equation}


where $u\left(\mathbf{r}\right)=e_{c}\Phi\left(\mathbf{r}\right)/K_{B}T$
and $magic=4\pi e_{c}^{2}/K_{B}T.$ For a diagonal dielectric tensor,
the resulting discretized linear PB equations at the nodes $u_{ijk}=u\left(x_{i},y_{j},z_{k}\right)$
for $1\leq i\leq N_{x}$, $1\leq j\leq N_{y}$ and $1\leq k\leq N_{z}$
reads 

\[
\left[\epsilon_{i-1/2,j,k}^{x}\frac{\left(h_{j-1}+h_{j}\right)\left(h_{k-1}+h_{k}\right)}{4h_{i-1}}+\epsilon_{i+1/2,j,k}^{x}\frac{\left(h_{j-1}+h_{j}\right)\left(h_{k-1}+h_{k}\right)}{4h_{i}}+\right.\]


\[
\epsilon_{i,j-1/2,k}^{y}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{k-1}+h_{k}\right)}{4h_{j-1}}+\epsilon_{i,j+1/2,k}^{y}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{k-1}+h_{k}\right)}{4h_{j}}+\]


\[
\epsilon_{i,j,k-1/2}^{k}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{j-1}+h_{j}\right)}{4h_{k-1}}+\epsilon_{i,j,k+1/2}^{k}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{j-1}+h_{j}\right)}{4h_{k}}+\]


\[
\left.\kappa_{ijk}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{j-1}+h_{j}\right)\left(h_{k-1}+h_{k}\right)}{8}\right]u_{ijk}+\]


\[
\left[-\epsilon_{i-1/2,j,k}^{x}\frac{\left(h_{j-1}+h_{j}\right)\left(h_{k-1}+h_{k}\right)}{4h_{i-1}}\right]u_{i-1jk}+\left[-\epsilon_{i+1/2,j,k}^{x}\frac{\left(h_{j-1}+h_{j}\right)\left(h_{k-1}+h_{k}\right)}{4h_{i}}\right]u_{i+1jk}+\]


\[
\left[-\epsilon_{i,j-1/2,k}^{y}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{k-1}+h_{k}\right)}{4h_{j-1}}\right]u_{ij-1k}+\left[-\epsilon_{i,j+1/2,k}^{y}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{k-1}+h_{k}\right)}{4h_{j}}\right]u_{ij+1k}+\]


\[
\left[-\epsilon_{i,j,k-1/2}^{k}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{j-1}+h_{j}\right)}{4h_{k-1}}\right]u_{ijk-1}+\left[-\epsilon_{i,j,k+1/2}^{k}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{j-1}+h_{j}\right)}{4h_{k}}\right]u_{ijk+1}=\]


\begin{equation}
magic\frac{\left(h_{i-1}+h_{i}\right)\left(h_{j-1}+h_{j}\right)\left(h_{k-1}+h_{k}\right)}{8}f_{ijk}\label{eq:two}\end{equation}


in which  

\[
h_{i}=x_{i+1}-x_{i},\: h_{j}=y_{j+1}-y_{j}\: h_{k}=z_{k+1}-z_{k}\]


The delta functions appearing in the right hand side of the starting
equations are approximated with linear B-splines (spl0) which spread
the point like charge along the nearest neighborhood. The resulting
$f_{ijk}$ represent the smearing of the point charges along the grid
points. 

For more details, including used unit system, please refer to the
Michel Holst's thesis and the APBS user guide online. To visualize
more clearly the problem, let's explicitly write the first equations
for a cubic grid of 5x5x5 containing general coefficients

\[
a_{222}u_{222}+a_{122}u_{122}+a_{322}u_{322}+a_{212}u_{212}+a_{232}u_{232}+a_{221}u_{221}+a_{223}u_{223}=f_{222}\]


\[
a_{322}u_{322}+a_{222}u_{222}+a_{422}u_{422}+a_{312}u_{312}+a_{332}u_{332}+a_{321}u_{321}+a_{323}u_{323}=f_{322}\]


\[
a_{422}u_{422}+a_{122}u_{122}+a_{322}u_{322}+a_{212}u_{212}+a_{232}u_{232}+a_{221}u_{221}+a_{223}u_{223}=f_{422}\]


\[
a_{232}u_{232}+a_{132}u_{132}+a_{332}u_{332}+a_{222}u_{222}+a_{242}u_{242}+a_{231}u_{231}+a_{233}u_{233}=f_{232}\]


\[
\ldots\]


in which the nodes are arranged using the natural ordering 

\[
U=[u_{111},u_{211},..,u_{N_{x}11,}u_{121},..,u_{221},u_{321},..,u_{N_{x}21}...,u_{N_{x}N_{y}N_{z}}]^{T}\]


Note that the prescribed values of nodes $u_{1jk},u_{N_{x},j,k},u_{i,1,k},u_{i,N_{y},k},u_{ij1}$
and $u_{ijN_{z}}$ along the faces of the box coming from the Dirichlet
boundary conditions will have their corresponding elements removed
in such a way that only equations for the interior nodes remain. In
other words, we will only consider the following set of unknown nodes

\[
U=[u_{222},u_{322},..,u_{N_{x}-1,22,}u_{232},..,u_{332},u_{432},..,u_{N_{x}-2,32}...,u_{N_{x}-1,N_{y}-1,N_{z}-1}]^{T}\]


in such a way that the previous equations become 

\[
a_{222}u_{222}+a_{322}u_{322}+a_{232}u_{232}+a_{223}u_{223}=f_{222}-a_{122}u_{122}-a_{212}u_{212}-a_{221}u_{221}\equiv b_{222}\]


\[
a_{322}u_{322}+a_{222}u_{222}+a_{422}u_{422}+a_{332}u_{332}+a_{323}u_{323}=f_{322}-a_{312}u_{312}-a_{321}u_{321}\equiv b_{322}\]


\[
a_{422}u_{422}+a_{322}u_{322}+a_{232}u_{232}+a_{223}u_{223}=f_{422}-a_{122}u_{122}-a_{212}u_{212}-a_{221}u_{221}\equiv b_{422}\]


\[
a_{232}u_{232}+a_{332}u_{332}+a_{222}u_{222}+a_{242}u_{242}+a_{233}u_{233}=f_{232}-a_{132}u_{132}-a_{231}u_{231}\equiv b_{232}\]


in which the boundary $u's$ are conveniently brought to the right-hand-side
of the equations. The resulting left-hand side equations can be written
in compact form in term of matrix vector product as follows

\[
Au=b\]


in which

\[
u\left(p\right)=u_{ijk},\qquad b\left(p\right)=b_{ijk},\qquad p=(k-2)(N_{x}-2)(N_{y}-2)+(j-2)(N_{x}-2)+i-1\]


\[
i=2,..,N_{x}-2,\quad j=2,..,N_{y}-2,\quad k=2,..,N_{z}-2\]


%
\begin{figure}
\includegraphics[scale=0.75]{Amatrix}

\caption{A Matrix representation}

\end{figure}


and $A$ is a (seven banded block tri-diagonal form) $(N_{x}-2)(N_{y}-2)(N_{z}-2)$
by $(N_{x}-2)(N_{y}-2)(N_{z}-2)$ squared symmetric positive definite
matrix containing the following nonzero elements (see figure 1):

\begin{itemize}
\item The main diagonal elements
\end{itemize}
\[
d_{0}(p)=\left[\epsilon_{i-1/2,j,k}^{x}\frac{\left(h_{j-1}+h_{j}\right)\left(h_{k-1}+h_{k}\right)}{4h_{i-1}}+\epsilon_{i+1/2,j,k}^{x}\frac{\left(h_{j-1}+h_{j}\right)\left(h_{k-1}+h_{k}\right)}{4h_{i}}+\right.\]


\[
\epsilon_{i,j-1/2,k}^{y}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{k-1}+h_{k}\right)}{4h_{j-1}}+\epsilon_{i,j+1/2,k}^{y}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{k-1}+h_{k}\right)}{4h_{j}}+\]


\[
\epsilon_{i,j,k-1/2}^{k}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{j-1}+h_{j}\right)}{4h_{k-1}}+\epsilon_{i,j,k+1/2}^{k}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{j-1}+h_{j}\right)}{4h_{k}}+\]


\begin{equation}
\left.\kappa_{ijk}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{j-1}+h_{j}\right)\left(h_{k-1}+h_{k}\right)}{8}\right]\label{eq:three}\end{equation}


\begin{itemize}
\item The Next upper band diagonal, which is shifted in one column to the
left from the first column, contains the following elements
\end{itemize}
\begin{equation}
\left[d_{1}(p)=-\epsilon_{i+1/2,j,k}^{x}\frac{\left(h_{j-1}+h_{j}\right)\left(h_{k-1}+h_{k}\right)}{4h_{i}}\right]\label{eq:four}\end{equation}


\begin{itemize}
\item The second upper band diagonal which is shifted $N_{x}-2$ columns
from the first column 
\end{itemize}
\begin{equation}
d_{2}(p)=\left[-\epsilon_{i,j+1/2,k}^{y}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{k-1}+h_{k}\right)}{4h_{j}}\right]\label{eq:five}\end{equation}


\begin{itemize}
\item The third upper band diagonal which is shifted $(N_{x}-2)(N_{y}-2)$
columns from the first column
\end{itemize}
\begin{equation}
d_{3}(p)=\left[-\epsilon_{i,j,k+1/2}^{k}\frac{\left(h_{i-1}+h_{i}\right)\left(h_{j-1}+h_{j}\right)}{4h_{k}}\right]\label{eq:six}\end{equation}


The remaining elements of the upper triangular squared matrix A are
set equal to zero. By symmetry we obtain the lower triagonal elements
of the matrix A. Because the matrix A is sparse and large, we can
implement efficient methods that optimally solve the linear system
for U. Specifically, we use the biconjugate gradient stabilized method
combined with the inexact LU decomposition of the matrix A. Having
the numerical values for the nodes in the interior of the box, we
finally add the previously removed prescribed values along the six
faces to get the solution over the complete set of grid points.


\subsection*{Computational algorithm}

\begin{enumerate}
\item We read the input file .inm to get the APBS input files (shifted dielectric
coefficients, kappa function and pqr data file) as well as the number
of grid points, the box Lengths, and temperature.
\item By using linear B-splines, we discretize the charge density to get
$f_{ijk}$ for $i=1,..,N_{x},\quad j=1,..,N_{y},\quad k=1,..,N_{z}$.
We also calculate the Dirichlet boundary condition along the six faces
of the box $u_{1jk},u_{N_{x},j,k},u_{i,1,k},u_{i,N_{y},k},u_{ij1}$
and $u_{ijN_{z}}$ using the temperature, the value of the bulk dielectric
coefficient (usually water) and ionic strength.
\item By using the expressions (\ref{eq:three}), (\ref{eq:four}),(\ref{eq:five}),
and (\ref{eq:six}), we evaluate the nonzero components of the matrix
$A$, e.g., the diagonal elements $d_{0}(p),d_{1}(p),d_{2}(p),$ and
$d_{3}(p),$ for $p=(k-2)(N_{x}-2)(N_{y}-2)+(j-2)(N_{x}-2)+i-1$ and
$i=2,..,N_{x}-1,\quad j=2,..,N_{y}-1,\quad k=2,..,N_{z}-1$. The values
for the shifted dielectric coefficients and kappa function elements
are obtained from the APBS input files. The values of the mesh size
$h_{i},h_{j}$ and $h_{k}$ are obtained from the number of grid points
and the Length of the box. Next, we built the sparse upper triangular
matrix A by filling with zeros the remaining elements of the matrix
A. Next, we obtain the lower triangular elements of the matrix A by
using the following symmetry property $A_{pq}=A_{qp}$ for $q=1,..,(N_{x}-2)(N_{y}-2)(N_{z}-2)$
and $p=q,..,(N_{x}-2)(N_{y}-2)(N_{z}-2)$. 
\item By using the values obtained for the discretized charge density $f_{ijk}$
and the values of the Dirichlet boundary elements multiplied by the
appropriate shifted dielectric coefficient values, we evaluate the
elements of $b_{ijk}$. We use the natural ordering $p=(k-2)(N_{x}-2)(N_{y}-2)+(j-2)(N_{x}-2)+i-1$
and $i=2,..,N_{x}-1,\quad j=2,..,N_{y}-1,\quad k=2,..,N_{z}-1$ to
construct the corresponding vector $b(p)$ (one index) from the data
array structure (three indices) $b_{ijk}$.
\item We use the inexact $LU$ decomposition of the matrix $A$. The default
tolerance value is set equal to 0.25 which provides a fast evaluation
of the matrices $L$ and $U$. 
\item The resulting $L$ and $U$ matrices, the matrix $A$ and the vector
$b$ are used to approximately solve $Au=b$ for the vector $u$ using
the biconjugate gradient stabilized method. The default accuracy is
set equal to 10\textasciicircum{}-9 and the maximum number of iteration
equal to 800. 
\item We use the natural ordering relationship to convert the resulting
vector $u(p)$ to data array structure to get the numerical solution
for $u_{ijk}$ for $i=2,..,N_{x}-1,\quad j=2,..,N_{y}-1,\quad k=2,..,N_{z}-1$.
\item Finally we add the previously removed values of the nodes at the faces
of the box to obtain the solution for the nodes $u_{ijk}$ over the
complete set of grid points, namely for $i=1,..,N_{x},\quad j=1,..,N_{y},\quad k=1,..,N_{z}$.
\item The electrostatic potential $u_{ijk}$ and the charge $f_{ijk}$ maps
are saved in dx format files.
\item The electrostatic potential surface $u_{ij(N_{z}+1)/2}$ is saved
in tiff and fig format files for visualization purpose. 
\end{enumerate}

\end{document}