%==================
\section{
\label{sec:algorithm}
Algorithm}
%==================
In this section, we present our recursive algorithm for the computation of a class of three- or four-electron integrals of arbitrary angular momentum.
The present recursive algorithm is based on a late-contraction scheme inspired by the HGP algorithm \cite{HGP} following a BOVVVVCCCCHHHH path.
The general skeleton of the algorithm is shown in Fig.~\ref{fig:algo} for the trident operator $C_{12}G_{13}G_{14}$.
We will use this example to illustrate each step.

Based on the shell data, the first step of the algorithm (step B) is to decide whether or not a given class of integrals is significant or negligible.
If the integral class is found to be significant by the screening algorithm depicted in Fig.~\ref{fig:scheme1}, an initial set of FIs is computed (step O).

Starting with these FIs, angular momentum is then built up over the different bra centres $\bA_1$, $\bA_2$, $\bA_3$ and $\bA_4$ using the VRRs derived in the previous section.
To minimise the computational cost, one has to think carefully how to perform this step.
Indeed, the cost depends on the order in which this increase in angular momentum is performed.
This is illustrated in Fig.~\ref{fig:graph}, where we have represented the various possible pathways for the 3-chain operator $C_{12}G_{23}$ (left) and the trident operator $C_{12}G_{13}G_{14}$ (right).
The red path corresponds to the path generating the least intermediates (i.e.~requiring the smallest number of classes in order to compute a given class).
Different paths are compared in Table \ref{tab:RRcount} for various two-, three- and four-electron operators, where we have reported the number of intermediates generated by each path for various integral classes.

Taking the 3-chain operator $C_{12}G_{23}$ as an example, one can see that, to compute a $\sexpval{ppp}$ class, it is more advantageous to build momentum over center $\bA_3$, then over centres $\bA_2$, and finally over center $\bA_1$ using VRRs with 4, 6 and 6 terms, respectively.
The alternative path corresponding to building momentum over $\bA_3$, $\bA_1$, and then $\bA_2$ with 4-, 5- and 7-term VRRs is slightly more expensive for a $\sexpval{ppp}$ class but becomes affordable for high angular momentum classes.
For both paths, using the TRR instead of the last VRR implies a large increase in the number of intermediates.

For the trident operator, we successively build angular momentum over $\bA_4$, $\bA_3$, $\bA_1$ and $\bA_2$ using VRRs with 4, 6, 8 and 7 terms.
The pathway using VRRs with 4, 6, 6, and 9 terms is more expensive due to the large number of terms of the VRR building up momentum over the last center.
Again, using the TRR instead of the last VRR significantly increases the number of intermediates.

The path involving the minimal number of intermediates is given in Table \ref{tab:RRcount} for various two-, three- and four-electron operators. 
It is interesting to point out that it is never beneficial to use the TRR derived in Eq.~\eqref{eq:TRR}.

One can easily show that, for operators involving the Coulomb operator, the number of intermediates required to compute a $n$-electron integral class $\sexpval{a \ldots a}$ increases as $\order{a^{n+1}}$ for the VRR-only paths (see Table \ref{tab:RRcount}).
This number is reduced to $\order{a^{n}}$ if one uses the TRR to build up angular momentum on the last center. 
However, the prefactor is much larger and the crossover happens for extremely high angular momentum for three- and four-electron integrals.
For ``pure''  GG operators, such as $G_{12}$ or $G_{13}G_{23}$, the number of intermediates required to compute a class $\sexpval{a \ldots a}$ increases as $\order{a^n}$ for any type of paths.

Finally, we note that the optimal path for the trident $C_{12}G_{13}G_{14}$ and the 4-chain $C_{12}G_{13}G_{34}$ is similar, thanks to their similar structure.
Indeed, these two operators can be seen as two ``linked'' GGs ($G_{13}G_{14}$ or $G_{13}G_{34}$) interacting with the Coulomb operator $C_{12}$ (see Fig.~\ref{fig:tree}), while the other 4-chain operator $C_{12}G_{14}G_{23}$ can be seen as two ``unlinked'' GGs ($G_{14}$ and $G_{23}$) interacting with the Coulomb operator.

When angular momentum has been built over all the bra centres, following the HGP algorithm \cite{HGP}, we contract $\sbraket{\ba_1 \ba_2 \ba_3  \ba_4}{\bo \bo \bo \bo}$ to form $\braket{\ba_1 \ba_2 \ba_3  \ba_4}{\bo \bo \bo \bo}$ (step CCCC).
We can perform the contraction at this point because all of the subsequent RRs are independent of the contraction coefficients and exponents.
More details about this contraction step can be found in Ref.~\cite{Gill94b}.

The last step of the algorithm (step HHHH) shifts momentum from the bra center $\bA_1$, $\bA_2$, $\bA_3$ and $\bA_4$ to the ket centres  $\bB_1$, $\bB_2$, $\bB_3$ and $\bB_4$ using the two-term HRRs given by Eq.~\eqref{eq:HRR}.

%%% FIGURE 3 %%%
\begin{figure}
	\centering
	\includegraphics[width=0.6\linewidth]{../Chapter4/fig/fig3}
\caption{
\label{fig:algo}
PRISM representation \cite{PRISM91} of the recursive algorithm used to compute a four-electron integral class $\braket{a_1 a_2 a_3 a_4}{b_1 b_2 b_3 b_4}$ over the trident operator $C_{12}G_{13}G_{14}$.
In our algorithm, we consider the (orange) BOVVVVCCCCHHHH path.}
\end{figure}
%%%

%%% FIGURE 4 %%%
\begin{figure*}
	\centering
	\includegraphics[height=0.38\textheight]{../Chapter4/fig/fig4a}
	\includegraphics[height=0.38\textheight]{../Chapter4/fig/fig4b}
\caption{
\label{fig:graph}
Graph representation of the VRRs for the 3-chain operator $C_{12}G_{23}$ (left) and trident operator $C_{12}G_{13} G_{14}$ (right).
The edge label gives the number of terms in the corresponding VRR.
The red path corresponds to the algorithm generating the smallest number of intermediates.}
\end{figure*}
%%%
