Higher Order Response in ${\cal O}(N)$ by Perturbed Projection
Val\'ery Weber
Perturbed projection for linear scaling solution of the coupled-perturbed self-consistent-field equations [Weber, Niklasson and Challacombe, Phys.\ Rev.\ Lett. {92}, 193002 (2004)] is extended to the computation of higher order static response properties. Although generally applicable, perturbed projection is developed here in the context of the self-consistent first and second electric hyperpolarizabilities of three dimensional water clusters at the Hartree-Fock level of theory. Non-orthogonal, density matrix analogues of Wigner's $2n+1$ rule are given up to fourth order. Linear scaling and locality of the higher order response densities under perturbation by a global electric field are demonstrated.
Response3.tex
% NOTE TO AIP TYPSETERS: TO CONVERT FROM TWO-COL TO PREPRINT, SWITCH
% COMMENTOUT COMMAND FROM A TO B IE. USE:
%
 \documentclass[prl,twocolumn,showpacs,twocolumngrid,superbib]{revtex4}
 \newcommand{\commentoutB}[1]{}
 \newcommand{\commentoutA}[1]{#1}
%
% INSTEAD OF THE FOLLOWING:
%
%\newcommand{\commentoutB}[1]{#1}
%\newcommand{\commentoutA}[1]{}
%\documentclass[prl,aps,preprint,showpacs,superbib]{revtex4}
%
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{bm}
\usepackage{alltt}
\usepackage{fancyhdr}

\pagestyle{fancy}

\def\Tr{{\rm Tr}}
\def\F{\mathcal{F}}
\def\D{\mathcal{D}}
\def\X{\mathcal{X}}
\def\E{\mathcal{E}}

\begin{document}

\title{Higher Order Response in ${\cal O}(N)$ by Perturbed Projection} 

\author{Val\'ery Weber}
\email{valery.weber@unifr.ch}
\affiliation{Department of Chemistry, University of Fribourg, 1700 Fribourg, Switzerland.}
\author{Anders M. N. Niklasson}
\author{Matt Challacombe}
\affiliation{Los Alamos National Laboratory, Theoretical Division, Los Alamos 87545, New Mexico, USA.}

\date{\today}

\begin{abstract}
Perturbed projection  for linear scaling solution of the coupled-perturbed 
self-consistent-field equations 
[Weber, Niklasson and  Challacombe, Phys.\ Rev.\ Lett. {\bf 92}, 193002 (2004)] 
is extended to the computation of higher order static response properties.
Although generally applicable, perturbed projection is developed here 
in the context of the self-consistent first and second electric hyperpolarizabilities 
of three dimensional water clusters at the Hartree-Fock level of theory. 
Non-orthogonal, density matrix analogues of Wigner's $2n+1$ rule are given up to fourth order.
Linear scaling and locality of the higher order response densities under perturbation 
by a global electric field are demonstrated.  
\end{abstract}

\keywords{Perturbed Projection, density matrix, perturbation theory, purification, linear scaling, 
          CPHF, Hartree-Fock, electric polarization, hyperpolarizabilty}

\pacs{02.70.-c, 71.15.Dx,31.15.Ar,31.15.Md,31.15.Ne,33.15.Kr,36.40.Cg}

\maketitle

\footnotetext[1]{LA-UR-04-5219} 

\section{Introduction}
First principles electronic structure theory has traditionally been limited 
to the study of small systems with a limited number of nonequivalent atoms. 
Despite the tremendous increase in computational power of digital computers this 
has remained the case, until the advent of reduced complexity algorithms over the
last decade \cite{GGalli96,DBowler97,SGoedecker99,POrdejon00,VGogonea01,SWu02}. In the 
best case, these reduced complexity algorithms scale only linearly with system size, $N$, 
allowing simulation capabilities to keep pace with hardware improvements.
These linear scaling algorithms exploit the quantum locality (or nearsightedness) of 
non-metallic systems,  manifested in the approximate exponential decay of density matrix elements 
with atom-atom separation through the effective use of sparse matrix methods. For small systems,
linear scaling methods may be inefficient due to overhead.  However, for large, complex systems
these methods hold the promise of major impact across materials science, chemistry and biology. 

So far, a majority of work in linear scaling electronic structure theory 
has focused on methods and calculations involving the ground state, with little 
attention devoted to the problem of response properties.  The calculation of static
response within  Hartree-Fock or Density Functional Theory may be  obtained through solution 
of the Coupled-Perturbed Self-Consistent-Field (CPSCF) equations, which yield properties 
such as the electric polarizability and hyperpolarizability \cite{HSekino86,SKarna91}, 
the Born-effective charge, the nuclear magnetic shielding tensor \cite{KWolinski90}, 
indirect spin-spin coupling
constant \cite{CPennington91,OMalkin96}, geometric derivatives (i.e.~higher order
analytic force constants) \cite{RAmos89} and polarizability derivatives such as the Raman 
intensity \cite{MLazzeri03,OQuinet01}, to name but a few.

Conventional approaches to solution of the CPSCF equations 
\cite{HSekino86,SKarna91,JPople79} are based on perturbation of the wave 
function, requiring an $N^3$-scaling eigensolve which may need to be followed by an ${\cal O}(N^5)$ 
transformation of two-electron integrals, depending on the method. 
In addition to the formal scaling of these conventional methods, they do not admit exploitation 
of quantum locality through the effective use of sparse matrix algebra.  
More recently, schemes with the potential for reduced complexity have been been put forward.
Ochsenfeld and Head-Gordon proposed a scheme based on the Li-Nunes-Vanderbilt 
density matrix minimization \cite{Ochsenfeld97}.  Later, Larsen {\em et al.} \cite{HLarsen01a} 
proposed iterative solution of the CPSCF equations involving equations derived from 
unitary operations and approximations to the matrix exponential.    In both of these approaches, 
a linear system of equations containing commutation relations is obtained, which {\em implicitly}
determines the response function.  However, the method of solution for these equations
is not discussed, and computational results are not presented.  Recently though, 
with an apparent reformulation of Ref.~[\onlinecite{Ochsenfeld97}],
Ochsenfeld, Kussmann and Koziol~\cite{COchsenfeld04} have achieved near linear scaling computation
of NMR chemical shifts for one-dimensional alcanes at the GIAO-HF/6-31G* level of theory,
but likewise provide no details on their method of solution.  These implicit commutation relations are Sylvester-like
and may be solved with a number of approaches \cite{JBrandts01}, the particulars of which are of interest. 

In contrast, Perturbed Projection \cite{VWeber04} is a recently developed alternative 
for $N$-scaling solution of the CPSCF equations that is simple and explicit.  Based on a recently developed density 
matrix perturbation theory \cite{ANiklasson04}, Perturbed Projection exploits the explicit 
relationship between the density matrix $\mathcal{D}$ and the effective Hamiltonian or 
Fockian $\mathcal{F}$ via spectral projection; $\mathcal{D}=\theta(\tilde{\mu}I-\mathcal{F})$, 
wherein $\theta$ is the Heaviside step function (spectral projector)  and the chemical potential 
$\tilde{\mu}$ determines occupied states via Aufbau filling.   Spectral projection can be carried 
out in a number of ways \cite{ANiklasson02A,ANiklasson03,RMcWeeny60,WClinton69,APalser98,GBeylkin99,KNemeth00,AHolas01}.
Of special interest here are recursive polynomial expansions of the projector, including the second 
order trace correcting (TC2) \cite{ANiklasson02A} and fourth order trace resetting (TRS4) 
\cite{ANiklasson03} purification algorithms.  These new methods (TC2 and TRS4) have convergence 
properties that depend only weakly on the band gap, do not require knowledge of the chemical potential 
and perform well for all occupation to state ratios. 
Perhaps most important to the current contribution, these methods converge rapidly to smooth, 
monotone projectors.

Prior to Ref.~[\onlinecite{COchsenfeld04}], Perturbed Projection demonstrated linear scaling in 
computation of the first electric polarizability for three-dimensional water clusters with the 
Hartree-Fock model \cite{VWeber04}. Also, in a preceding paper, we outlined a non-orthogonal density matrix
perturbation theory \cite{ANiklasson05a} for response to a change in basis (i.e.~as occurs in 
the evaluation of higher order geometric energy derivatives \cite{RAmos89}). 
In this article, the Perturbed Projection method is extended to higher orders
in the electric polarizability, up to fourth order in the total energy.  

This paper is organized as follows: First we describe the perturbation expansion and 
the computation of response properties through solution of the CPSCF equations.
Then we present extension of Perturbed Projection through higher orders and the 
computation of properties using a density matrix analogue of Wigner's $2 n+1$ rule.
Next,  we present several examples of calculated higher order response properties.
We show the saturation of hyperpolarizabilities up to  fourth order (i.e. up to the 
second hyperpolarizability $\gamma$) for a series water chains.
We also  demonstrate linear scaling complexity for the solution of the higher order CPSCF
equations and an approximate exponential decay in elements of higher order response functions
for 3D water clusters.  Finally, we discuss these results and present our conclusions. 

\section{The Coupled Perturbed Self-Consistent-Field Equations}

The Coupled-Perturbed Self-Consistent-Field (CPSCF) equations yield
static response functions and properties in models including both the 
Hartree-Fock (HF) and Density Functional Theory (DFT).  In the following
we develop Perturbed Projection for solution of the CPSCF equations in the framework of polarization and 
Hartree-Fock theory.  In many cases, the extension of Perturbed Projection to
the computation of other static perturbations is straightforward.  In
the case of model chemistries that involve DFT, an extra programming effort is required \cite{Lee_1994,PSalek02}.
Also, in the case of properties such as the NMR chemical shift and geometric derivatives 
(force constants), perturbation of the non-orthogonal basis requires additional 
considerations that we have detailed in a preceding paper \cite{ANiklasson05a}.

\subsection{Notation}

Superscripts and subscripts refer to perturbation order and 
self-consistent cycle count respectively. The symbols $\mathcal{D},\mathcal{F},\dots$
are matrices in an orthogonal representation, while
$D,F,\dots$ are the corresponding matrices in a non-orthogonal basis.
The transformation between orthogonal and non-orthogonal 
representations is carried out in ${\cal O}(N)$ using
congruence transformations \cite{JWilkinson65,GStewart73} provided 
by the approximate inverse (AINV) algorithm for computing  sparse 
approximate inverse Cholesky factors with a computational complexity
scaling linearly with the system size \cite{MBenzi95,MBenzi96,MBenzi01}.

\subsection{Response expansions}

Within HF theory, the total electronic energy $E_{\rm tot}$ of 
a molecule in a static electric field $\mathcal{E}$ is
\begin{equation}
  \begin{split}\label{totalE}
   E_{\rm tot}(\E)&=Tr \{ D \cdot (h^0+\mu\E) \}+\frac{1}{2}Tr \{ D \cdot (J[D]+K[D]) \} \\
                  &=Tr \{ D \cdot F[D] \}-\frac{1}{2}Tr \{D \cdot (J[D]+K[D]) \},
  \end{split}
\end{equation}
where $D \equiv D[\E]$ is the density matrix in the electric field $\mathcal{E}$, 
$h^0$ is the core Hamiltonian, $\mu$ is the dipole moment matrix, 
$J[D]$ is the Coulomb matrix, $K[D]$ the exact HF exchange matrix
and 
\begin{equation}
F \equiv F[\E]=h^0+\mu\E+J[D(\E)]+K[D(\E)]
\end{equation}
is the Fockian.
The total energy of a molecule in a homogeneous electric field may 
be developed in a Taylor expansion series around $\E = 0$ as
\begin{equation}
  \begin{split}
    E_{\rm tot}(\E)= E_{\rm tot}(0) 
    &-\sum_a\mu_a\E^a\\
    &-\frac{1}{2!}\sum_{ab}\alpha_{ab}\E^a\E^b\\
    &-\frac{1}{3!}\sum_{abc}\beta_{abc}\E^a\E^b\E^c\\
    &-\frac{1}{4!}\sum_{abcd}\gamma_{abcd}\E^a\E^b\E^c\E^d\\
    &+\dots,
  \end{split}
\end{equation}
 where $\alpha_{ab}$ is the polarizability, and $\beta_{abc}$ and 
 $\gamma_{abcd}$ are the first and second 
 hyperpolarizabilities, respectively, $\mu_a$ is the dipole 
 moment, and $\E^a$ is the electric field in direction $a$. 
 The polarizability $\alpha_{ab}$ is the second order response 
 of the total energy with respect to variation in the electric field 
 while the higher derivatives, $\beta_{abc}$ and $\gamma_{abcd}$, give 
 rise to the first and second hyperpolarizabilities \cite{HSekino86,SKarna91} where 
 \begin{subequations}\label{pol}
   \begin{gather}
     \alpha_{ab}=
     -\frac{\partial^2 E_{\rm tot}}{\partial \mathcal{E}^a\partial \mathcal{E}^b}
     \bigg|_{\mathcal{E}=0}=
     -2Tr[D^a\mu_b],\\
     \beta_{abc}=
     -\frac{\partial^3 E_{\rm tot}}{\partial \mathcal{E}^a\partial \mathcal{E}^b\partial \mathcal{E}^c}
     \bigg|_{\mathcal{E}=0}=
     -4Tr[D^{ab}\mu_c],\\
%     -2Tr[D^{ab}\mu_c],\\
     \gamma_{abcd}=
     -\frac{\partial^4 E_{\rm tot}}{\partial\mathcal{E}^a\partial\mathcal{E}^b\partial \mathcal{E}^c\partial \mathcal{E}^d}
     \bigg|_{\mathcal{E}=0}=
     -12Tr[D^{abc}\mu_d].
%     -2Tr[D^{abc}\mu_d].
   \end{gather}
\end{subequations}
Here $D^{a\ldots}$ denotes a density matrix derivative with respect to a field in directions $a\ldots$ 
at $\mathcal{E} = 0$.  The  density matrix derivative or ``response function'' is given by 
\begin{equation}
% \displaystyle\D^{a\ldots}=N!
 \displaystyle\D^{a\ldots}=
 \frac{\partial^N}{\partial\E^{a\ldots}}\theta(\tilde{\mu}I-
 \F(\E))\bigg|_{\E=0} \label{DDeriv1}.
\end{equation}
 The Fockian may also be expanded order by order in the perturbation to yield
\begin{equation}\label{FockianTaylor}
  \begin{split}
    \F(\E)=\F^{0} & +\sum_a \F^{a}\E^{a}\\
    &+\frac{1}{2!}\sum_{ab} \F^{ab}\E^{a}\E^{b}\\
    &+\frac{1}{3!}\sum_{abc} \F^{abc}\E^{a}\E^{b}\E^{c}+\dots ~,
  \end{split}
\end{equation}
where $\F^{a}$ stands for $\partial\F(\E)/\partial\E^{a}$,
$\F^{ab}=\partial^2\F(\E)/\partial\E^{a}\partial\E^{b}$,
and so on for the higher order terms.
A similar expansion also holds for the density matrix $\D(\E)$.

Within HF theory, the unperturbed Fockian $F^0$ in the non-orthogonal basis is 
\begin{equation}
F^0=h^0+J(D^0)+K(D^0), \label{fockian0}
\end{equation}
while the first variation of the Fockian is 
\begin{equation}
F^a=\mu_a+J(D^a)+K(D^a)
\end{equation}
and the higher terms are given by 
\begin{equation}
F^{ab\ldots}=J(D^{ab\ldots})+K(D^{ab\ldots}). \label{fockianN}
\end{equation}
In computation of the unperturbed Fockian, the Coulomb matrix $J$ may be computed in ${\cal O}(N{\rm lg}N)$ 
with the Quantum Chemical Tree Code (QCTC) \cite{MChallacombe97} and the
exchange matrix $K$ computed in ${\cal O}(N)$ with the ${\cal O}(N)$-exchange  (ONX) algorithm 
that exploits quantum locality of the density matrix $D^0$ \cite{ESchwegler97}.
Likewise the Fockian derivatives, $F^{a\ldots}$, may be computed 
with the same algorithms in linear scaling time if elements of 
$D^{a\ldots}$ manifest an approximate exponential decay with atom-atom separation, 
similar to the decay properties of $D^0$. 

While the expansions above are given explicitly for Hartree-Fock Theory, similar expressions 
hold also for Kohn-Sham and hybrid HF/DFT,  which involve variation of the  exchange-correlation 
matrix $V_{xc}^{a\ldots}(D^0,D^a,\ldots)$ \cite{Lee_1994,PSalek02}.

\subsection{Conditions for self-consistency}\label{SelfConsistency}

The derivative density matrices and derivative Fockians depend on 
each other implicitly, and must be solved for self-consistently
via the CPSCF equations.
The necessary and sufficient criteria for convergence of the 
CPSCF equations involve generalized self-consistence conditions \cite{Furche_2001},
\begin{gather}
    [\F^{0} ,\D^{0}]=0,\label{eq:commutators1}\\
    [\F^{a} ,\D^{0}]+[\F^{0},\D^{a}]=0,\label{eq:commutators2}\\
  \begin{split}
    [\F^{ab},\D^{0}]&+\frac{1}{2}[\F^{a},\D^{b}]+\frac{1}{2}[\F^{b},\D^{a}] \\
    &+[\F^{0},\D^{ab}]=0,\label{eq:commutators3}\\
  \end{split}\\
  \begin{split}
    [\F^{abc},\D^{0}]&+\frac{1}{3}[\F^{ab},\D^{c}]+\frac{1}{3}[\F^{ac},\D^{b}]\\
    &+\frac{1}{3}[\F^{ab},\D^{c}]+\frac{1}{3}[\F^{a},\D^{bc}]\\
    &+\frac{1}{3}[\F^{b},\D^{ac}]+\frac{1}{3}[\F^{c},\D^{ab}]\\
    &+[\F^{0},\D^{abc}]=0,\label{eq:commutators4}\\
  \end{split}
\end{gather}
%\begin{subequations}\label{eq:commutators}
%  \begin{gather}
%    [\F^{0} ,\D^{0}]=0,\\
%    [\F^{a} ,\D^{0}]+[\F^{0},\D^{a}]=0,\\
%    [\F^{ab},\D^{0}]+\frac{1}{2}[\F^{a},\D^{b}]+\frac{1}{2}[\F^{b},\D^{a}]+[\F^{0},\D^{ab}]=0,\\
%    [\F^{abc},\D^{0}]
%+\frac{1}{3}[\F^{ab},\D^{c}]+\frac{1}{3}[\F^{ac},\D^{b}]+\frac{1}{3}[\F^{ab},\D^{c}]
%+\frac{1}{3}[\F^{a},\D^{bc}]+\frac{1}{3}[\F^{b},\D^{ac}]+\frac{1}{3}[\F^{c},\D^{ab}]
%+[\F^{0},\D^{abc}]=0,\\
%  \end{gather}
%\end{subequations}
in addition to generalized idempotency-like constraints \cite{Furche_2001},
\begin{gather}
  \D^{0} =\D^{0} \D^{0},\label{eq:anticommutators1}\\
  \D^{a} =\{\D^{a},\D^{0}\},\label{eq:anticommutators2}\\
  \D^{ab}=\{\D^{ab},\D^{0}\}+\frac{1}{2}\{\D^{a},\D^{b}\}\label{eq:anticommutators3}\\
  \begin{split}
    \D^{abc}=\{\D^{abc},\D_0\}&+\frac{1}{3}\{\D^{ab},\D^{c}\}+\frac{1}{3}\{\D^{ac},\D^{b}\}\\
    &+\frac{1}{3}\{\D^{bc},\D^{a}\}\label{eq:anticommutators4},
  \end{split}
\end{gather}
%\begin{subequations}\label{eq:anticommutators}
%  \begin{gather}
%    \D^{0} =\D^{0} \D^{0},\\
%    \D^{a} =\{\D^{a},\D^{0}\},\\
%    \D^{ab}=\{\D^{ab},\D^{0}\}+\frac{1}{2}\{\D^{a},\D^{b}\}\\
%    \D^{abc}=\{\D^{abc},\D_0\}+
%        \frac{1}{3}\left(\{\D^{ab},\D^{c}\}+\{\D^{ac},\D^{b}\}
%	                 +\{\D^{bc},\D^{a}\}\right),
%  \end{gather}
%\end{subequations}
where the anti-commutator notation $\{A,B\} = AB+BA$ has been used.

\section{Solving the higher order CPSCF equations with Perturbed Projection}

In solution of the CPSCF equations, it is first necessary to determine the ground state density matrix $\mathcal{D}^0$.  This may 
be accomplished in ${\cal O}(N)$ using a purification algorithm such as Niklasson's \cite{ANiklasson02A} 
second order trace correcting scheme (TC2) in conjunction with sparse atom-blocked linear algebra 
\cite{ANiklasson03,MChallacombe00B}.  Linear scaling is achieved for insulating systems through 
the dropping (filtering) of atom-atom blocks with Frobenious norm below a numerical threshold 
($\tau \sim 10^{-4}-10^{-7}$).  At SCF convergence the TC2 algorithm generates a polynomial sequence 
defining the ground state projector, from which the derivative density matrices are directly obtained.

Having solved the groundstate SCF equations, solution of the CPSCF equations commences with a guess at the 
derivative densities, followed by computation of derivative Fockians.  At the $r^{\rm th}$ 
CPSCF cycle, the $n^{\rm th}$ order derivative Fockians are 
\begin{equation}
    F^{a\ldots}_{r}= \left\{
    \begin{array}{ll}
      \mu_a+J(D^{a}_r)+K(D^{a}_r), & n=1\label{FockBuild}\\
      J(D^{a\ldots}_r)+K(D^{a\ldots}_r), & n>1 \,.
    \end{array}\right.
\end{equation}
After construction of the derivative Fockians, response functions through 
$D^{a\ldots}_{r+1}$ are computed, constituting one cycle in solution of the CPSCF.  
As described in Section \ref{ResponseFunctions},  these response functions are 
obtained directly through variation of the occupied subspace projector,  
\begin{equation}
%    \displaystyle\D^{a\ldots}_{r+1}=n!
    \displaystyle\D^{a\ldots}_{r+1}=
    \frac{\partial^n}{\partial\E^{a\ldots}}\theta(\tilde{\mu}I-
    \F_r(\E))\bigg|_{\E=0} \, , \label{DDeriv}
\end{equation}
which is accomplished via the Niklasson and Challacombe density matrix 
perturbation theory \cite{ANiklasson04}.  

After a few CPSCF cycles, the approach to self-consistency may be accelerated with 
Weber and Daul's DDIIS algorithm \cite{VWeber03}, 
\begin{equation}
    \displaystyle\widetilde{\F}^{a\ldots}_{r}=\sum_{k=r-s}^{r}c_k \F^{a\ldots}_{k} \label{DDIISEq} \\
\end{equation}
in which the $c_k$ coefficients are chosen to minimize the 
$n^{\rm th}$ order commutation 
relations, as in Eqs.~(\ref{eq:commutators2})-(\ref{eq:commutators4}). The application of the
DDIIS algorithm to acceleration of higher order CPSCF equations is developed further in Section \ref{DDIIS}.

At self-consistency, the conditions given in Section \ref{SelfConsistency} have been met, and it is then 
appropriate to compute response properties.   In general, we can use the expectation value 
\begin{equation}
E^{(\gamma)} = 2 (\gamma-1)! Tr(D^{(\gamma-1)} h^{(1)}), \label{Np1Rule}
\end{equation}
where in the case of polarization we have the (hyper)polarizabilities $\alpha_{ab}=-2\Tr[D^a\mu_b]$, 
$\beta_{abc}=-4\Tr[D^{ab}\mu_c]$ or $\gamma_{abcd}=-12\Tr[D^{abc}\mu_d]$.
Alternatively, it is possible to construct density matrix analogues of 
Wigner's $2 n+1$ rule, which allows the evaluation of response properties up to order $2 n+1$ from response 
functions of order $n$.  These analogues are given in Section \ref{Wigner2Np1} through $4^{th}$ order in
the energy.

\subsection{Perturbed Projection}\label{ResponseFunctions}

Although a number of analytic, asymptotically discontinuous representations exist for the Heaviside 
step function $\theta$, direct representation (and variation) of these forms as in Eq.~(\ref{DDeriv}) 
is problematic. Polynomial expansion of the step function is a alternative choice, but demands a very high 
order and can be costly.  Specifically, polynomial expansion of $\theta$ with a $p$'th order polynomial 
incurs a cost that is at best ${\cal O}(\sqrt{p})$ \cite{WLiang04}.  Polynomial expansion techniques,
such as those based on the the Tchebychev polynomials, may also be plagued by Gibbs oscillations \cite{AVoter96},
which are high order ripples in the approximate $\theta$ due to incompleteness.
Alternatively, recursive purification methods achieve high order representation in 
${\cal O}( \log p )$ \cite{ANiklasson03}.  Also, purification methods (such as TC2 and TRS4) 
yield projectors that are smooth and strictly monotonic,  

Each Perturbed Projection sequence is based on a corresponding  purification scheme or generator, such as TC2 \cite{ANiklasson02A}.   
The Perturbed Projection sequence is obtained by collecting terms of the response order by 
order upon perturbative expansion of its generator \cite{ANiklasson04}.
Perturbed Projection provides explicit, recursive formulae 
for the construction of response functions, retaining the convergence properties,  smoothness and 
montonicity of the generating sequence.   These explicit formulae stand in contrast to methods where the 
density matrix derivatives are implicitly defined as solutions to equations of Sylvester type \cite{Ochsenfeld97,HLarsen01a,COchsenfeld04}.

Sufficient to compute fourth order properties using the $2 n+1$ rule presented in Section \ref{Wigner2Np1}, 
Perturbed Projection is outlined in the following for computation of the second order response function.  
The Perturbed Projection sequence is started with the  $\X^{a\ldots}_{0}$, which are 
prepared from the Fockians $\F^0$, $\F^a$, and $\F^{ab}$ by  compressing their spectrum into the domain of 
convergence \cite{ANiklasson02A} using
\begin{equation}
    \X^0_{0}=\frac{\F_{max}-\F^0}{\F_{max}-\F_{min}} 
\end{equation}
and 
\begin{equation}
    \X^{a\ldots}_{0}=\frac{\F^{a\ldots}_{n}}{\F_{min}-\F_{max}},
\end{equation}
where $\F_{min}$ and $\F_{max}$ are upper and lower bounds to the eigenvalues of $\F^0$.  

While Perturbed Projection can be formulated within any purification scheme, we focus here on the
simple and efficient TC2 method \cite{ANiklasson02A}.  Briefly, TC2 
constructs a ground state projector through a series of trace correcting projections;  
when the trace is larger than $N_{\rm e}$, $x^2$ is used to reduce the trace, and 
when the trace is less than  $N_{\rm e}$, $2 x-x^2$ is used to increase the trace.  
The resulting sequence of correcting projections yields a step at the correct chemical potential. 
Within this framework, the second order TC2 Perturbed Projection sequence is 
\begin{equation}\label{PP1}
\left.
\begin{array}{ll}
\X^{ab}_{i+1}&=\{\X^{ab}_{i},\X^0_{i}\}+\frac{1}{2}\{\X^a_{i},\X^b_{i}\}\\
\X^a_{i+1}&=\{\X^a_{i},\X^0_{i}\}\\
\X^b_{i+1}&=\{\X^b_{i},\X^0_{i}\}\\
\X^0_{i+1}&=(\X^0_{i})^2 
\end{array} 
\right\}\,  {\rm Tr}[\mathcal{X}^0_{i}]\ge N_e 
\end{equation}
or 
\begin{equation}\label{PP2}
\left.
\begin{array}{ll}
      \X^{ab}_{i+1}&=2 \X^{ab}_{i}-(\{\X^{ab}_{i},\X^0_{i}\}+\frac{1}{2}\{\X^a_{i},\X^b_{i}\})\\
      \X^b_{i+1}&=2 \X^b_{i}-\{\X^b_{i},\X^0_{i}\} \\
      \X^a_{i+1}&=2 \X^a_{i}-\{\X^a_{i},\X^0_{i}\} \\
      \X^0_{i+1}&=2 \X^0_{i}-(\X^0_{i})^2
\end{array} 
\right\}\, {\rm Tr}[\mathcal{X}^0_{i}]< N_e.
\end{equation}
As with the  TC2 generator, the $n^{\rm th}$ order response functions
\begin{equation}
 \D^{a...} = n!\lim_{i\rightarrow\infty} \X_i^{a...},
\end{equation}
converge quadratically, reaching  convergence when either 
the error $\varepsilon = |\Tr[\X^0_{i}]-N_e| + |\Tr[\X^{a}_{i}]| + \ldots$, or 
the maximum element in the change $\delta \X^{a\ldots} = |\X^{a\ldots}_{i+1}-\X^{a\ldots}_{i}|$ 
falls below $\tau$, the atom-atom block drop tolerance described in Ref.~[\onlinecite{ANiklasson03}].
As described more completely in Ref.~[\onlinecite{ANiklasson05a}], when the solution gets 
close to convergence, i.e.~$|\Tr[\X^0_{i}]-N_e|< \epsilon$ with 
$\epsilon \approx 10^{-1}-10^{-3}$, we alternate the projection at each step,
which protects the convergence under the incomplete sparse linear algebra.

\subsection{Derivative DIIS}\label{DDIIS}

Direct inversion in the iterative subspace (DIIS), introduced
some time ago by Pulay \cite{Pulay80,Pulay82}, accelerates convergence toward 
self-consistency.   DIIS employs information accumulated during preceding 
iterations to construct an effective Fockian $\widetilde \F_{k}$ 
at the $k$-th SCF cycle, which minimizes the commutation error between the Fockian
and the density matrix. The effective Fockian is then used instead of $\F_{k}$
to generate an improved density matrix.  

Recently, Weber and Daul have developed the Derivative DIIS (DDIIS) scheme for accelerating 
convergence of the CPSCF equations \cite{VWeber03}.  Like DIIS, DDIIS is based on 
minimization of the Frobenious norm of an error matrix 
\begin{equation}
  \widetilde e_r^{a\ldots}=\sum_{i=r-s}^{r}c_i e_i^{a\ldots},
\end{equation}
where the $e_i^{a\ldots}$'s are just the $n$-th order commutator relation
of Eqs.~(\ref{eq:commutators1}-\ref{eq:commutators4}) (e.g. the first order error matrix 
is given by $e_i^{a}=[\F^{a}_i ,\D^{0}]+[\F^{0},\D^{a}_i]$).
The optimal coefficients $c_i$ are solutions to the quadratic programming problem
\begin{equation}
  \inf \left \{-\frac{1}{2}\sum_{i,j=r-s}^nc_iB_{ij}c_j,\, \sum_{i=r-s}^r c_i=1 \right \},
\end{equation}
where elements of the $\mathbf{B}$ matrix are given by 
$B_{ij}=\Tr[e_i^{a\ldots}(e_j^{a\ldots})^T]$.
A working equation is then obtained through the associated Euler-Lagrange equation
\begin{equation}\label{eq:diismatrix}
 \left ( \begin{array}{cc}
     \mathbf{B}     & \mathbf{1} \\
     \mathbf{1}^{T} & 0 
   \end{array}\right )
 \cdot \left ( \begin{array}{c}
     \mathbf{c}     \\
     \lambda  
   \end{array}\right )
  =  \left ( \begin{array}{c}
     \mathbf{0}      \\
         1  
   \end{array}\right ) \, ,
\end{equation}
 where $\mathbf{0}=(0,\ldots,0)^{T}$ and $\mathbf{1}=(1,\ldots,1)^{T}$ are
 vectors whose components are 0 and 1 respectively 
 and $\lambda$ is the Lagrange multiplier of the constraint 
 $\sum_{i=n-m}^{n}c_{i}=1$. The set of linear equations is solved
 by inverting the left-hand side matrix.  In the event of a singular or near singular
 matrix, the rank of Eq.~(\ref{eq:diismatrix}) is reduced by discarding the oldest entries (increasing $s$)
 until the linear system stabilizes.

\subsection{Density matrix formulation of Wigner's $2n+1$ rule }\label{Wigner2Np1}

Wigner's $2n+1$ rule, traditionally predicated on derivatives of the wavefunction,
yields order $2n+1$ in the energy response from $n$-th order derivatives \cite{SKarna91,SEpstein74}. 
A density matrix analogue of Wigner's $2n+1$ rule for a single perturbation parameter
was given to third order by McWeeny \cite{RMcWeeny62} and up to fourth order by 
Niklasson and Challacombe \cite{ANiklasson04} in the orthogonal representation.  
For completeness, we present non-orthogonal generalizations up to fourth order, 
which may require less memory due to the less dense structure of non-orthogonal matrix intermediates.

The first and second order energy corrections are well known, corresponding simply to expectation
values as in Eq.~(\ref{Np1Rule}).  Beyond second order, the $2 n+1$ rule offers a valuable alternative. 
The third order non-orthogonal contribution is 
\begin{align}
%  \begin{split}
    E^{abc}
%&=\frac{1}{6}P(a,b,c)\Tr([\D^a,\D^0]\D^b\F^c)\label{eq:2n+1 third order ortho}\\
%    &=\frac{1}{6}\sum_{P(a,b,c)}\Tr[[D^a,D^0]_SSD^bF^c]\label{eq:2n+1 third order}
    &=2\sum_{P(a,b,c)}\Tr[[D^a,D^0]_SSD^bF^c]\label{eq:2n+1 third order}
%  \end{split}
\end{align}
where $P(a,b,c)$ stands for the permutation operator such that all
permutations of $a$, $b$ and $c$ are made (e.g. $P(a,b,c)$ generates the sum of
all the six terms: $(a,b,c)$, $(a,c,b)$, $(b,a,c)$, $(b,c,a)$, $(c,a,b)$ and $(c,b,a)$)
and $[A,B]_S=ASB-BSA$ where $S$ is the overlap matrix.  Similarly,  the fourth order non-orthogonal 
contribution is 
\begin{align}
  \begin{split}
%    E^{abcd} =&\frac{1}{24}\sum_{P(a,b,c,d)}\Tr[[D^{ab},D^0]_S SD^cF^d + \\ 
    E^{abcd} =&\sum_{P(a,b,c,d)}\Tr[[D^{ab},D^0]_S SD^cF^d + \\ 
    &[D^a,D^0]_S S(D^{bc}F^d +D^bF^{cd})]\label{eq:2n+1 fourth order}
  \end{split}
\end{align}

For the orthogonal case $S=I$ and $D^a,F^c,\ldots$ are replaced by $\D^a,\F^c,\ldots$.
In most cases the complexity of these equations
can be reduced by taking advantage of indicial symmetry; 
$a,b,c,$ and $d$ represent the Cartesian directions $x,y,z$
so that terms with indecies in the same direction simplify. 
For example, $E^{aaaa}$ reduces to only one term requiring only 15 matrix multiplications.
In the worst case, where all the directions are different, i.e. $E^{aabc}$ 
(or any other permutation of $(a,a,b,c)$), the relation \eqref{eq:2n+1 fourth order} 
reduces to include only 12 terms with 180 matrix multiplications. 
Similar reductions of the computational cost also apply to 
Eq. \eqref{eq:2n+1 third order}. The number of matrix-matrix multiplies can be further reduced
if one uses an orthogonal representation, but this typically involves matrix-matrix multiplies 
with more dense intermediates.


\commentoutA{

\begin{table}[h]
  \centering
  \caption{\protect
    The longitudinal polarizability, $\alpha_{zz}$, for water chains at the RHF/6-31G level
    of theory, computed with {\sc MondoSCF} using {\tt GOOD} and {\tt TIGHT} numerical thresholds, 
   and also with the {\sc GAMESS} quantum chemistry package \cite{gamess}.
  }\label{tab:Alpha_1D_Values}
  \begin{tabular}{cccc}
    \toprule
    $N_{\rm H_2 O}$ &\multicolumn{1}{c}{{\sc GAMESS}}
        &\multicolumn{1}{c}{{\tt GOOD}}
        &\multicolumn{1}{c}{{\tt TIGHT}}\\
    \hline
    1 & 5.8136 & 5.813620 & 5.813588 \\
    2 & 6.3448 & 6.345037 & 6.344822 \\
    3 & 6.5844 & 6.584658 & 6.584435 \\
    4 & 6.7276 & 6.727905 & 6.727672 \\
    5 & 6.8226 & 6.823290 & 6.822857 \\
   10 & 7.0308 & 7.031056 & 7.030858 \\
   15 & 7.1047 & 7.104904 & 7.104770 \\
   20 & 7.1424 & 7.142580 & 7.142422 \\
%#ALPHA/N
%#    LOOSE    GOOD     TIGHT    VERYTIGHT      
% 1 5.817223 5.813620 5.813588 5.813593
% 2 6.350891 6.345037 6.344822 6.344813
% 3 6.588345 6.584658 6.584435 6.584420
% 4 6.738352 6.727905 6.727672 6.727651
% 5 6.836304 6.823290 6.822857 6.822602
%10 7.031962 7.031056 7.030858 7.030828
%15 7.104774 7.104904 7.104770 7.104738
%20 7.141731 7.142580 7.142422 7.142388
    \botrule
  \end{tabular}
\end{table}

\begin{table}[h]
  \centering
  \caption{\protect
    The longitudinal first hyperpolarizability, $\beta_{zzz}$,
    for water chains at the RHF/6-31G level of theory, computed with 
    {\sc MondoSCF} using {\tt GOOD} and {\tt TIGHT} numerical thresholds, 
    and also with the {\sc GAMESS} quantum chemistry package \cite{gamess}.
  }\label{tab:Beta_1D_Values}
  \begin{tabular}{cccccc}
    \toprule
    $N_{\rm H_2 O}$ &\multicolumn{1}{c}{{\sc GAMESS}}
    &\multicolumn{1}{c}{{\tt GOOD}}
    &\multicolumn{1}{c}{{\tt GOOD}\footnote[1]{The density matrix based 2n+1 rule has been used.}}
    &\multicolumn{1}{c}{{\tt TIGHT}}
    &\multicolumn{1}{c}{{\tt TIGHT}$^a$} \\
    \hline
     1 & -30.6125 & -30.611029 & -30.612627 & -30.612163 & -30.612256 \\
     2 & -29.5444 & -29.547427 & -29.548504 & -29.544907 & -29.544994 \\
     3 & -25.3696 & -25.372208 & -25.373615 & -25.370297 & -25.370381 \\
     4 & -22.1411 & -22.143436 & -22.145040 & -22.141494 & -22.141603 \\
     5 & -19.8925 & -19.902088 & -19.904449 & -19.896462 & -19.897141 \\
    10 & -14.8063 & -14.807075 & -29.617990 & -14.806973 & -14.807119 \\
    15 & -12.9713 & -12.969238 & -12.972227 & -12.971940 & -12.972124 \\
    20 & -12.0334 & -12.028709 & -12.033633 & -12.034014 & -12.034238 \\
%#BETA/N
%#        LOOSE                GOOD               TIGHT              VERYTIGHT      
%#     SCF       2n+1      SCF       2n+1      SCF       2n+1      SCF       2n+1
% 1 30.465445 30.541530 30.611029 30.612627 30.612163 30.612256 30.612280 30.612298
% 2 29.457319 29.533744 29.547427 29.548504 29.544907 29.544994 29.544586 29.544608
% 3 25.346540 25.422709 25.372208 25.373615 25.370297 25.370381 25.369861 25.369874
% 4 22.188342 22.251189 22.143436 22.145040 22.141494 22.141603 22.140943 22.140951
% 5 20.029374 20.053896 19.902088 19.904449 19.896462 19.897141 19.892449 19.892476
%10 15.004856 15.043761 14.807075 29.617990 14.806973 14.807119 14.806291 14.806326
%15 13.099761 13.226662 12.969238 12.972227 12.971940 12.972124 12.971260 12.971298
%20 12.115774 12.299637 12.028709 12.033633 12.034014 12.034238 12.033362 12.033402
    \botrule
  \end{tabular}
\end{table}

\begin{table}[h]
  \centering
  \caption{\protect
    The longitudinal second hyperpolarizability, $\gamma_{zzzz}$,
    for water chains at the RHF/6-31G level of theory, computed with 
    {\sc MondoSCF} using {\tt GOOD} and {\tt TIGHT} numerical thresholds, 
    and also with the {\sc GAMESS} quantum chemistry package \cite{gamess}.
  }\label{tab:Gamma_1D_Values}
  \begin{tabular}{rccccc}
    \toprule
    $N_{\rm H_2 O}$ &\multicolumn{1}{c}{{\sc GAMESS}}
    &\multicolumn{1}{c}{{\tt GOOD}}
    &\multicolumn{1}{c}{{\tt GOOD}\footnote[1]{The density matrix based 2n+1 rule has been used.}}
    &\multicolumn{1}{c}{{\tt TIGHT}}
    &\multicolumn{1}{c}{{\tt TIGHT}$^a$} \\
    \hline
     1 &  330.5753 & 330.54375 & 330.54319 & 330.57193 &  330.5724 \\
     2 &  820.1398 & 820.19231 & 820.19667 & 820.14775 &  820.1493 \\
     3 & 1008.5656 & 1008.5855 & 1008.6073 & 1008.5752 & 1008.5765 \\
     4 & 1103.4813 & 1103.5053 & 1103.5280 & 1103.4883 & 1103.4902 \\
     5 & 1168.9563 & 1169.2104 & 1169.2531 & 1169.0630 & 1169.0754 \\
    10 & 1324.2906 & 1325.1208 & 1324.6321 & 1324.2975 & 1324.2999 \\
    15 & 1381.8657 & 1383.0802 & 1382.2322 & 1381.8758 & 1381.8738 \\
    20 & 1411.4264 & 1414.4092 & 1411.8528 & 1411.4410 & 1411.4292 \\
%#GAMMA/N
%#            LOOSE                   GOOD                    TIGHT                  VERYTIGHT      
%#       SCF         2n+1        SCF         2n+1        SCF         2n+1        SCF         2n+1
% 1  326.453668  327.657510  330.543748  330.543185  330.571927  330.572391  330.575212  330.575172
% 2  811.011253  813.149259  820.192314  820.196669  820.147747  820.149334  820.143928  820.143860
% 3  997.747780 1000.161063 1008.585466 1008.607306 1008.575157 1008.576475 1008.567843 1008.567803
% 4 1092.976882 1095.290284 1103.505309 1103.528026 1103.488253 1103.490194 1103.478891 1103.478823
% 5 1168.090056 1168.020522 1169.210364 1169.253062 1169.063044 1169.075387 1168.954564 1168.954518
%10 1310.823281 1311.068128 1325.120848 1324.632059 1324.297519 1324.299857 1324.289143 1324.289082
%15 1375.184287 1364.961009 1383.080212 1382.232175 1381.875768 1381.873821 1381.864844 1381.864664
%20 1418.062288 1388.371139 1414.409224 1411.852777 1411.440956 1411.429186 1411.425500 1411.425048
    \botrule
  \end{tabular}
\end{table}



\begin{figure}[h]
  \caption{\protect
    Total CPU time of the fifth CPSCF iteration of fourth order for
    the water cluster sequence with the 6-31G and 6-31G** 
    basis sets and the {\tt GOOD} and {\tt TIGHT} 
    numerical thresholds (see text) controlling numerical
    precision of the result. The lines are fits to the 
    last three and four points, respectively.
  }\label{fig:Gamma_scaling}
  \resizebox*{3.6in}{!}{\includegraphics[angle=-90.00]{Gamma_h2o3D_6-31G_6-31Gss_G_T}}
\end{figure}

\begin{figure}[h]
  \caption{\protect
    QCTC CPU time of the fifth CPSCF iteration of fourth order for
    the water cluster sequence with the 6-31G and 6-31G** 
    basis sets and the {\tt GOOD} and {\tt TIGHT} 
    numerical thresholds (see text) controlling numerical
    precision of the result. The lines are fits to the 
    last three and four points, respectively.
  }\label{fig:Gamma_QCTC_Timing}
  \resizebox*{3.6in}{!}{\includegraphics[angle=-90.00]{Gamma_QCTC_T}}
\end{figure}

\begin{figure}[h]
  \caption{\protect
    TC2 CPU time of the fifth CPSCF iteration of fourth order for
    the water cluster sequence with the 6-31G and 6-31G** 
    basis sets and the {\tt GOOD} and {\tt TIGHT} 
    numerical thresholds (see text) controlling numerical
    precision of the result. The lines are fits to the 
    last three and four points, respectively.
  }\label{fig:Gamma_TC2R_Timing}
  \resizebox*{3.6in}{!}{\includegraphics[angle=-90.00]{Gamma_TC2R_T}}
\end{figure}

\begin{figure}[h]
  \caption{\protect
    ONX CPU time of the fifth CPSCF iteration of fourth order for
    the water cluster sequence with the 6-31G and 6-31G** 
    basis sets and the {\tt GOOD} and {\tt TIGHT} 
    numerical thresholds (see text) controlling numerical
    precision of the result. The lines are fits to the 
    last three and four points, respectively.
  }\label{fig:Gamma_ONX_Timing}
  \resizebox*{3.6in}{!}{\includegraphics[angle=-90.00]{Gamma_ONX_T}}
\end{figure}

\begin{figure}[h]
  \caption{Total CPU times with increasing order of the response for 
           the fifth CPSCF cycle computed as the $n+1$ expectation value, Eq.~(\ref{Np1Rule}).}\label{TimeWithOrder}
  \resizebox*{3.6in}{!}{\includegraphics[angle=-90.00]{Mix_h2o3D_6-31G_T.ps}}
\end{figure}

\begin{figure}[t]
  \caption{\protect
    Superposition of the magnitudes of the RHF/6-31G density matrix
    derivative elements $D_{cd}$, $D^{x}_{cd}$, $D^{xx}_{cd}$ and $D^{xxx}_{cd}$
    along the $x$ axis with the separation of basis function centers
    for $\rm (H_2O)_{150}$. The density matrix 
    derivatives have been converged to within {\tt TIGHT} (e.g. 
    a matrix threshold $\tau=10^{-6}$ $[a.u.]$).
  }\label{fig:Superposition_Decay}
  \resizebox*{3.6in}{!}{\includegraphics[angle=-90.00]{DMix_40pts}}
\end{figure}

}

\section{Results} \label{RESULTS}

We have implemented these methods in the MondoSCF suite of linear scaling quantum chemistry 
programs \cite{MondoSCF}.  The construction of the Fockian and derivative Fockian was carried
out using the linear scaling QCTC method for computation of the Coulomb matrix \cite{MChallacombe97,CTymczak04a} and the
ONX algorithm \cite{ESchwegler97,CTymczak04b} for computation of the Hartree-Fock exchange matrix. 
The CPSCF equations were solved in an entirely orthogonal representation,
with higher order properties evaluated using the $2 n+1$ rule with the non-orthogonal formulae given
by Eqs.~(\ref{eq:2n+1 third order}-\ref{eq:2n+1 fourth order}).  Two different levels of numerical 
accuracy have been used, {\tt GOOD} and {\tt TIGHT}.   Thresholds that define the {\tt GOOD} 
accuracy level include a matrix threshold $\tau=10^{-5}$, as well as other numerical thresholds 
detailed in Ref.~[\onlinecite{CTymczak04a}], which deliver 6 digits of relative accuracy in the total 
energy.  The {\tt TIGHT} option involves the matrix threshold $\tau=10^{-6}$ and delivers 8 
digits of relative accuracy in the total energy.  

Calculations were carried out on a single Intel Xeon 2.4GHz processor running RedHat Linux 8.0 and  
executables compiled with Portland Group Fortran Compiler pgf90 4.0-2 \cite{PGF90}.

Convergence of the CPSCF equations for the water systems described in the following 
are typically achieved in about 10 cycles, independent of cluster size, basis set,
matrix threshold or order of the response calculated.

All results are reported in atomic units.  Also, unless otherwise noted, all timings and values 
have been obtained by computing the $n^{\rm th}$ order response function and evaluation with the
$n+1$ rule (expectation value), Eq.~(\ref{Np1Rule}).
\subsection{One dimensional water chains}

Perturbed Projection has been used to compute the (hyper)polarizabilities $\alpha_{zz}$, 
$\beta_{zzz}$ and $\gamma_{zzzz}$ of linear water chains up to (H$_2$O)$_{20}$. 
These calculations have been carried out with {\sc MondoSCF} at the RHF/6-31G level of theory using 
both the {\tt GOOD} and {\tt TIGHT} thresholding parameters,  as well as with the conventional algorithms implemented
in the {\sc GAMESS} quantum chemistry package \cite{gamess}.  These static properties have 
been evaluated at the geometries given by Otto {\em et al.}~\cite{POtto99}, and the GAMESS
results are given to the number of digits provided by that program.
The {\sc MondoSCF} results have been obtained both as expectation values, given by Eq.~(\ref{Np1Rule}), 
and using the non-orthogonal density matrix $2 n+1$ rules given in 
Eqs.~(\ref{eq:2n+1 third order}-\ref{eq:2n+1 fourth order}).

As a benchmark, we have also carried out calculations with the linear chain (H$_2$O)$_{20}$ with 
the {\tt VERYTIGHT} numerical thresholding parameters, which employ a $10^{-7}$ drop tolerance and
aim to provide 10 digits of precision in the total energy.   These {\tt VERYTIGHT} calculations yield  
$\alpha_{zz}=7.142422\,a.u.$, $\beta_{zzz}=-12.033362\,a.u.$ and $\gamma_{zzzz}=1411.425500\,a.u.$. 

\subsection{Linear scaling: 3D water clusters}

Linear scaling computation of the RHF/6-31G and RHF/6-31G** second hyperpolarizability,
achieved with Perturbed Projection, is shown for three-dimensional water clusters 
in Fig~\ref{fig:Gamma_scaling}.  These timings are the total CPU time for the fifth CPSCF cycle, 
including build time for $\F^{abc}$ ({\sc ONX} and {\sc QCTC}), 
iterative construction of $\D^{abc}$ (Perturbed Projection via {\sc TC2}) and all intermediate
steps including the congruence transformation.
A breakdown of the dominant contributions to these totals are given
in Figs~\ref{fig:Gamma_QCTC_Timing}-\ref{fig:Gamma_ONX_Timing}, which show timings for Coulomb 
summation ({\sc QCTC}), Perturbed Projection ({\sc TC2}), and exact exchange ({\sc ONX}).

Figure \ref{TimeWithOrder} shows the increase in cost associated with computing  higher order response functions.
Corresponding to this increase in cost, Fig.~\ref{fig:Superposition_Decay} shows the magnitude of atom-atom blocks 
of density matrix derivatives up to third order as a function of atom-atom distance 
when perturbed by a static electric field.  The density response shows an approximate 
exponential decay as a function of internuclear distance with the rate of  decay
slowing slightly and the distribution shifted up with increasing order in the perturbation.

\section{Discussion}

In our current formulation, the increase in magnitude and reduction of locality in elements of the 
response function makes achieving linear scaling more difficult with increasing order in perturbation. 
Nevertheless, linear scaling has been achieved at the HF level of theory up to fourth order (i.e.~$\gamma$) 
in the total energy for three-dimensional systems and  non-trivial basis sets.  At fourth order,
Perturbed Projection and exact exchange were the dominant costs in solving the CPSCF equations,
as shown in Figs~\ref{fig:Gamma_TC2R_Timing} and \ref{fig:Gamma_ONX_Timing}.  For the fourth order Perturbed Projection, 
$N$-scaling is achieved 
between 70 to 110 water molecules, depending on $\tau$.  Despite a nearly dense $D^{abc}$, the 
dominant work in its construction always involves multiplication with matrices that are significantly 
more sparse, as $\X^{abc}_i\X^{0}_i$ or $\X^{ab}_i\X^{c}_i$.   Likewise, $N$-scaling is achieved between 
70 to 90 water molecules for construction of the Hartree-Fock exchange contribution.  In this case, 
the approximate decay of the density matrix still leads to linear scaling through ordered skip out lists,
as described in Ref.~[\onlinecite{ESchwegler97}].  In both cases, the increase in response function magnitude
is equivalent to tightening numerical thresholds, which increases the cost and delays the onset of linear
scaling. 

In tables \ref{tab:Alpha_1D_Values}-\ref{tab:Gamma_1D_Values} we find that 
a reduction of the drop tolerance by one order of magnitude leads to an increase in precision by
1-2 significant digits, with {\tt GOOD} and {\tt TIGHT} yielding approximatively 3-4 and 5-6 
correct digits  independent of the response order.  We further observe about 1 extra digit of accuracy 
when using the $2n+1$ rule.  This might be expected from the higher order error propagation resulting from 
products of lower order response functions, relative to evaluation with  Eq.~(\ref{Np1Rule}), 
which involves an error that is always linear in a higher order derivative density matrix.  

As shown in Fig.~\ref{TimeWithOrder}, computing the second order response is significantly cheaper 
than the third order response, and involves an earlier onset of linear scaling. Because evaluation 
of properties with the $2 n+1$ rule is of negligible cost relative to solving the CPSCF equations, the cost 
difference for evaluating $\gamma$ with the $2 n+1$ rule relative to the $n+1$ expectation is just the 
difference (roughly 2:3) between the computation of $\beta$ and $\gamma$ shown in Fig.~\ref{TimeWithOrder}.  

\section{Conclusions}

Linear scaling has been demonstrated for the computation of response properties beyond second
order in the total energy using Perturbed Projection for solution of the Coupled-Perturbed Self-Consistent-Field 
equations.
In addition, we have provided details of 
the computational method, used three-dimensional systems and non-trivial basis sets
to demonstrate linear scaling and provided a (preliminary) assessment of 
error control.  Perturbed Projection for the 
computation of higher order response functions is quadratically convergent, simple 
to implement through higher orders and numerically stable.  Perturbed Projection is 
not unique to the Hartree-Fock model, the TC2 generator or the {\sc MondoSCF} 
$N$-scaling algorithms, but can be straightforwardly extended to models that include 
exchange correlation (DFT), other other purification schemes such as TRS4 \cite{ANiklasson03}
as well as other electronic structure programs.  

We have shown that response functions (density matrix derivatives) through fourth 
order are local upon {\em global} electric perturbation, corresponding to an 
approximate exponential decay of matrix elements.  However, the magnitude of the 
corresponding response functions increases with increasing perturbation order, equivalent 
to tightening of the matrix drop tolerance, $\tau$.   While we have not attempted to 
work out a detailed analysis for the propagation of error, it may be possible 
to develop a more effective thresholding scheme for high orders. In addition 
to being somewhat more accurate, the $2 n+1$ rule also provides a significantly 
cheaper alternative to the computation of expectation values and an earlier onset 
of linear scaling. 

A similar exponential decay in the first order response corresponding to a {\em local} nuclear 
displacement has likewise been demonstrated by Ochsenfeld and Head-Gordon \cite{Ochsenfeld97}. 
This behavior is expected to hold generally for both local and global perturbations to 
insulating systems.  Thus, the potential exists for Perturbed Projection to achieve linear scaling
for a large class of static molecular properties within the HF, DFT and hybrid HF/DFT model 
chemistries. Of particular interest, the recently developed non-orthogonal density matrix perturbation 
theory put forward in a proceeding article \cite{ANiklasson05a} may enable linear scaling computation of analytic 
second derivatives, which are important in computation of the Hessian matrix.  

\begin{acknowledgments}
 This work has been supported by the US Department of Energy 
 under contract W-7405-ENG-36 and the ASCI project.  
 The Advanced Computing Laboratory of Los 
 Alamos National Laboratory is acknowledged.
 All the numerical computations have been
 performed on computing resources located at this facility.
\end{acknowledgments}

\bibliography{Response3}
%\bibliography{mondo_new}

\commentoutB{

\clearpage

\begin{center}
\bf  TABLES\\[1.cm]
\end{center}

\begin{table}[h]
  \centering
  \caption{\protect
    The longitudinal polarizability, $\alpha_{zz}$, for water chains at the RHF/6-31G level
    of theory, computed with {\sc MondoSCF} using {\tt GOOD} and {\tt TIGHT} numerical thresholds, 
   and also with the {\sc GAMESS} quantum chemistry package \cite{gamess}.
  }\label{tab:Alpha_1D_Values}
  \begin{tabular}{cccc}
    \toprule
    $N_{\rm H_2 O}$ &\multicolumn{1}{c}{{\sc GAMESS}}
        &\multicolumn{1}{c}{{\tt GOOD}}
        &\multicolumn{1}{c}{{\tt TIGHT}}\\
    \hline
    1 & 5.8136 & 5.813620 & 5.813588 \\
    2 & 6.3448 & 6.345037 & 6.344822 \\
    3 & 6.5844 & 6.584658 & 6.584435 \\
    4 & 6.7276 & 6.727905 & 6.727672 \\
    5 & 6.8226 & 6.823290 & 6.822857 \\
   10 & 7.0308 & 7.031056 & 7.030858 \\
   15 & 7.1047 & 7.104904 & 7.104770 \\
   20 & 7.1424 & 7.142580 & 7.142422 \\
    \botrule
  \end{tabular}
\end{table}

\clearpage

\begin{table}[h]
  \centering
  \caption{\protect
    The longitudinal first hyperpolarizability, $\beta_{zzz}$,
    for water chains at the RHF/6-31G level of theory, computed with 
    {\sc MondoSCF} using {\tt GOOD} and {\tt TIGHT} numerical thresholds, 
    and also with the {\sc GAMESS} quantum chemistry package \cite{gamess}.
  }\label{tab:Beta_1D_Values}
  \begin{tabular}{cccccc}
    \toprule
    $N_{\rm H_2 O}$ &\multicolumn{1}{c}{{\sc GAMESS}}
    &\multicolumn{1}{c}{{\tt GOOD}}
    &\multicolumn{1}{c}{{\tt GOOD}\footnote[1]{The density matrix based 2n+1 rule has been used.}}
    &\multicolumn{1}{c}{{\tt TIGHT}}
    &\multicolumn{1}{c}{{\tt TIGHT}$^a$} \\
    \hline
     1 & -30.6125 & -30.611029 & -30.612627 & -30.612163 & -30.612256 \\
     2 & -29.5444 & -29.547427 & -29.548504 & -29.544907 & -29.544994 \\
     3 & -25.3696 & -25.372208 & -25.373615 & -25.370297 & -25.370381 \\
     4 & -22.1411 & -22.143436 & -22.145040 & -22.141494 & -22.141603 \\
     5 & -19.8925 & -19.902088 & -19.904449 & -19.896462 & -19.897141 \\
    10 & -14.8063 & -14.807075 & -29.617990 & -14.806973 & -14.807119 \\
    15 & -12.9713 & -12.969238 & -12.972227 & -12.971940 & -12.972124 \\
    20 & -12.0334 & -12.028709 & -12.033633 & -12.034014 & -12.034238 \\
    \botrule
  \end{tabular}
\end{table}

\clearpage

\begin{table}[h]
  \centering
  \caption{\protect
    The longitudinal second hyperpolarizability, $\gamma_{zzzz}$,
    for water chains at the RHF/6-31G level of theory, computed with 
    {\sc MondoSCF} using {\tt GOOD} and {\tt TIGHT} numerical thresholds, 
    and also with the {\sc GAMESS} quantum chemistry package \cite{gamess}.
  }\label{tab:Gamma_1D_Values}
  \begin{tabular}{rccccc}
    \toprule
    $N_{\rm H_2 O}$ &\multicolumn{1}{c}{{\sc GAMESS}}
    &\multicolumn{1}{c}{{\tt GOOD}}
    &\multicolumn{1}{c}{{\tt GOOD}\footnote[1]{The density matrix based 2n+1 rule has been used.}}
    &\multicolumn{1}{c}{{\tt TIGHT}}
    &\multicolumn{1}{c}{{\tt TIGHT}$^a$} \\
    \hline
     1 &  330.5753 & 330.54375 & 330.54319 & 330.57193 &  330.5724 \\
     2 &  820.1398 & 820.19231 & 820.19667 & 820.14775 &  820.1493 \\
     3 & 1008.5656 & 1008.5855 & 1008.6073 & 1008.5752 & 1008.5765 \\
     4 & 1103.4813 & 1103.5053 & 1103.5280 & 1103.4883 & 1103.4902 \\
     5 & 1168.9563 & 1169.2104 & 1169.2531 & 1169.0630 & 1169.0754 \\
    10 & 1324.2906 & 1325.1208 & 1324.6321 & 1324.2975 & 1324.2999 \\
    15 & 1381.8657 & 1383.0802 & 1382.2322 & 1381.8758 & 1381.8738 \\
    20 & 1411.4264 & 1414.4092 & 1411.8528 & 1411.4410 & 1411.4292 \\
    \botrule
  \end{tabular}
\end{table}

\clearpage

\begin{figure}[h]
\begin{center}
\bf  FIGURES\\[1.cm]
\end{center}

  \caption{\protect
    Total CPU time of the fifth CPSCF iteration of fourth order for
    the water cluster sequence with the 6-31G and 6-31G** 
    basis sets and the {\tt GOOD} and {\tt TIGHT} 
    numerical thresholds (see text) controlling numerical
    precision of the result. The lines are fits to the 
    last three and four points, respectively.
  }\label{fig:Gamma_scaling}

  \caption{\protect
    QCTC CPU time of the fifth CPSCF iteration of fourth order for
    the water cluster sequence with the 6-31G and 6-31G** 
    basis sets and the {\tt GOOD} and {\tt TIGHT} 
    numerical thresholds (see text) controlling numerical
    precision of the result. The lines are fits to the 
    last three and four points, respectively.
  }\label{fig:Gamma_QCTC_Timing}

  \caption{\protect
    TC2 CPU time of the fifth CPSCF iteration of fourth order for
    the water cluster sequence with the 6-31G and 6-31G** 
    basis sets and the {\tt GOOD} and {\tt TIGHT} 
    numerical thresholds (see text) controlling numerical
    precision of the result. The lines are fits to the 
    last three and four points, respectively.
  }\label{fig:Gamma_TC2R_Timing}

  \caption{\protect
    ONX CPU time of the fifth CPSCF iteration of fourth order for
    the water cluster sequence with the 6-31G and 6-31G** 
    basis sets and the {\tt GOOD} and {\tt TIGHT} 
    numerical thresholds (see text) controlling numerical
    precision of the result. The lines are fits to the 
    last three and four points, respectively.
  }\label{fig:Gamma_ONX_Timing}

  \caption{Total CPU times with increasing order of the response for 
           the fifth CPSCF cycle computed as the $n+1$ expectation value, Eq.~(\ref{Np1Rule}).}\label{TimeWithOrder}

  \caption{\protect
    Superposition of the magnitudes of the RHF/6-31G density matrix
    derivative elements $D_{cd}$, $D^{x}_{cd}$, $D^{xx}_{cd}$ and $D^{xxx}_{cd}$
    along the $x$ axis with the separation of basis function centers
    for $\rm (H_2O)_{150}$. The density matrix 
    derivatives have been converged to within {\tt TIGHT} (e.g. 
    a matrix threshold $\tau=10^{-6}$ $[a.u.]$).
  }\label{fig:Superposition_Decay}


\end{figure}

\clearpage

\begin{center}
Figure 1, V.~Weber, A.~Niklasson,  and M.~Challacombe \\[1.cm]
\resizebox*{5in}{!}{\includegraphics[angle=-90.00]{figure1.ps}}
\end{center}

\clearpage

\begin{center}
Figure 2, V.~Weber, A.~Niklasson,  and M.~Challacombe \\[1.cm]
\resizebox*{5in}{!}{\includegraphics[angle=-90.00]{figure2.ps}}
\end{center}

\clearpage

\begin{center}
Figure 3, V.~Weber, A.~Niklasson,  and M.~Challacombe \\[1.cm]
\resizebox*{5in}{!}{\includegraphics[angle=-90.00]{figure3.ps}}
\end{center}

\clearpage

\begin{center}
Figure 4, V.~Weber, A.~Niklasson,  and M.~Challacombe \\[1.cm]
\resizebox*{5in}{!}{\includegraphics[angle=-90.00]{figure4.ps}}
\end{center}

\clearpage

\begin{center}
Figure 5, V.~Weber, A.~Niklasson,  and M.~Challacombe \\[1.cm]
\resizebox*{5in}{!}{\includegraphics[angle=-90.00]{figure5.ps}}
\end{center}

\clearpage

\begin{center}
Figure 6, V.~Weber, A.~Niklasson,  and M.~Challacombe \\[1.cm]
\resizebox*{6in}{!}{\includegraphics[angle=-90.00]{figure6.ps}}
\end{center}

}

\end{document}

Response3.bbl
\begin{thebibliography}{54}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname bibnamefont\endcsname\relax
  \def\bibnamefont#1{#1}\fi
\expandafter\ifx\csname bibfnamefont\endcsname\relax
  \def\bibfnamefont#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\providecommand*{\bibinfo}[2]{#2}
\providecommand*{\eprint}[1]{#1}
\providecommand*{\url}[1]{#1}
\begingroup\makeatletter
 \@temptokena{%
  \expandafter\ifx\csname citenamefont\endcsname\relax
   \DeclareRobustCommand\citenamefont{\@firstofone}%
   \global\let\citenamefont\citenamefont
   \global\expandafter\let\csname citenamefont \expandafter\endcsname\csname
  citenamefont \endcsname
  \fi
 }\if@filesw\immediate\write\@auxout{\the\@temptokena}\fi
\expandafter\endgroup\the\@temptokena

\bibitem[{\citenamefont{Galli}(1996)}]{GGalli96}
\bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Galli}}, \bibinfo{journal}{Cur.
  Op. Sol. State Mat. Sci.} \textbf{\bibinfo{volume}{1}}(\bibinfo{number}{6}),
  \bibinfo{pages}{864} (\bibinfo{year}{1996}).

\bibitem[{\citenamefont{Bowler} \emph{et~al.}(1997)\citenamefont{Bowler, Aoki,
  Goringe, Horsfield, and Pettifor}}]{DBowler97}
\bibinfo{author}{\bibfnamefont{D.~R.} \bibnamefont{Bowler}},
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Aoki}},
  \bibinfo{author}{\bibfnamefont{C.~M.} \bibnamefont{Goringe}},
  \bibinfo{author}{\bibfnamefont{A.~P.} \bibnamefont{Horsfield}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{D.~G.}
  \bibnamefont{Pettifor}}, \bibinfo{journal}{Mod. Sim. Mat. Sci. Eng.}
  \textbf{\bibinfo{volume}{5}}(\bibinfo{number}{3}), \bibinfo{pages}{199}
  (\bibinfo{year}{1997}).

\bibitem[{\citenamefont{Goedecker}(1999)}]{SGoedecker99}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Goedecker}},
  \bibinfo{journal}{Rev. Mod. Phys.} \textbf{\bibinfo{volume}{71}},
  \bibinfo{pages}{1085} (\bibinfo{year}{1999}).

