diff --git a/doc/plnwv_API.tex b/doc/plnwv_API.tex
new file mode 100644
index 0000000000..2d4f57e44c
--- /dev/null
+++ b/doc/plnwv_API.tex
@@ -0,0 +1,522 @@
+\documentstyle{article}
+\begin{document}
+
+\title{API for plane-waves in NWChem}
+
+\author{Arthur Smith}
+
+\date{\today}
+\maketitle
+
+\section{General outline}
+There are two general features of plane-wave basis sets which justify
+a substantially different treatment from, for example, Gaussian basis
+sets. First, the basis set is not linked to the actual atomic positions,
+other than to the size of the periodic cell in which the atoms are contained,
+(in fact, for plane waves each individual basis element extends throughout
+this cell). Second, the basis sets tend to be very large while calculations
+involving a single basis element (overlaps for example) are essentially
+trivial, so the computational challenges differ considerably. For plane-wave
+density-functional theory with separable non-local pseudopotentials,
+all calculations can be accomplished in a time linear (or $N \log(N)$) in
+the basis-set size (there are further dependences on the number of atoms
+and orbitals, of course).
+
+Both of these features (independence of atomic positions, and large
+size with mostly trivial overlaps) can also be said to apply to
+other types of basis sets that we may wish to implement in the future: 
+pure real-space grid methods (the plane-wave technique already does
+much of the computation on a real-space grid), multi-grid methods,
+mixed-basis techniques, or even wavelets. I believe the API should
+be designed in such a way that programmers will need to make only minimal
+modifications to a higher-level application code in order to utilize
+a new basis set when/if one becomes available.
+
+\section{Pre-requisites}
+
+Plane-wave basis sets are only directly applicable to periodic systems -
+solids in particular, although other systems can be readily treated
+by super-cell methods. Thus the application must have previously
+specified a periodically repeated cell of some sort.
+
+In addition,
+particularly for relative small periodic cells, a collection of
+special {\bf k}-points must be specified (selected according to the
+cell symmetry) which will be used to approximate the integration of
+quantities over the Brillouin zone of the cell with a weighted sum
+over these {\bf k}-points.
+
+For large systems or for molecules/clusters treated in the supercell
+approximation, it is often appropriate to just use the single point
+{\bf k = 0}. In any case, the code should assume the {\bf k}-points
+have been previously specified.
+
+\section{Handling the basis sets}
+
+There is then a slightly different plane-wave basis set for each
+{\bf k}-point, and another (usually larger one) for the charge density.
+These are specified by the dimensions of the periodic
+cell of the system, and by some sort of cutoff - generally the requirement
+that:
+\begin{equation}
+\hbar^2|{\bf k + G}|^2/2m < E_{cut}
+\end{equation}
+for the wavefunction basis sets and:
+\begin{equation}
+\hbar^2 |{\bf G}|^2/2m < d \cdot E_{cut}
+\end{equation}
+for the charge density, where the parameter $d$ is between 1 and 4 and
+determines the degree to which the charge density completely reproduces 
+all oscillations in the wavefunctions.
+
+The reciprocal lattice vectors {\bf G} are simply integer linear combinations
+\begin{equation}
+{\bf G}^{\alpha} = n^{\alpha}_1 {\bf b}_1 + n^{\alpha}_2 {\bf b}_2 +
+			n^{\alpha}_3 {\bf b}_3
+\end{equation}
+of the fundamental reciprocal lattice vectors ${\bf b}_i$, which satisfy
+\begin{equation}
+{\bf b}_i \cdot {\bf a}_j = 2 \pi \delta_{ij}
+\end{equation}
+and where the ${\bf a}_j$ are the fundamental lattice vectors specifying
+the periodic cell of the system.
+
+The {\bf G} values used (and the corresponding $n^{\alpha}_i$ triplets)
+are specified by spherical (ellipsoidal) regions surrounding the origin for
+the charge density and {\bf k = 0} wavefunctions, or by similar regions
+surrounding points close to the origin for ${\bf k} \neq {\bf 0}$.
+
+In order to apply Fast-Fourier-Transform (FFT) techniques to allow
+transformations between the plane-wave basis set and a real-space grid on
+which charge density and potential calculations are performed,
+the triplets $n^{\alpha}_i$ must be mapped to corresponding points in
+a three-dimensional integer rectangular grid $N_1 \times N_2 \times N_3$.
+Each of the $N_i$ must be greater than twice the largest absolute value
+of all the $n_i$, and for most library FFTs must be a power of two. After
+applying the FFT, the resulting data is associated with a real-space grid
+at the points:
+\begin{equation}
+{\bf r}_{ijk} = {i \over N_1} {\bf a}_1 + {j \over N_2} {\bf a}_2 +
+		{k \over N_3} {\bf a}_3
+\end{equation}
+
+Note that the point i, j, k in the $N_1 \times N_2 \times N_3$ reciprocal
+space array is associated with a triplet $n^{\alpha}_l$ close to the
+origin as follows:
+\begin{equation}
+{\rm if\ \ } i < N_1/2 {\rm\ \  then\ \ } n^{\alpha}_1 = i
+\end{equation}
+\begin{equation}
+{\rm if\ \ } i > N_1/2 {\rm\ \  then\ \ } n^{\alpha}_1 = i - N_1
+\end{equation}
+and similarly for j and k.
+
+
+\subsection{Data distribution}
+Utilizing these FFT methods, the application of the Hamiltonian to
+a wavefunction involves diagonal operators or products of diagonal
+operators (see the section on applying H to $\psi_i$ below). One
+nice aspect of a diagonal operator on a parallel computer is that
+however the vectors involved are distributed on the machine, application
+of the operator involves no communication (as long as the components
+of the operator itself are either distributed in the same way as the
+vectors, or are calculated when needed). In other words, calculation of
+\begin{equation}
+	g_i ({\bf r}) = f({\bf r}) \psi_i({\bf r})
+\end{equation}
+is a local operation, no matter where the components $\psi_i({\bf r})$
+actually reside.
+
+However, communication {\it is}\ required during the FFT operations
+if $\psi_i({\bf r})/\psi_i({\bf G})$ is distributed in {\bf r}/{\bf G},
+and extensive communication
+is also required during orthogonalization and similar inner-product
+operations if $\psi_i({\bf r})/\psi_i({\bf G})$ is distributed in the
+index $i$. In addition, since we are thinking of data-decomposition and not
+task-decomposition here, steps must be taken to ensure that every
+process carries roughly the same number of coefficients $\psi_i({\bf G})$,
+which is best accomplished by some form of interleaving. Note that
+$\psi_i({\bf r})$ will naturally be evenly distributed in {\bf r} under most
+decompositions because the {\bf r} lie on a rectangular grid, while
+the relevant {\bf G} values lie within a sphere, for which block
+decompositions may leave some processors with few or no {\bf G}-vectors.
+
+A further complication in the {\bf k = 0} situation is that the
+coefficients for {\bf G} and {\bf -G} are complex conjugates, and therefore
+nearly a factor of two in memory and speed can be gained by keeping only
+one of the two coefficients. But that means in the data decomposition,
+{\bf G} and {\bf -G} should be covered by the same process, requiring
+a reversed interleaving for negative values of the $n_i$.
+
+\subsection{The initialization API}
+Thus the fundamental information needed to initialize a particular plane-wave
+basis set is:
+\begin{enumerate}
+\item The periodic cell vectors ${\bf a}_i$
+\item The {\bf k}-vectors used.
+\item The energy cutoff $E_{cut}$
+\item The charge density basis parameter $d$.
+\item Some information about number of processes, preferred data
+decomposition, and possibly process communication topology, in order
+to determine the distribution of G and R vectors among nodes.
+\item Preferred values of $N_1, N_2, N_3$ (optional)
+\item The vectors ${\bf b}_i$ (optional - might be useful
+for modularization).
+\end{enumerate}
+
+Returned from the initialization will be:
+\begin{enumerate}
+\item The size of the basis set for each k-point, and for the charge
+density, including total size, and number allocated to this process.
+\item Some form of listing of the G-vectors: at least the $n^{\alpha}_i$,
+preferably also the actual vectors {\bf G} and their (squared?) magnitudes.
+\item The output values of $N_1, N_2, N_3$ for the whole space and
+for the sub-blocks (or other arrangement) allocated to this process.
+\item Optionally, the ${\bf b}_i$ vectors.
+\end{enumerate}
+
+\subsection{The allocation API}
+Typically, we'll be allocating some number $n_{orbitals}$ of
+complex coefficients of the $ng({\bf k})$ plane waves for each of
+the {\bf k}. In addition, there will be several auxiliary arrays
+individually of approximate size $ng$, including arrays to handle
+the charge density and components of the potential. And we also need
+a few complex and real arrays on the $N_1\times N_2 \times N_3$
+real-space grid.
+
+These would presumably be allocated with GA routines or something
+similar - for example in our hybrid decomposition approach, the first
+set of coefficient arrays would be distributed across all the processors. The
+second set, of auxiliary arrays, would be distributed across the processors
+within a group (the spatial decomposition subgroups) but duplicated (sometimes
+with identical data) in all such groups, in order to eliminate all
+communications except during FFT and orthogonalization (and global operation)
+routines.
+
+Inputs then to an allocator would be:
+\begin{enumerate}
+\item A tag specifying the basis set (previously initialized).
+\item The data type of the coefficients (real or complex, or something else?)
+\item The number of ($ng$ arrays) of this type to be
+allocated ($n_{orbitals}$ for example.
+\item Optionally: A flag indicating whether the data should be
+initialized with random values or with zeroes.
+\end{enumerate}
+
+Output would be a tag specifying the coefficient matrix, similar to
+GA matrix allocation, maybe?
+
+\subsection{API for basis set operations}
+
+Once we have the basis set specifications, there are a number of
+operations we'll need to do anything with them:
+\begin{enumerate}
+\item Copy
+\item Scale or AXPY
+\item application of a diagonal operator (GBMV?)
+\item Inner product
+\item Send and receive, or get and put to a different process
+\item Transformation between basis sets
+\item Computation of basis set overlaps for different basis sets
+\end{enumerate}
+
+Most of these are straightforward: probably only the last two
+will require much thought. We'll want to do more complicated
+things with the charge density (at least in calculation of
+the exchange-correlation energy), but that is probably a special
+case that can be dealt with separately.
+
+There are a couple of transformations between basis sets we
+might want to consider implementing initially: within a G-space
+representation, changing form one cutoff energy to a different
+one is trivial in principle (just chop, or zero-pad) but
+not necessarily in practice where reading in a starting configuration
+from an old basis set may misplace coefficients if you are not
+careful (we have not yet successfully implemented this using
+our current code, for example). Transformations between G-space and
+real-space are essential to computation of the charge density and
+evaluation of the real-space potentials, and this can be viewed
+as one kind of general basis-set transformation (implemented
+through fast Fourier transforms).
+
+I think computation of overlaps will not be necessary until
+mixed Gaussian and plane-wave bases are being considered, but
+it might be worth-while thinking about how to do it now.
+
+\section{Applying H to $\psi_i$}
+
+We can use various minimization techniques to isolate the
+lowest $n_{occupied}$ eigenstates of H. Because $ng$ is so much
+larger, we never want to construct the full $ng\times ng$ matrix,
+so most of these techniques start from the residual
+\begin{equation}
+H \psi_i - \epsilon_i \psi_i.
+\end{equation}
+and the most computationally intensive part of the code is
+going to be evaluation of the Hamiltonian-vector product H$\psi_i$.
+
+\subsection{Kinetic energy}
+
+The kinetic energy operator is evaluated in the reciprocal space
+representation:
+\begin{equation}
+T \psi_i({\bf G}) = {\hbar^2|{\bf k + G}|^2 \over 2m} \psi_i({\bf G})
+\end{equation}
+
+The values of $|{\bf k + G}|^2$ should be pre-stored (possibly
+returned from the initialization routine above) and so this is
+a straightforward operation. However, since this is peculiar to
+the plane-wave basis, it might be a good idea to define a kinetic
+energy API (for example, for pure real-space grids, the kinetic energy
+could be evaluated using finite-difference techniques). For such
+an API the input is a basis-set specifier and the coefficient matrix
+(with an integer specifying how many orbitals to work on at once)
+and output is the matrix of result vectors (placed in a work area provided).
+
+\subsection{The real-space potentials}
+
+These are comprised of the local pseudopotential, the Hartree
+(electron-electron Coulomb) potential and  the exchange-correlation potential.
+The potentials themselves are evaluated as follows:
+\begin{equation}
+V_{local}({\bf G}) = \sum_{is} u_{is}({\bf G}) S_{is}({\bf G})
+\end{equation}
+where $is$ is an index counting the atomic species in the system,
+$u_{is}$ is the potential from a single atom of this species, and
+$S_{is}$ is a structure factor for the species:
+\begin{equation}
+S_{is}({\bf G}) = \sum_{ia} \exp(i {\bf G}\cdot {\bf r}_{ia})
+\end{equation}
+where $ia$ runs over all atoms of species $is$. $u_{is}$ and $S_{is}$
+would be distributed in the same fashion as the {\bf G}'s, but with
+all species $is$ together (not distributed), and so duplicated
+in the same fashion as $V_{local}({\bf G})$ itself: among processes with the
+same {\bf G}'s but different collections of orbitals. So this computation
+involves no communication, except for the FFT required in the end
+to produce a real-space potential.
+
+The electron-electron direct Coulomb interaction is easily evaluated:
+\begin{equation}
+V_{Hartree}({\bf G}) = 4\pi \rho(G)/G^2
+\end{equation}
+except for ${\bf G} = {\bf 0}$ where the corresponding term is
+a constant that is added to the constant electrostatic energy of the
+periodic system usually evaluated by Ewald summation. Similarly, the
+only communication required here is the final FFT, although there
+will be duplication if the orbitals are distributed.
+
+Both $V_{local}$ and $V_{Hartree}$ are then FFT'd to real space, and
+can be added to the exchange correlation potential:
+\begin{equation}
+V_{xc}({\bf r}) = V_{xc}(\rho({\bf r}), \nabla \rho({\bf r}), ...)
+\end{equation}
+which requires no communication at all, except in the
+evaluation of the charge density and gradients:
+
+The charge density $\rho({\bf r})$ is calculated as the sum
+\begin{equation}
+\rho({\bf r}) = \sum_{i} |\psi_i({\bf r})|^2
+\end{equation}
+and the gradient is:
+\begin{equation}
+\nabla\rho({\bf r}) = FFT \bigl\{ i {\bf G} \rho({\bf G}) \bigr\}
+\end{equation}
+and similarly for higher derivatives.
+
+The duplication discussed could be eliminated at the expense of
+more communication: we should probably study this at some point. In
+particular since the gradient-dependent exchange correlation potential 
+cannot be done by table-lookup (because it depends on more than one
+parameter, unlike the LDA) the computation involved may be sufficient
+to justify distributing it.
+
+An API for a set of subroutine calls to evaluate these potentials
+seems appropriate - what is needed on input is quite a variety of
+things: for the charge density calculation we probably want to
+pass the basis set information and the G-space coefficient matrix,
+and the output would be either the charge density on the real-space
+grid, or else the fitted charge-density in G-space (or both).
+For the local potential we need to pass in geometry information and
+information about the pseudopotential, as well as the basis set - output
+is the G-space representation of that potential. For the Hartree potential
+we pass in the charge density and basis set info, output is the G-space
+potential. For the gradients calculation pass in the G-space charge
+density, return the real-space charge density and gradients to be
+passed to an exchange-correlation subroutine, or else just return the
+real-space exchange-correlation potential.
+
+We also need an initialization API for the pseudopotential to
+generate the $u_{is}({\bf G})$, and
+for the exchange-correlation potentials (at least to allow
+choices between the possibilities, and also possibly to allow
+setting some of their parameters). More on the pseudopotential below...
+
+\subsection{The non-local pseudopotential}
+Using the Kleinman-Bylander separable decomposition, the non-local
+pseudopotential (for each species $is$!) can be evaluated in
+reciprocal space as follows (from memory, so bear with me if I
+missed a minus sign or something):
+\begin{eqnarray}
+V_{nl}^{lm} \psi_i ({\bf G}) = C^l
+		&\sum_{ia} U_{nl}^l({\bf k+G})Y_{lm}^*({\bf k + G})
+		\exp(i {\bf G}\cdot{\bf r}_{ia}) \cdot \cr
+		&\sum_{{\bf G}'} U_{nl}^l({\bf k + G}')Y_{lm}({\bf k + G}')
+		\exp(-i {\bf G}'\cdot{\bf r}_{ia})
+		\psi_i({\bf G}')\cr
+\end{eqnarray}
+where $ia$ is an index over all the atoms of species $is$, and $l$, $m$ are
+standard angular momentum quantum numbers. The $U_{nl}^l$ and $Y_{lm}$
+can be computed at initialization and are $O(N)$ in size so not too bad.
+This would then be summed over $l$, $m$, and species $is$.
+
+One big potential problem data and
+communication-wise here is the $\exp(i {\bf G}\cdot {\bf r}_{ia})$ matrix,
+which if precomputed and stored is $O(N^2)$ in size, and yet needs to be present
+for {\bf all} the orbitals (note that non-local calculation is an $O(N^3)$
+operation and thus one of the major time-consuming parts of the computation
+for large systems). What we do is evaluate the matrix only for {\bf G}
+vectors along the principal axes ${\bf b}_i$, and then compute the product
+of three complex numbers to get the final value when needed - this means
+we need to only store a matrix $O(N^{4/3})$ in size, but we need to
+ensure the correct sort of data-locality for that matrix.
+
+There are other ways of doing the nonlocal pseudopotentials: real-space only
+evaluations for example, or taking advantage of possible redundancy if a part
+of the system has a shorter periodicity (for example, the bulk portion of
+a surface problem) or other symmetry, or using complicated hybrids like
+the ultra-soft Vanderbilt pseudopotentials, and presumably other methods
+will come along.
+
+So I think two (or possibly three) subroutine calls would be involved
+that could be standardized on: initialization of the pseudopotential
+(calculation of the $U_{nl}({\bf k+ G})$ and local $u({\bf G})$ values)
+and of other needed quantities ($Y_{lm}({\bf k + G})$ and the
+$\exp(i {\bf G}\cdot {\bf r}_{ia})$ matrix); and application of
+the pseudopotential to a set of orbitals (either with an ion force
+calculation, or with forces calculated separately).
+
+\subsection{Dealing with the pseudopotentials}
+
+There may be several types of pseudopotentials we would want to
+handle - for example, it might be worth testing the Gaussian-derived
+ECP's with plane waves (maybe not?) Anyway, we must initially
+read in some representation of the pseudopotentials about each species
+of atom in the system, and the input format we currently use most
+often is a simple radial grid. These then need to be converted to
+the local $u({\bf G})$ and nonlocal $U_{nl}({\bf k+G})$ forms needed in
+the actual application of the pseudopotentials.
+
+Thus there needs to be some specification of the input representation
+of the pseudopotentials, and an initialization routine for converting
+that representation to one appropriate for the basis set. The associated
+API would have as input the input representation and specification
+of the (plane-wave) basis set, and as output the converted pseudopotential
+values.
+
+\section{Achieving self-consistency}
+
+With plane-wave basis sets the central problem of converging to
+a self-consistent charge-density (represented by a collection of
+occupied Kohn-Sham eigenfunctions) is fundamentally the same as for
+density functional theory with other basis sets. However, the principal
+constraint imposed by the large size of the plane-wave basis set is
+that we never know the full matrix $H$, only matrix-vector products
+$H \psi_i$, meaning that direct diagonalization cannot be performed,
+and instead iterative approaches are needed. In addition, of course,
+$H$ itself changes with each iteration because the charge density
+is changing. Among a wide variety of possibilities there are basically
+two iterative schemes that have become adopted in plane-wave basis-set
+methods: Car-Parrinello ``molecular-dynamics'' type iterations, and
+Teter-Payne-Allen preconditioned gradient iterations (often
+misnamed ``conjugate gradient'' iterations because the first implementation
+happened to use a multiple conjugate gradients updating scheme that
+has since proved inefficient compared to just using a single preconditioned
+gradient).
+
+The Car-Parrinello approach is to treat the energy gradient $H\psi$ as
+a force on the coefficients of $\psi$, so that in a ``time-step'' $\tau$,
+\begin{equation}
+\psi_i(t + \tau) = \psi_i(\tau) - {\tau^2 \over 2 m_{eff}} H\psi_i(tau),
+\end{equation}
+after which the new $\psi_i$ are orthogonalized to one another.
+When doing real dynamics a velocity term may also be included - as stated
+this is effectively just a steepest-descents algorithm. As such it
+is simple but rather inefficient for getting self-consistency, but
+on the other hand it is rather robust (always converging to the lowest
+occupied subspace, if you wait long enough) and in fact when coupled
+with atomic dynamics can be competitive with the Born-Oppenheimer
+approach required when using the preconditioned gradient methods.
+
+The preconditioned gradient methods have the following iterative step:
+\begin{enumerate}
+\item Generate residual vectors: $ R_i = (H - \epsilon_i) \psi_i$
+\item Apply a preconditioner: $B_i = P R_i$. As with the
+Davidson method, the preconditioner should be close to an inverse
+of the Hamiltonian.  In practice, a rational function of the kinetic
+energy (approaching 1/K.E. for high energy) works well with a
+plane-wave basis set.
+\item Orthogonalize the $B$'s to one another and to the $\psi's$.
+\item Generate $ H B_i $ and the subspace matrix  $ <V_i|H|V_j>$ where
+$V_i$ is $\psi_i$ for $i < n$ and $B_{i - n}$ for $i >= n$.
+\item Diagonalize within the subspace, and let the new $\psi$ vectors
+be the lowest $n$ eigenvectors.
+\item Generate the new charge density (generally by mixing in some
+of the old, otherwise instabilities arise) and the resulting new
+Hamiltonian, and repeat.
+\end{enumerate}
+These methods (there are variants based on whether you let $n = 1$ and
+have an outer loop over bands, or else let $n =$ all states, or else
+with different preconditioners) tend to be much faster than the
+steepest-descent approach to obtain a self-consistent set of electronic
+states. On the other hand they require more computation per step
+(twice as many Hamiltonian-vector operations, for example) and experience
+shows they tend to be less robust.
+
+Both methods could be implemented with no knowledge of the
+underlying basis sets or data decomposition, as long as the
+following subroutine calls were available:
+\begin{enumerate}
+\item Application of $H$ to $\psi$-type vectors
+\item Application of a preconditioner to $\psi$-type vectors
+\item AXPY operations on $\psi$-type vectors and $\rho$-type vectors
+\item Inner products on $\psi$-type vectors (at least for
+generating the $<V|H|V>$ matrix, possibly also for orthogonalization unless
+that was a separate subroutine call)
+\item Generating $\rho$ from $\psi_i$, and modifying $V_{xc}$ and $V_{Hartree}$
+accordingly, to complete the self-consistent loop.
+\end{enumerate}
+In addition, a diagonalization routine for $2 n\times 2 n$ matrices is
+needed, and probably should be done in parallel in accordance with
+the orbital decomposition of the data.
+
+\section{Energy and forces}
+
+At the highest level, an API that calls on the electronic
+energy minimization to generate the energy and forces for a particular
+configuration of atoms might be workable, and would certainly be
+useful. This would facilitate rapid comparison of different methods
+for a particular problem, in addition to providing a natural interface
+for standard molecular dynamics codes.
+
+Input to such a high level API would include:
+\begin{enumerate}
+\item Geometry - atomic configuration
+\item Basis set, pseudopotential specifications (data decomp?)
+\item (Optionally) Trial set of orbitals
+\item Convergence parameters: iteration method, maximum number of
+iterations, tolerance on the final numbers
+\end{enumerate}
+
+Output would of course be energy and forces, as well as
+(optionally) a final set of orbitals for restarting, and a flag
+specifying whether the convergence tolerance was actually reached.
+
+In addition, I think it would be important to have callable
+routines that estimate (or at least keep track of) memory usage
+and time per iteration, to allow analysis of resource utilization
+by the calling program. To do the estimate means having some kind
+of performance model of the code - we do already have such a thing
+for the memory usage, but have never tried it for timing...
+
+\end{document}
+\bye