\documentstyle[11pt]{report}
\pagestyle{empty}
\textwidth = 16 cm
\textheight = 25 cm
\topmargin =  -1.0 cm
\oddsidemargin = 0.0 cm
\begin{document}
{\bf Generalized problem}.

\par\noindent
The problem: minimize
\begin{equation}
f_j={1\over 2} <\psi_j|H|\psi_j> 
\end{equation}
under the constraints:
\begin{equation}
<\psi_j|S|\psi_j>=1, \qquad <\psi_j|S|\psi_k>=0,\quad k<j
\end{equation}
%The alghoritm used to solve the problem is exactly the same as the one
%used for "normal" constrains, apart from
%the insertion of the S matrix in the right places.  
%
% l'algoritmo e` lo stesso perche` appunto era gia` implementato con il
% precondizionamento che tiene conto degli accorgimenti particolari da
% fare sull'ortogonalizzazione, pero` la scrittura di questa prima parte
% non e` vero che e` la stessa del caso senza overlap
%
%
{\bf bisogna fare come per il problema con il precondizionamento perche`
$<\psi^0_j|{\bf S^2}|\psi_k> \ne 0$ }
...This is equivalent to solve the Euler-Lagrange problem of minimising
\begin{equation}
\widetilde f_j={1\over 2} <\psi_j|H|\psi_j> -
  \epsilon_j (<\psi_j|S|\psi_j>-1) -
  \sum_{k<j} \lambda_{jk} <\psi_j|S|\psi_k>
\end{equation}
If we consider the gradient of the expression above we have:
\begin{equation}\label{prima}
|g_j>=H|\psi_j> - \epsilon_j S|\psi_j> - 
       \sum_{k<j} \lambda_{jk} S|\psi_k>
\end{equation}
but we are not able to find a simple expression for $\epsilon_j$ and
$\lambda_{jk}$ in order to have $<g_j|S|\psi_k>=0$ because
$<\psi_j|S^2|\psi_k> \ne 0$.  \\
So we use a slightly different approach:
we firstly consider only the first constraint. The problem is
equivalent to the minimization of

\begin{equation}
\widetilde f_j={1\over 2} <\psi_j|H|\psi_j> -
  \epsilon_j (<\psi_j|S|\psi_j>-1) 
\end{equation}
{\bf First step.}
For any given index $j$ we assume a starting vector 
$|\psi^0_j>$ such that
\begin{equation}
<\psi^0_j|S|\psi^0_j>=1,\qquad <\psi^0_j|S|\psi_k>=0, \quad k<j
\end{equation}
The gradient $|g^0_j>$ is
\begin{equation}
|g^0_j>=H|\psi^0_j> - \epsilon^{0^*}_j S|\psi^0_j>  
\end{equation}
where
\begin{equation}
\epsilon^{0^*}_j = { <\psi^0_j|SH|\psi^0_j> \over <\psi^0_j|S^2|\psi^0_j>}
\end{equation}
ensures that $<g^0_j|S|\psi^0_j>=0$. We also define:
\begin{equation}
\epsilon^{0}_j = <\psi^0_j|H|\psi^0_j>
\end{equation}
which is an approximation for the true eigenvalue $\epsilon_j$ of H.
The second set of constraints is accounted requiring that
\begin{equation}
<g^0_j|S|\psi_k>=0
\end{equation}
This is obtained by an explicit Gram-Shmidt orthogonalization, with the
overlap matrix S, of $|g^0_j>$
to $|\psi_k>, k<j $
\begin{equation} \label{seconda}
|g^0_j> = H|\psi^0_j> - \epsilon^{0^*}_j S|\psi^0_j> - 
\sum_{k<j}\lambda^0_{kj} |\psi_k>
\end{equation}
where
\begin{equation}
\lambda^0_{kj}=<g^0_j|S|\psi_k>
\end{equation}
Note that equation \ref{seconda} is very much different from equation
\ref{prima}, the last term of the right hand side is 
$ \sum_{k<j}\lambda^0_{kj} |\psi_k> $ instead of 
$\sum_{k<j}\lambda^0_{kj} S|\psi_k>$.\\
Now: let us consider the conjugate direction $|h^0_j>=|g^0_j>$ and the 
normalized direction $|\widetilde h^0_j>$:
\begin{equation}
|\widetilde h^0_j> = {|h^0_j>\over <h^0_j|S|h^0_j>^{1/2}}
\end{equation}
By construction
\begin{equation}
<\widetilde h^0_j|S|\widetilde h^0_j>=1,  \quad 
<\widetilde h^0_j|S|\psi^0_j>=0,  \quad 
<\widetilde h^0_j|S|\psi_k>=0  \quad 
\end{equation}
Consider the vector
\begin{equation}
|\psi^1_j(\theta)>= \cos\theta |\psi^0_j> + \sin\theta|\widetilde h^0_j>
\end{equation}
By construction,
\begin{equation}
<\psi^1_j(\theta)|S|\psi^1_j(\theta)>=1, \qquad 
<\psi^1_j(\theta)|S|\psi_k>=0,\quad k<j 
\end{equation}
We look for the value of $\theta$ that minimizes 
\begin{equation}
\epsilon^1_j(\theta) = <\psi^1_j(\theta)|H|\psi^1_j(\theta)>
\end{equation}
The solution is
\begin{equation}
\theta^0_j={1\over 2} \mbox{arctag}\left({a^0_j\over\epsilon^0_j-b^0_j}\right)
\end{equation}
where
\begin{equation}
a^0_j=< \psi^0_j|H|\widetilde h^0_j> + c.c., \quad
b^0_j=< \widetilde h^0_j|H|\widetilde h^0_j> 
\end{equation}
and 
\begin{equation}
\epsilon^1_j = {\epsilon^0_j+b^0_j - \sqrt{(\epsilon^0_j-b^0_j)^2
                                     +(a^0_j)^2}\over 2}
\end{equation}
which is a new approximation for $\epsilon_j$.\\
{\bf Following steps.} We construct the new gradient
$|g^n_j>$ as above and the new conjugate direction
$|h^n_j>$ as in the conjugate gradient method. We define
\begin{equation}
|u^n_j>=|g^n_j>+\gamma^{n-1}_j|h^{n-1}_j>
\end{equation}
\begin{equation}
\gamma^{n-1}_j={<g^n_j|S|g^n_j>\over<g^{n-1}_j|S|g^{n-1}_j>}
\simeq{<g^n_j-g^{n-1}_j|S|g^n_j>\over<g^{n-1}_j|S|g^{n-1}_j>}
\end{equation}
By construction,
\begin{equation}
<g^n_j|S|\psi_k>=0 \longrightarrow <u^n_j|S|\psi_k>=0
\end{equation}
but
\begin{equation}
\beta^n_j= <u^n_j|S|\psi^n_j>=\gamma^{n-1}_j<h^{n-1}_j|S|\psi^n_j>
                =\gamma^{n-1}_j\sin\theta^{n-1}_j
                     <h^{n-1}_j|S|\widetilde h^{n-1}_j> \neq 0
\end{equation}
so that explicit orthogonalization of $|\psi^n_j>$ to $|u^n_j>$ is
needed (it is a consequence of the constraint of unitary norm):
\begin{equation}
|h^n_j>=|u^n_j> - \beta^n_j |\psi^n_j>
\end{equation}
Then one proceeds as above, finds $\epsilon^n_j$ and check if
$\epsilon^n_j - \epsilon^{n-1}_j < tol$, where $tol$ is a small number
which defines the precision of $\epsilon_j$. 
\par\bigskip
{\bf Preconditioning.} We define a (diagonal) preconditioning matrix
$P$ and auxiliary functions $|y>$ :
\begin{equation}
|y_j>=P^{-1}|\psi_j>
\end{equation}
One has to solve the equivalent problem of minimizing
\begin{equation}
f_j={1\over 2} <y_j|(PHP)|y_j>
\end{equation}
under the constraints that
\begin{equation}
<y_j|(PSP)|y_j>=1, \qquad <y_j|PS|\psi_k>=0,\quad k<j
\end{equation}
Again we firstly consider only the first constraint. The problem is equivalent
to the minimization of
\begin{equation}
\widetilde f_j={1\over 2} <y_j|(PHP)|y_j> - 
               \widetilde\epsilon_j (<y_j|(PSP)|y_j>-1)
\end{equation}
The gradient $|g_j>$ is now given by
\begin{equation}
|g_j>=(PHP)|y_j> - \widetilde\epsilon_j(PSP)|y_j>
\end{equation}
where
\begin{equation}
\widetilde\epsilon_j = {<y_j|(PSP)(PHP)|y_j> \over <y_j|(PSP)^2| y_j>}
\end{equation}
ensures that
\begin{equation}
<g_j|(PSP)|y_j> = 0
\end{equation}
Now we have to care about the second set of constraints. We want that
\begin{equation}
<g_j|PS|\psi_k> = 0
\end{equation}
And again this is obtained by orthogonalizing, with the overlap matrix $S$,
$|Pg_j>$ to $|\psi_k> ,k<j$:
\begin{equation}
|Pg^0_j> = |Pg^0_j> - \sum_{k<j} <g^0_j|PS|\psi_k>|\psi_k>
\end{equation}

{\bf Modified preconditioned algorithm}. Starting from an initial
guess $|\psi^0_j>$ (where $<\psi_k|S|\psi^0_j>=0$, 
$<\psi^0_j|S|\psi^0_j>=1$), one searches for new vector as:
\begin{equation}
|\psi^n_j(\theta)>= \cos\theta |\psi^{n-1}_j> + 
                    \sin\theta |P\widetilde h^{n-1}_j>
\end{equation}
where
\begin{equation}
|P\widetilde h^n_j> = {|Ph^n_j>\over <Ph^n_j|S|Ph^n_j>^{1/2}}
\end{equation}
and $|Ph^n_j>$ is the conjugate gradient. We define
\begin{equation}
|Pu^n_j>=|Pg^n_j>+\gamma^{n-1}_j|Ph^{n-1}_j>,
\end{equation}
\begin{equation}
\gamma^{n-1}_j={<Pg^n_j|P^{-2}S|Pg^n_j>
\over<Pg^{n-1}_j|P^{-2}S|Pg^{n-1}_j>}
\simeq{<Pg^n_j-Pg^{n-1}_j|P^{-2}S|Pg^n_j>
\over<Pg^{n-1}_j|P^{-2}S|Pg^{n-1}_j>}
\end{equation}
with $\gamma^{-1}_j=0$.
The gradient $|Pg^n_j>$ is given by
\begin{equation}
|Pg^n_j>=P^2H|\psi^n_j> - \widetilde\epsilon^n_j SP^2|\psi^n_j>- 
       \sum_{k<j}\lambda_{jk}^n|\psi_k>
\end{equation}
where
\begin{equation}
\lambda_{jk}^n = <\psi_k|S(P^2H-\widetilde\epsilon^n_j SP^2)|\psi^n_j>,
\qquad
\widetilde\epsilon_j^n = {<\psi^n_j| SP^2 H|\psi^n_j> \over <\psi^n_j|
S^2P^2| \psi^n_j>} 
\end{equation}
$|Pu^n_j>$ is by construction orthogonal to 
$|\psi_k>$ but not to  $|\psi^n_j>$:
\begin{equation}
\beta^n_j= <Pu^n_j|S|\psi^n_j>=
\gamma^{n-1}_j<Ph^{n-1}_j|S|\psi^n_j>
                            =\gamma^{n-1}_j\sin\theta^{n-1}_j
                            <Ph^{n-1}_j|S|Ph^{n-1}_j> \neq 0
\end{equation}
One gets the conjugate direction $|Ph^n_j>$ from
\begin{equation}
|Ph^n_j>=|Pu^n_j> - \beta^n_j |\psi^n_j>
\end{equation}
Note that only $P^2$ (or $P^{-2}$) is actually used.

\end{document}