\bibitem[{\citenamefont{Ordejon}(2000)}]{POrdejon00}
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Ordejon}},
  \bibinfo{journal}{Phys. Status Solidi B}
  \textbf{\bibinfo{volume}{217}}(\bibinfo{number}{1}), \bibinfo{pages}{335}
  (\bibinfo{year}{2000}).

\bibitem[{\citenamefont{Gogonea} \emph{et~al.}(2001)\citenamefont{Gogonea,
  Suarez, van~der Vaart, and Merz}}]{VGogonea01}
\bibinfo{author}{\bibfnamefont{V.}~\bibnamefont{Gogonea}},
  \bibinfo{author}{\bibfnamefont{D.}~\bibnamefont{Suarez}},
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{van~der Vaart}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{K.~W.} \bibnamefont{Merz}},
  \bibinfo{journal}{Curr. Opin. Struct. Biol.}
  \textbf{\bibinfo{volume}{11}}(\bibinfo{number}{2}), \bibinfo{pages}{217}
  (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Wu and Jayanthi}(2002)}]{SWu02}
\bibinfo{author}{\bibfnamefont{S.~Y.} \bibnamefont{Wu}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{C.~S.} \bibnamefont{Jayanthi}},
  \bibinfo{journal}{Phys. Rep.} \textbf{\bibinfo{volume}{358}},
  \bibinfo{pages}{1} (\bibinfo{year}{2002}).

\bibitem[{\citenamefont{Sekino and Bartlett}(1986)}]{HSekino86}
\bibinfo{author}{\bibfnamefont{H.}~\bibnamefont{Sekino}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{R.~J.} \bibnamefont{Bartlett}},
  \bibinfo{journal}{J. Chem. Phys.} \textbf{\bibinfo{volume}{85}},
  \bibinfo{pages}{976} (\bibinfo{year}{1986}).

\bibitem[{\citenamefont{Karna and Dupuis}(1991)}]{SKarna91}
\bibinfo{author}{\bibfnamefont{S.~P.} \bibnamefont{Karna}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Dupuis}},
  \bibinfo{journal}{J. Comput. Chem.} \textbf{\bibinfo{volume}{12}},
  \bibinfo{pages}{487} (\bibinfo{year}{1991}).

\bibitem[{\citenamefont{Wolinski} \emph{et~al.}(1990)\citenamefont{Wolinski,
  Hinton, and Pulay}}]{KWolinski90}
\bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Wolinski}},
  \bibinfo{author}{\bibfnamefont{J.~F.} \bibnamefont{Hinton}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Pulay}},
  \bibinfo{journal}{jacs} \textbf{\bibinfo{volume}{112}}, \bibinfo{pages}{8251}
  (\bibinfo{year}{1990}).

\bibitem[{\citenamefont{Pennington and Slichter}(1991)}]{CPennington91}
\bibinfo{author}{\bibfnamefont{C.~H.} \bibnamefont{Pennington}}
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{C.~P.}
  \bibnamefont{Slichter}}, \bibinfo{journal}{Phys. Rev. Lett.}
  \textbf{\bibinfo{volume}{66}}, \bibinfo{pages}{381} (\bibinfo{year}{1991}).

\bibitem[{\citenamefont{Malkina} \emph{et~al.}(1996)\citenamefont{Malkina,
  Salahub, and Malkin}}]{OMalkin96}
\bibinfo{author}{\bibfnamefont{O.~L.} \bibnamefont{Malkina}},
  \bibinfo{author}{\bibfnamefont{D.~R.} \bibnamefont{Salahub}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{V.~G.}
  \bibnamefont{Malkin}}, \bibinfo{journal}{J. Chem. Phys.}
  \textbf{\bibinfo{volume}{105}}, \bibinfo{pages}{8793} (\bibinfo{year}{1996}).

\bibitem[{\citenamefont{Amos and Rice}(1989)}]{RAmos89}
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Amos}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.~E.} \bibnamefont{Rice}},
  \bibinfo{journal}{Comp. Phys. Rep.} \textbf{\bibinfo{volume}{10}},
  \bibinfo{pages}{147} (\bibinfo{year}{1989}).

\bibitem[{\citenamefont{Lazzeri and Mauri}(2003)}]{MLazzeri03}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Lazzeri}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Mauri}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{90}},
  \bibinfo{pages}{36401} (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{Quinet and Champagne}(2001)}]{OQuinet01}
\bibinfo{author}{\bibfnamefont{O.}~\bibnamefont{Quinet}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Champagne}},
  \bibinfo{journal}{J. Chem. Phys.} \textbf{\bibinfo{volume}{115}},
  \bibinfo{pages}{6293} (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Pople} \emph{et~al.}(1979)\citenamefont{Pople,
  Krishnan, Schlegel, and Binkley}}]{JPople79}
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Pople}},
  \bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Krishnan}},
  \bibinfo{author}{\bibfnamefont{H.~B.} \bibnamefont{Schlegel}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{J.~S.}
  \bibnamefont{Binkley}}, \bibinfo{journal}{ijqcs}
  \textbf{\bibinfo{volume}{S13}}, \bibinfo{pages}{225} (\bibinfo{year}{1979}).

\bibitem[{\citenamefont{Ochsenfeld and Head-Gordon}(1997)}]{Ochsenfeld97}
\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Ochsenfeld}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Head-Gordon}},
  \bibinfo{journal}{Chem. Phys. Lett.} \textbf{\bibinfo{volume}{270}},
  \bibinfo{pages}{399} (\bibinfo{year}{1997}).

\bibitem[{\citenamefont{Larsen} \emph{et~al.}(2001)\citenamefont{Larsen,
  Helgaker, Olsen, and Jorgensen}}]{HLarsen01a}
\bibinfo{author}{\bibfnamefont{H.}~\bibnamefont{Larsen}},
  \bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Helgaker}},
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Olsen}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Jorgensen}},
  \bibinfo{journal}{J. Chem. Phys.} \textbf{\bibinfo{volume}{115}},
  \bibinfo{pages}{10344} (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Ochsenfeld}
  \emph{et~al.}(2004)\citenamefont{Ochsenfeld, Kussmann, and
  Koziol}}]{COchsenfeld04}
\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Ochsenfeld}},
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Kussmann}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Koziol}},
  \bibinfo{journal}{Angewandte Chemie} \textbf{\bibinfo{volume}{43}},
  \bibinfo{pages}{4485} (\bibinfo{year}{2004}).

\bibitem[{\citenamefont{Brandts}(2001)}]{JBrandts01}
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Brandts}},
  \bibinfo{journal}{Lect. Notes Comp. Sci} \textbf{\bibinfo{volume}{2179}},
  \bibinfo{pages}{462} (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Weber} \emph{et~al.}(2004)\citenamefont{Weber,
  Niklasson, and Challacombe}}]{VWeber04}
\bibinfo{author}{\bibfnamefont{V.}~\bibnamefont{Weber}},
  \bibinfo{author}{\bibfnamefont{A.~M.~N.} \bibnamefont{Niklasson}},
  \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Challacombe}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{92}},
  \bibinfo{pages}{193002} (\bibinfo{year}{2004}).

\bibitem[{\citenamefont{Niklasson and Challacombe}(2004)}]{ANiklasson04}
\bibinfo{author}{\bibfnamefont{A.~M.~N.} \bibnamefont{Niklasson}}
  \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Challacombe}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{92}},
  \bibinfo{pages}{193001} (\bibinfo{year}{2004}).

\bibitem[{\citenamefont{Niklasson}(2002)}]{ANiklasson02A}
\bibinfo{author}{\bibfnamefont{A.~M.~N.} \bibnamefont{Niklasson}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{66}},
  \bibinfo{pages}{155115} (\bibinfo{year}{2002}).

\bibitem[{\citenamefont{Niklasson} \emph{et~al.}(2003)\citenamefont{Niklasson,
  Tymczak, and Challacombe}}]{ANiklasson03}
\bibinfo{author}{\bibfnamefont{A.~M.~N.} \bibnamefont{Niklasson}},
  \bibinfo{author}{\bibfnamefont{C.~J.} \bibnamefont{Tymczak}},
  \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Challacombe}},
  \bibinfo{journal}{J. Chem. Phys.}
  \textbf{\bibinfo{volume}{118}}(\bibinfo{number}{19}), \bibinfo{pages}{8611}
  (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{McWeeny}(1960)}]{RMcWeeny60}
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{McWeeny}},
  \bibinfo{journal}{Rev. Mod. Phys.} \textbf{\bibinfo{volume}{32}},
  \bibinfo{pages}{335} (\bibinfo{year}{1960}).

\bibitem[{\citenamefont{Clinton} \emph{et~al.}(1969)\citenamefont{Clinton,
  Galli, and Massa}}]{WClinton69}
\bibinfo{author}{\bibfnamefont{W.~L.} \bibnamefont{Clinton}},
  \bibinfo{author}{\bibfnamefont{A.~J.} \bibnamefont{Galli}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{L.~J.} \bibnamefont{Massa}},
  \bibinfo{journal}{Phys. Rev.}
  \textbf{\bibinfo{volume}{177}}(\bibinfo{number}{1}), \bibinfo{pages}{7}
  (\bibinfo{year}{1969}).

\bibitem[{\citenamefont{Palser and Manolopoulos}(1998)}]{APalser98}
\bibinfo{author}{\bibfnamefont{A.~H.~R.} \bibnamefont{Palser}}
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{D.~E.}
  \bibnamefont{Manolopoulos}}, \bibinfo{journal}{Phys. Rev. B}
  \textbf{\bibinfo{volume}{58}}, \bibinfo{pages}{12704} (\bibinfo{year}{1998}).

\bibitem[{\citenamefont{Beylkin} \emph{et~al.}(1999)\citenamefont{Beylkin,
  Coult, and Mohlenkamp}}]{GBeylkin99}
\bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Beylkin}},
  \bibinfo{author}{\bibfnamefont{N.}~\bibnamefont{Coult}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.~J.} \bibnamefont{Mohlenkamp}},
  \bibinfo{journal}{J. Comp. Phys.}
  \textbf{\bibinfo{volume}{152}}(\bibinfo{number}{1}), \bibinfo{pages}{32}
  (\bibinfo{year}{1999}).

\bibitem[{\citenamefont{Nemeth and Scuseria}(2000)}]{KNemeth00}
\bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Nemeth}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{G.~E.} \bibnamefont{Scuseria}},
  \bibinfo{journal}{J. Chem. Phys.}
  \textbf{\bibinfo{volume}{113}}(\bibinfo{number}{15}), \bibinfo{pages}{6035}
  (\bibinfo{year}{2000}).

\bibitem[{\citenamefont{Holas}(2001)}]{AHolas01}
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Holas}},
  \bibinfo{journal}{Chem. Phys. Lett.}
  \textbf{\bibinfo{volume}{340}}(\bibinfo{number}{5--6}), \bibinfo{pages}{552}
  (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Niklasson} \emph{et~al.}(2005)\citenamefont{Niklasson,
  Weber, and Challacombe}}]{ANiklasson05a}
\bibinfo{author}{\bibfnamefont{A.~M.~N.} \bibnamefont{Niklasson}},
  \bibinfo{author}{\bibfnamefont{V.}~\bibnamefont{Weber}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Challacombe}},
  \emph{\bibinfo{title}{Density matrix perturbation theory for non-orthogonal
  representations}} (\bibinfo{year}{2005}), \bibinfo{note}{submitted to the
  Journal of Chemical Physics}.

\bibitem[{\citenamefont{Lee and Colwell}(1994)}]{Lee_1994}
\bibinfo{author}{\bibfnamefont{A.~M.} \bibnamefont{Lee}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{S.~M.} \bibnamefont{Colwell}},
  \bibinfo{journal}{J. Chem. Phys.} \textbf{\bibinfo{volume}{101}},
  \bibinfo{pages}{9704} (\bibinfo{year}{1994}).

\bibitem[{\citenamefont{Salek} \emph{et~al.}(2002)\citenamefont{Salek, Vahtras,
  Helgaker, and {\AA}rgen}}]{PSalek02}
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Salek}},
  \bibinfo{author}{\bibfnamefont{O.}~\bibnamefont{Vahtras}},
  \bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Helgaker}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{H.}~\bibnamefont{{\AA}rgen}},
  \bibinfo{journal}{J. Chem. Phys.} \textbf{\bibinfo{volume}{117}},
  \bibinfo{pages}{9630} (\bibinfo{year}{2002}).

\bibitem[{\citenamefont{Wilkinson}(1965)}]{JWilkinson65}
\bibinfo{author}{\bibfnamefont{J.~H.} \bibnamefont{Wilkinson}},
  \emph{\bibinfo{title}{The aglebraic Eigenvalue Problem}}
  (\bibinfo{publisher}{Clarendon Press}, \bibinfo{address}{Oxford},
  \bibinfo{year}{1965}).

\bibitem[{\citenamefont{Stewart}(1973)}]{GStewart73}
\bibinfo{author}{\bibfnamefont{G.~W.} \bibnamefont{Stewart}},
  \emph{\bibinfo{title}{Introduction to Matrix Computations}}
  (\bibinfo{publisher}{Acadmeic Press}, \bibinfo{address}{London},
  \bibinfo{year}{1973}).

\bibitem[{\citenamefont{Benzi and Meyer}(1995)}]{MBenzi95}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Benzi}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{C.~D.} \bibnamefont{Meyer}},
  \bibinfo{journal}{SIAM J. Sci. Comput.}
  \textbf{\bibinfo{volume}{16}}(\bibinfo{number}{5}), \bibinfo{pages}{1159}
  (\bibinfo{year}{1995}).

\bibitem[{\citenamefont{Benzi} \emph{et~al.}(1996)\citenamefont{Benzi, Meyer,
  and T{\accent23 u}ma}}]{MBenzi96}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Benzi}},
  \bibinfo{author}{\bibfnamefont{C.~D.} \bibnamefont{Meyer}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{T{\accent23 u}ma}},
  \bibinfo{journal}{SIAM J. Sci. Comput.}
  \textbf{\bibinfo{volume}{17}}(\bibinfo{number}{5}), \bibinfo{pages}{1135}
  (\bibinfo{year}{1996}).

\bibitem[{\citenamefont{Benzi} \emph{et~al.}(2001)\citenamefont{Benzi, R, and
  Tuma}}]{MBenzi01}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Benzi}},
  \bibinfo{author}{\bibfnamefont{R.~K.} \bibnamefont{R}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Tuma}},
  \bibinfo{journal}{Comp. Meth. App. Mech. Eng.}
  \textbf{\bibinfo{volume}{190}}(\bibinfo{number}{49-50}),
  \bibinfo{pages}{6533} (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Challacombe and Schwegler}(1997)}]{MChallacombe97}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Challacombe}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{E.}~\bibnamefont{Schwegler}},
  \bibinfo{journal}{J. Chem. Phys.} \textbf{\bibinfo{volume}{106}},
  \bibinfo{pages}{5526} (\bibinfo{year}{1997}).

\bibitem[{\citenamefont{Schwegler} \emph{et~al.}(1997)\citenamefont{Schwegler,
  Challacombe, and Head-Gordon}}]{ESchwegler97}
\bibinfo{author}{\bibfnamefont{E.}~\bibnamefont{Schwegler}},
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Challacombe}},
  \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Head-Gordon}},
  \bibinfo{journal}{J. Chem. Phys.} \textbf{\bibinfo{volume}{106}},
  \bibinfo{pages}{9708} (\bibinfo{year}{1997}).

\bibitem[{\citenamefont{Furche}(2001)}]{Furche_2001}
\bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Furche}}, \bibinfo{journal}{J.
  Chem. Phys.} \textbf{\bibinfo{volume}{114}}, \bibinfo{pages}{5982}
  (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Challacombe}(2000)}]{MChallacombe00B}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Challacombe}},
  \bibinfo{journal}{Comput. Phys. Commun.} \textbf{\bibinfo{volume}{128}},
  \bibinfo{pages}{93} (\bibinfo{year}{2000}).

\bibitem[{\citenamefont{Weber and Daul}(2003)}]{VWeber03}
\bibinfo{author}{\bibfnamefont{V.}~\bibnamefont{Weber}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Daul}},
  \bibinfo{journal}{Chem. Phys. Lett.} \textbf{\bibinfo{volume}{370}},
  \bibinfo{pages}{99} (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{Liang} \emph{et~al.}(2004)\citenamefont{Liang, Baer,
  Saravanan, Shao, Bell, and Head-Gordon}}]{WLiang04}
\bibinfo{author}{\bibfnamefont{W.~Z.} \bibnamefont{Liang}},
  \bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Baer}},
  \bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Saravanan}},
  \bibinfo{author}{\bibfnamefont{Y.~H.} \bibnamefont{Shao}},
  \bibinfo{author}{\bibfnamefont{A.~T.} \bibnamefont{Bell}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Head-Gordon}},
  \bibinfo{journal}{J. Comp. Phys.} \textbf{\bibinfo{volume}{194}},
  \bibinfo{pages}{575} (\bibinfo{year}{2004}).

\bibitem[{\citenamefont{Voter} \emph{et~al.}(1996)\citenamefont{Voter, Kress,
  and Silver}}]{AVoter96}
\bibinfo{author}{\bibfnamefont{A.~F.} \bibnamefont{Voter}},
  \bibinfo{author}{\bibfnamefont{J.~D.} \bibnamefont{Kress}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{R.~N.} \bibnamefont{Silver}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{53}},
  \bibinfo{pages}{12733} (\bibinfo{year}{1996}).

\bibitem[{\citenamefont{Pulay}(1980)}]{Pulay80}
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Pulay}},
  \bibinfo{journal}{Chem. Phys. Let.}
  \textbf{\bibinfo{volume}{73}}(\bibinfo{number}{2}), \bibinfo{pages}{393}
  (\bibinfo{year}{1980}).

\bibitem[{\citenamefont{Pulay}(1982)}]{Pulay82}
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Pulay}}, \bibinfo{journal}{J.
  Comput. Chem.} \textbf{\bibinfo{volume}{3}}(\bibinfo{number}{4}),
  \bibinfo{pages}{556} (\bibinfo{year}{1982}).

\bibitem[{\citenamefont{Epstein}(1974)}]{SEpstein74}
\bibinfo{author}{\bibfnamefont{S.~T.} \bibnamefont{Epstein}},
  \emph{\bibinfo{title}{The Variation Method in Quantum Chemistry}}
  (\bibinfo{publisher}{Academic Press}, \bibinfo{address}{San Fransisco},
  \bibinfo{year}{1974}).

\bibitem[{\citenamefont{McWeeny}(1962)}]{RMcWeeny62}
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{McWeeny}},
  \bibinfo{journal}{Rev. Mod. Phys.} \textbf{\bibinfo{volume}{126}},
  \bibinfo{pages}{1028} (\bibinfo{year}{1962}).

\bibitem[{\citenamefont{Schmidt} \emph{et~al.}(1993)\citenamefont{Schmidt,
  Baldridge, Boatz, Elbert, Gordon, Jensen, Koseki, Matsunaga, Nguyen, Su,
  Windus, Dupuis} \emph{et~al.}}]{gamess}
\bibinfo{author}{\bibfnamefont{M.~W.} \bibnamefont{Schmidt}},
  \bibinfo{author}{\bibfnamefont{K.~K.} \bibnamefont{Baldridge}},
  \bibinfo{author}{\bibfnamefont{J.~A.} \bibnamefont{Boatz}},
  \bibinfo{author}{\bibfnamefont{S.~T.} \bibnamefont{Elbert}},
  \bibinfo{author}{\bibfnamefont{M.~S.} \bibnamefont{Gordon}},
  \bibinfo{author}{\bibfnamefont{J.~H.} \bibnamefont{Jensen}},
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Koseki}},
  \bibinfo{author}{\bibfnamefont{N.}~\bibnamefont{Matsunaga}},
  \bibinfo{author}{\bibfnamefont{K.~A.} \bibnamefont{Nguyen}},
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Su}},
  \bibinfo{author}{\bibfnamefont{T.~L.} \bibnamefont{Windus}},
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Dupuis}}, \emph{et~al.},
  \bibinfo{journal}{J. Comp. Chem.} \textbf{\bibinfo{volume}{14}},
  \bibinfo{pages}{1347} (\bibinfo{year}{1993}).

\bibitem[{\citenamefont{Challacombe}
  \emph{et~al.}(2001)\citenamefont{Challacombe, Schwegler, Tymczak, Gan,
  Nemeth, Weber, Niklasson, and Henkelman}}]{MondoSCF}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Challacombe}},
  \bibinfo{author}{\bibfnamefont{E.}~\bibnamefont{Schwegler}},
  \bibinfo{author}{\bibfnamefont{C.~J.} \bibnamefont{Tymczak}},
  \bibinfo{author}{\bibfnamefont{C.~K.} \bibnamefont{Gan}},
  \bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Nemeth}},
  \bibinfo{author}{\bibfnamefont{V.}~\bibnamefont{Weber}},
  \bibinfo{author}{\bibfnamefont{A.~M.~N.} \bibnamefont{Niklasson}},
  \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Henkelman}},
  \emph{\bibinfo{title}{{\sc MondoSCF} v1.0$\alpha$9, A program suite for
  massively parallel, linear scaling SCF theory and {\it ab initio} molecular
  dynamics.}} (\bibinfo{year}{2001}),
  \urlprefix\url{http://www.t12.lanl.gov/home/mchalla/},
  \bibinfo{note}{\mbox{L}os Alamos National Laboratory (LA-CC 01-2), Copyright
  University of California.}

\bibitem[{\citenamefont{Tymczak and Challacombe}(2004)}]{CTymczak04a}
\bibinfo{author}{\bibfnamefont{C.~J.} \bibnamefont{Tymczak}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Challacombe}},
  \emph{\bibinfo{title}{Linear scaling computation of the Fock matrix. VIII.
  Periodic Density Functional Theory at the $\Gamma$-point}}
  (\bibinfo{year}{2004}), \bibinfo{note}{cond-mat/0405500, To appear in
  J.~Chem.~Phys.}

\bibitem[{\citenamefont{Tymczak} \emph{et~al.}(2004)\citenamefont{Tymczak,
  Weber, Schwegler, and Challacombe}}]{CTymczak04b}
\bibinfo{author}{\bibfnamefont{C.~J.} \bibnamefont{Tymczak}},
  \bibinfo{author}{\bibfnamefont{V.}~\bibnamefont{Weber}},
  \bibinfo{author}{\bibfnamefont{E.}~\bibnamefont{Schwegler}},
  \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Challacombe}},
  \emph{\bibinfo{title}{Linear scaling computation of the Fock matrix. IX.
  Periodic boundaries for exact exchange at the $\Gamma$-point}}
  (\bibinfo{year}{2004}), \bibinfo{note}{cond-mat/0406094, To appear in
  J.~Chem.~Phys.}

\bibitem[{\citenamefont{\mbox{The Portland Group}}(2002)}]{PGF90}
\bibinfo{author}{~\bibnamefont{\mbox{The Portland Group}}},
  \emph{\bibinfo{title}{pgf90 v4.2}} (\bibinfo{year}{2002}),
  \urlprefix\url{http://www.pgroup.com/}.

\bibitem[{\citenamefont{Otto} \emph{et~al.}(1999)\citenamefont{Otto, Gu, and
  Ladik}}]{POtto99}
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Otto}},
  \bibinfo{author}{\bibfnamefont{F.~L.} \bibnamefont{Gu}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Ladik}}, \bibinfo{journal}{J.
  Chem. Phys.} \textbf{\bibinfo{volume}{110}}, \bibinfo{pages}{2717}
  (\bibinfo{year}{1999}).

\end{thebibliography}
