\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{CJKutf8} 
\usepackage[utf8]{inputenc}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[english]{babel}



\usepackage{float}






%%% Your header encountered errors when converting with LaTeXML
%%% Please uncomment it incrementally and check for problems
%%%   If you need assistance, contact us at help@authorea.com
%%% 
%
%
%
%\iflatexml
%\documentclass[12pt]{article}
%\fi
%\setlength{\oddsidemargin}{0in}  \setlength{\textwidth}{6.5in}    \setlength{\topmargin}{-0.5in}    \setlength{\textheight}{9in}     \newenvironment{wileykeywords}{\textsf{Keywords:}\hspace{\stretch{1}}}{\hspace{\stretch{1}}\rule{1ex}{1ex}}
%\usepackage{amsmath,amssymb,bm}
%\usepackage{amsmath,amssymb,bm}
%\usepackage{amsmath,amssymb,bm}
%
%\usepackage{graphicx}
%% Include figure files
%\usepackage{color}
%% Include colors for document elements
%\usepackage{dcolumn}
%% Align table columns on decimal point
%\usepackage{bm}
%% bold math
%\usepackage[numbers,super,comma,sort&compress]{natbib}
%\usepackage{multirow}
%
%\usepackage{setspace}
%
%% \usepackage[version=3]{mhchem} % Not (yet) supported by LaTeXML.
%
%% \usepackage{natmove} % Not (yet) supported by LaTeXML.
%
%\usepackage{float}
%
%\usepackage{placeins}
%
%\usepackage{graphics}
%
%\usepackage{rotating}
%
%% \usepackage{sidecap} % Not (yet) supported by LaTeXML.
%
%\usepackage{wrapfig}
%
%\usepackage{verbatim}
%
%% \usepackage{varioref} % Not (yet) supported by LaTeXML.
%
%\usepackage{subfigure}
%
%% \usepackage{ifthen} % Not (yet) supported by LaTeXML.
%
%% \usepackage{commath} % Not (yet) supported by LaTeXML.
%
%\usepackage{enumerate}
%
%\usepackage{natbib}
%
%% \usepackage[T1]{fontenc} % Not (yet) supported by LaTeXML.
%       % Use modern font encodings
%% \usepackage{axodraw4j} % Not (yet) supported by LaTeXML.
%
%\usepackage{tikz}
%
%% First macro on next line not (yet) supported by LaTeXML: 
%% \usetikzlibrary{backgrounds}
%% First macro on next line not (yet) supported by LaTeXML: 
%% \usetikzlibrary{shapes,arrows}
%\usepackage{pifont}
%% http://ctan.org/pkg/pifont
%% First macro on next line not (yet) supported by LaTeXML: 
%% \tikzstyle{decision} = [diamond, draw, fill=blue!20, text width=4.5em, text badly centered, node distance=3cm, inner sep=0pt]
%% First macro on next line not (yet) supported by LaTeXML: 
%% \tikzstyle{block} = [rectangle, draw, fill=blue!20, text width=5em, text centered, rounded corners, minimum height=4em]
%% First macro on next line not (yet) supported by LaTeXML: 
%% \tikzstyle{line} = [draw, -latex']
%% First macro on next line not (yet) supported by LaTeXML: 
%% \tikzstyle{cloud} = [draw, ellipse,fill=red!20, node distance=3cm,minimum height=2em]
%\definecolor{mydarkgreen}{rgb}{0.0, 0.5, 0.0}
%\definecolor{myred}{rgb}{0.8, 0.0, 0.0}
%\definecolor{myamber}{rgb}{1.0, 0.49, 0.0}
%\definecolor{myauburn}{rgb}{0.43, 0.21, 0.1}
%\definecolor{mybrown}{rgb}{0.65, 0.16, 0.16}
%\definecolor{mycapri}{rgb}{0.0, 0.75, 1.0}
%\definecolor{mytang}{rgb}{1.0, 0.6, 0.4}
%\definecolor{myrose}{rgb}{1.0, 0.01, 0.24}
%\definecolor{myballblue}{rgb}{0.13, 0.67, 0.8}
%\definecolor{mycobalt}{rgb}{0.0, 0.28, 0.67}
%\newcommand*{\mycommand}[1][]{\texttt{\emph{#1}}}
%
%\newcommand{\dE}[][]{$\Delta E$}
%
%\newcommand{\kcal}[][]{kcal mol$^{\textrm{-}1}$}
%
%\newcommand{\cm}[][]{cm$^{\textrm{-}1}$}
%
%\newcommand{\cxz}[][]{\hbox{cc-pCV$X$Z}}
%
%\newcommand{\xz}[][]{\hbox{cc-pV$X$Z}}
%
%\newcommand{\dz}[][]{\hbox{cc-pVDZ}}
%
%\newcommand{\tz}[][]{\hbox{cc-pVTZ}}
%
%\newcommand{\qz}[][]{\hbox{cc-pVQZ}}
%
%\newcommand{\atz}[][]{\hbox{aug-cc-pVTZ}}
%
%\newcommand{\aqz}[][]{\hbox{aug-cc-pVQZ}}
%
%\newcommand{\apz}[][]{\hbox{aug-cc-pV5Z}}
%
%\newcommand{\pz}[][]{\hbox{cc-pV5Z}}
%
%\newcommand{\hz}[][]{\hbox{cc-pV6Z}}
%
%\newcommand{\cdz}[][]{\hbox{cc-pCVDZ}}
%
%\newcommand{\ctz}[][]{\hbox{cc-pCVTZ}}
%
%\newcommand{\cqz}[][]{\hbox{cc-pCVQZ}}
%
%\newcommand{\cpz}[][]{\hbox{cc-pCV5Z}}
%
%\newcommand{\poptzd}[][]{\hbox{6-311G(d)}}
%
%\newcommand{\poptzdp}[][]{\hbox{6-311G(d,p)}}
%
%\newcommand{\cas}[][]{CASSCF}
%
%\newcommand{\mr}[][]{MRMP2}
%
%\newcommand{\mrcas}[][]{MRMP2//CASSCF}
%
%\newcommand{\mrcasbs}[][]{MRMP2/6-311G(d,p)//CASSCF/6-311G(d,p)}
%
%\newcommand{\cc}[][]{CCSD(T)}
%
%\newcommand{\ccdft}[][]{CCSD(T)//B3LYP}
%
%\newcommand{\ccdftbs}[][]{CCSD(T)/6-311G(d,p)//B3LYP/6-311G(d,p)}
%
%\newcommand{\dt}[][]{\mbox{\large\bfseries $\cdot$}}
%\definecolor{background-color}{gray}{0.98}
%\renewcommand\@[][]{}
%\renewcommand{\baselinestretch}[][]{1.5}
%\providecommand{\noopsort}[1][]{}
%\providecommand{\singleletter}[1][]{#1}


\begin{document}
\begin{CJK}{UTF8}{gbsn}

\title{Molint 1.0: An Application Programming Interface Framework for the
Computation of Molecular Integrals and Their First Derivatives for the
Density-Fitted Methods}



\author[1]{Ugur Bozkaya}%
\affil[1]{Hacettepe Universitesi}%


\vspace{-1em}



  
  \date{April 22, 2026}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
The efficient computation of molecular integrals and their derivatives is a crucial step in molecular property evaluation in modern quantum chemistry. As an integral tensor decomposition technique,  the density-fitting (DF) approach becomes a popular tool to reduce the memory and disk requirements for the electron repulsion integrals. In this study, an application programming interface (API) framework, denoted (MFW),  for the computation of molecular integrals and their first derivatives, over contracted Gaussian functions, for the density-fitted methods is reported. The MFW is free software and it includes overlap, dipole, kinetic, potential, metric, and 3-index integrals, and their first derivatives. Furthermore, the MFW provides a smooth approach to build the Fock matrix and evaluate analytic gradients for the density-fitted methods. The MFW is a {\sc C++/Fortran} hybrid code, which can take advantage of shared-memory parallel programming techniques. Our results demonstrate that the MFW is an efficient and user-friendly API for the computation of molecular integrals and their first derivatives.%
\end{abstract}%



\sloppy


\begin{wileykeywords}
Electron repulsion integrals, Density-Fitting, Analytic Gradients

\end{wileykeywords}

\clearpage

\makeatletter
biblabel[1]{#1.}
\makeatother

\normalsize

\clearpage

\section*{\sffamily \Large Introduction} % Not needed for rapid communications
Computation of molecular integrals is a crucial step of any molecular property evaluation in modern quantum chemistry.\hyperref[csl:1]{(\bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Helgaker}}, n.d.}; \hyperref[csl:2]{\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Dupuis}}, n.d.}; \hyperref[csl:3]{\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Rys}}, n.d.}; \hyperref[csl:4]{\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Dupuis}} \bibnamefont{and}, n.d.}; \hyperref[csl:5]{\bibinfo{author}{\bibfnamefont{L.~E.} \bibnamefont{McMurchie}}, n.d.}; \hyperref[csl:6]{\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Obara}} \bibnamefont{and}, n.d.}; \hyperref[csl:7]{\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Obara}} \bibnamefont{and}, n.d.}; \hyperref[csl:8]{\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Head-Gordon}} \bibnamefont{and}, n.d.}; \hyperref[csl:9]{\bibinfo{author}{\bibfnamefont{P.~M.~W.} \bibnamefont{Gill}}, n.d.}; \hyperref[csl:10]{\bibinfo{author}{\bibfnamefont{T.~P.} \bibnamefont{Hamilton}} \bibnamefont{and}, n.d.}; \hyperref[csl:11]{\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Lindh}}, n.d.}; \hyperref[csl:12]{\bibinfo{author}{\bibfnamefont{P.~M.~W.} \bibnamefont{Gill}} \bibnamefont{and}, n.d.}; \hyperref[csl:13]{\bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Ishida}}, 5198}; \hyperref[csl:14]{\bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Ishida}}, n.d.}; \hyperref[csl:15]{\bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Ishida}}, n.d.}; \hyperref[csl:16]{\bibinfo{author}{\bibfnamefont{J.~P.} \bibnamefont{Kenny}}, n.d.}; \hyperref[csl:17]{\bibinfo{author}{\bibfnamefont{N.}~\bibnamefont{Flocke}} \bibnamefont{and}, n.d.}; \hyperref[csl:18]{\bibinfo{author}{\bibfnamefont{P.~M.~W.} \bibnamefont{Gill}}, n.d.}; \hyperref[csl:19]{\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Reine}}, n.d.}; \hyperref[csl:20]{\bibinfo{author}{\bibfnamefont{B.~P.} \bibnamefont{Pritchard}}, n.d.)} Especially, the computation of electron repulsion integrals (ERIs) requires a special attention. To overcome memory and disk problems associated with the conventional 4-index ERIs, tensor decomposition of ERIs has attracted significant  interest.\hyperref[csl:21]{(\bibinfo{author}{\bibfnamefont{J.~L.} \bibnamefont{Whitten}}, n.d.}; \hyperref[csl:22]{\bibinfo{author}{\bibfnamefont{B.~I.} \bibnamefont{Dunlap}}, n.d.}; \hyperref[csl:23]{\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Feyereisen}}, n.d.}; \hyperref[csl:24]{\bibinfo{author}{\bibfnamefont{O.}~\bibnamefont{Vahtras}}, n.d.}; \hyperref[csl:25]{\bibinfo{author}{\bibfnamefont{A.~P.} \bibnamefont{Rendell}} \bibnamefont{and}, n.d.}; \hyperref[csl:26]{\bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Weigend}}, n.d.}; \hyperref[csl:27]{\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Sodt}}, n.d.}; \hyperref[csl:28]{\bibinfo{author}{\bibfnamefont{H.-J.} \bibnamefont{Werner}}, n.d.}; \hyperref[csl:29]{\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Sch{\"u}tz}} \bibnamefont{and}, n.d.}; \hyperref[csl:30]{\bibinfo{author}{\bibfnamefont{H.-J.} \bibnamefont{Werner}} \bibnamefont{and}, n.d.}; \hyperref[csl:31]{\bibinfo{author}{\bibfnamefont{N.~H.~F.} \bibnamefont{Beebe}} \bibnamefont{and}, n.d.}; \hyperref[csl:32]{\bibinfo{author}{\bibfnamefont{I.}~\bibnamefont{Roeggen}} \bibnamefont{and}, n.d.}; \hyperref[csl:33]{\bibinfo{author}{\bibfnamefont{H.}~\bibnamefont{Koch}}, n.d.}; \hyperref[csl:34]{\bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Aquilante}}, n.d.}; \hyperref[csl:35]{\bibinfo{author}{\bibfnamefont{A.~E.} \bibnamefont{DePrince}} \bibnamefont{and}, n.d.}; \hyperref[csl:36]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 1179}; \hyperref[csl:37]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, n.d.}; \hyperref[csl:38]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 4415)} The most popular ERI factorization technique is the density fitting (DF) approach.\hyperref[csl:21]{(\bibinfo{author}{\bibfnamefont{J.~L.} \bibnamefont{Whitten}}, n.d.}; \hyperref[csl:22]{\bibinfo{author}{\bibfnamefont{B.~I.} \bibnamefont{Dunlap}}, n.d.}; \hyperref[csl:23]{\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Feyereisen}}, n.d.}; \hyperref[csl:24]{\bibinfo{author}{\bibfnamefont{O.}~\bibnamefont{Vahtras}}, n.d.}; \hyperref[csl:25]{\bibinfo{author}{\bibfnamefont{A.~P.} \bibnamefont{Rendell}} \bibnamefont{and}, n.d.}; \hyperref[csl:26]{\bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Weigend}}, n.d.}; \hyperref[csl:27]{\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Sodt}}, n.d.}; \hyperref[csl:28]{\bibinfo{author}{\bibfnamefont{H.-J.} \bibnamefont{Werner}}, n.d.}; \hyperref[csl:35]{\bibinfo{author}{\bibfnamefont{A.~E.} \bibnamefont{DePrince}} \bibnamefont{and}, n.d.}; \hyperref[csl:36]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 1179}; \hyperref[csl:37]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, n.d.}; \hyperref[csl:38]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 4415)} In the DF technique, the four-dimensional ERI tensor is replaced with the three-dimensional tensors, expanding ERIs in terms of a predefined auxiliary basis set.

In addition to the energy, analytic gradients for electron correlation methods with the DF approach have been reported for MP2 with conventional HF reference (DF-MP2),\hyperref[csl:39]{(\bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Weigend}} \bibnamefont{and}, n.d.}; \hyperref[csl:40]{\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{H\"{a}ttig}}, n.d.}; \hyperref[csl:41]{\bibinfo{author}{\bibfnamefont{R.~A.} \bibnamefont{Distasio}}, n.d.)} as well as with the DF-HF reference,\hyperref[csl:42]{(\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 1241)} the second-order coupled cluster (CC2),\hyperref[csl:43]{(\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{H\"{a}ttig}}, n.d.}; \hyperref[csl:44]{\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{K{\"o}hn}} \bibnamefont{and}, n.d.)}  local MP2 (DF-LMP2),\hyperref[csl:45]{(\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Sch{\"u}tz}}, n.d.)} complete active space second-order perturbation theory (DF-CASPT2),\hyperref[csl:46]{(\bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{Gy{\"o}rffy}}, n.d.)}  orbital-optimized MP2, MP2.5, and MP3 (DF-OMP2, DF-OMP2.5, DF-OMP3),\hyperref[csl:47]{(\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 4389}; \hyperref[csl:48]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, n.d.)} local time-dependent coupled cluster response theory (DF-TD-LCC2),\hyperref[csl:49]{(\bibinfo{author}{\bibfnamefont{L.}~\bibnamefont{Lederm{\"u}ller}}, n.d.)} the coupled-cluster singles and doubles and CCD methods (DF-CCSD and DF-CCD),\hyperref[csl:50]{(\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}} \bibnamefont{and}, n.d.)} and recently for the coupled-cluster singles and doubles with perturbative triples [DF-CCSD(T)] method.\hyperref[csl:51]{(\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}} \bibnamefont{and}, n.d.)}

One may use the existent 4-index integrals code to generate 3-index DF integrals. However, a better approach is using more specific algorithms for 3-index integrals. In a 2004 study,  Ahlrichs\hyperref[csl:52]{(\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Ahlrichs}}, n.d.)} presented an efficient algorithm to compute 3-center integrals based on the Obara-Saika (OS)  scheme. In 2017 and 2018 studies, \citeauthor{Kallay17:3index}\hyperref[csl:53]{(\bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Samu}} \bibnamefont{and}, n.d.}; \hyperref[csl:54]{\bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Samu}} \bibnamefont{and}, n.d.)} implemented various approaches, such as OS, McMurchie and Davidson(MD), and Rys quadrature, for computations of 3-index ERIs and their first derivatives. Efficient implementations of 3-index ERIs are available in {\sc Libint} library of Valeev,\hyperref[csl:55]{(\bibinfo{note}{E. F. Valeev, Libint: A Library For the Evaluation of Molecular 		Integrals of Many-Body Operators Over Gaussian Functions, 		\url{http://libint.valeyev.net/.} Last retrieved 13 August 2020.}, 2020}; \hyperref[csl:16]{\bibinfo{author}{\bibfnamefont{J.~P.} \bibnamefont{Kenny}}, n.d.}; \hyperref[csl:56]{\bibinfo{author}{\bibfnamefont{E.~F.} \bibnamefont{Valeev}} \bibnamefont{and}, n.d.)} which is available in {\sc Psi4}\hyperref[csl:57]{(\bibinfo{author}{\bibfnamefont{D.~G.~A.} \bibnamefont{Smith}}, n.d.)}, in {\sc Mrcc} program,\hyperref[csl:58]{(\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{K{\'{a}}llay}}, n.d.)} and in some commercial software.

In this research, an application programming interface (API) framework, denoted {\sc Molint} (MFW),  for the computation of molecular integrals and their first derivatives, over contracted Gaussian functions, for the density-fitted methods is reported. The MFW is an efficient and ``user-friendly" API framework, which provides a straightforward way to access molecular integrals and their derivatives. The MFW includes overlap, dipole, kinetic, potential, metric, and 3-index integrals as well as their first derivatives. Furthermore, the MFW provides a smooth approach to build the Fock matrix and compute the molecular gradients. The MFW is a {\sc C++/Fortran} hybrid code, which can take advantage of shared-memory parallel programming techniques. The features and efficiency of the MFW are demonstrated in illustrative applications.

\section*{\sffamily \Large Spherical Gaussians}
The general form of unnormalized primitive Gaussian type orbitals (GTOs) centered at $R_A$ can be written as

\begin{eqnarray}
G_{ijk}(a,{\bf |r - R_A|}) & = &  \big( x - A_x\big)^{i} \big( y - A_y\big)^{j} \big( z - A_z\big)^{k}   e^{ -a |r - R_A|^{2} } ,

\end{eqnarray}
which may be further simplified to:

\begin{eqnarray}
G_{ijk}( {\bf r }, a, {\bf A}) & = &  x_A^{i} y_A^{j} z_A^{k}   e^{ -a r_A^{2} } .

\end{eqnarray}

\begin{comment}
The full set of GTOs with the same angular momentum number $l = i + j + k$ and the same exponent $a$ are said to constitute a \emph{shell} . The number of GTOs in a shell with angular momentum $l$ is

\begin{eqnarray}
N_l^c& = &  \frac{(l+1)(l+2)}{2} .

\end{eqnarray}

\end{comment}

A real valued spherical GTOs with quantum numbers $l$ and $m$, with exponent $a$, centered at ${\bf A}$ is given by,

\begin{eqnarray}
G_{lm}( {\bf r }, a, {\bf A}) & = &  S_{lm}(x_A, y_A, z_A)   e^{ -a r_A^{2} } ,

\end{eqnarray}
where $S_{lm}({\bf r_A})$ is the real solid-harmonics.\hyperref[csl:1]{(\bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Helgaker}}, n.d.)}

Contracted Gaussian type orbitals (CGTOs) can be written as

\begin{eqnarray}
G_{\mu lm}( {\bf r }, a, {\bf A}) & = & \sum_{p}  d_{\mu p}  G_{lm}( {\bf r }, a_p, {\bf A}) ,

\end{eqnarray}
where $d_{\mu p} $ are the contraction coefficients.

\subsection*{Normalization of Contracted Cartesian GTOs}
The general form of unnormalized contracted Cartesian Gaussian type orbitals (CCGTOs) centered at the $R_A$ can be written as

\begin{eqnarray}
\chi_{\mu}({\bf r_A})  & = &  \sum_{p}^{n_p} d_{p \mu} G_{p}(a_{p \mu}, {\bf r_A}) ,

\end{eqnarray}
where $d_{p \mu} $ are the contraction coefficients and $n_p$ is the number of PGOs, and:

\begin{eqnarray}
G_{p}(a_{p \mu}, {\bf r_A }) & = &   x_A^{i}  y_A^{j}  z_A^{k}  e^{ -a_{p \mu}r_A^{2} } .

\end{eqnarray}

The contraction coefficients $d_{p \mu}$ already include normalization constants so that the resulting combination is properly normalized. Published contraction coefficients $\tilde{d}_{p \mu}$ are linear coefficients for normalized primitives, hence the normalization-including contraction coefficients $d_{p \mu}$ have to be computed from them. Hence, a normalized CGTO can be written as follows:

\begin{eqnarray}
\chi_{\mu}({\bf r_A})  & = &  N \sum_{p}^{n_p} \tilde{d}_{p \mu} \tilde{G}_{p} ,

\end{eqnarray}
where $\tilde{d}_{p \mu}$ are linear coefficients for normalized primitives $\tilde{G}_{p}$. Hence,

\begin{eqnarray}
\chi_{\mu}({\bf r_A})  & = &  N \sum_{p}^{n_p} \tilde{d}_{p \mu} N_p G_{p} .

\end{eqnarray}
Now let us define (apply primitive normalization):

\begin{eqnarray}
d_{p \mu}  & = &  N_p \tilde{d}_{p \mu} ,

\end{eqnarray}
then,

\begin{eqnarray}
\chi_{\mu}({\bf r_A})  & = &  N \sum_{p}^{n_p} d_{p \mu} G_{p} .

\end{eqnarray}

Finally, we employ normalization scheme of the {\sc Libint} package.\hyperref[csl:55]{(\bibinfo{note}{E. F. Valeev, Libint: A Library For the Evaluation of Molecular 		Integrals of Many-Body Operators Over Gaussian Functions, 		\url{http://libint.valeyev.net/.} Last retrieved 13 August 2020.}, 2020}; \hyperref[csl:16]{\bibinfo{author}{\bibfnamefont{J.~P.} \bibnamefont{Kenny}}, n.d.)} In this scheme, Cartesian Gaussian functions in a shell of angular momentum $l$ have the same normalization factor $N$, which is determined such that
the Cartesian functions $x^l$, $y^l$, and $z^l$ are normalized to unity.

\begin{eqnarray}
\pi^{3/2} \frac{ \big( 2l-1\big)!! }{ 2^l   } \sum_{p,q}^{n_p,n_q} \frac{d_{p \mu} d_{q \mu} }{\big(a_{p \mu} + b_{q \mu}\big)^{l+3/2}}  & = &  \frac{1}{N^2} ,

\end{eqnarray}
the final coefficients are:

\begin{eqnarray}
d^\prime_{p \mu} & = &  N d_{p \mu} .

\end{eqnarray}

\section*{\sffamily \Large Hermite Gaussians}
The Hermite Gaussian with exponent $p$ centered on ${\bf P}$ is defined as follows

\begin{eqnarray}
\Lambda_{tuv}({\bf r}, p, {\bf P}) & = &  \Big( \dpd{}{P_x} \Big)^{t} \Big( \dpd{}{P_y} \Big)^{u} \Big( \dpd{}{P_z} \Big)^{v}  e^{ -p r_P^{2} } ,

\end{eqnarray}
where $ {\bf r_P} = {\bf r} - {\bf P}$.

This functions are can be factorized as in the case of Cartesian GTOs

\begin{eqnarray}
\Lambda_{tuv}({\bf r}, p, {\bf P}) & = &  \Lambda_{t}(x, p, P_x) \Lambda_{u}(y, p, P_y) \Lambda_{v}(x, p, P_z) .

\end{eqnarray}
The $x$-component for example

\begin{eqnarray}
\Lambda_{t}(x, p, P_x) & = &  \Big( \dpd{}{P_x} \Big)^{t} \Lambda_{0}(x, p, P_x) ,

\end{eqnarray}
where

\begin{eqnarray}
\Lambda_{0}(x, p, P_x) & = &  e^{ -p x_P^{2} } .

\end{eqnarray}

Hermite functions have the following recursion relation:

\begin{eqnarray}
x_P \Lambda_{t}  & = &  \frac{1}{2p} \Lambda_{t+1} + t \Lambda_{t-1} .
\label{HG:rec2} 

\end{eqnarray}

\section*{\sffamily \Large Evaluation of The Boys Function}
Before continue to the evaluation of molecular integrals let us consider the Boys function, which is central for the computation of potential and electron repulsion integrals.

\begin{eqnarray}
F_{n}(x) & = & \int_{0}^{1} t^{2n} e^{ -x t^{2} }  \,dt ,

\end{eqnarray}
where $F_{n}(x)$ is the $n$-th order Boys function. Asymptotic values of the Boys function are

\begin{eqnarray}
F_{n}(0) & = & \int_{0}^{1} t^{2n}  \,dt \ = \ \frac{1}{2n+1} ,
\label{boys:xzero} 

\end{eqnarray}
and at large $x$:

\begin{eqnarray}
F_{n}(x) & \approx & \frac{(2n-1)!!}{2^{n+1}} \sqrt{\frac{\pi}{x^{2n+1}}} .
\label{boys:func:large} 

\end{eqnarray}

The Boys function satisfies the upward recursion:

\begin{eqnarray}
F_{n+1}(x) & = &\frac{(2n+1) F_{n}(x) - exp(-x)}{2x} ,

\end{eqnarray}
and the downward recursion:

\begin{eqnarray}
F_{n-1}(x) & = &\frac{2x F_{n}(x) + exp(-x)}{2n-1} .

\end{eqnarray}
The downward recursion is computationally more stable for small $x$ values.

Saunders\hyperref[csl:59]{(\bibinfo{author}{\bibfnamefont{V.~R.} \bibnamefont{Saunders}}, n.d.)} suggest the following series expansion for the boys function:

\begin{eqnarray}
F_{n}(x) & = & e^{-x} \sum_{k=0}^{\infty} \frac{(2x)^k}{(2n+1) (2n+3) \dots (2n+2k+1)}, 
\label{boys:func:saunders} 

\end{eqnarray}
which may be written as

\begin{eqnarray}
F_{n}(x) & = & e^{-x} \sum_{k=0}^{\infty} \frac{(2x)^k}{g_{nk} } ,
\label{boys:func:saunders2} 

\end{eqnarray}
where

\begin{eqnarray}
g_{nk} & = & \prod_{i=0}^{k} (2n+2i+1) .

\end{eqnarray}
Saunders note that this series is unconditionally convergent. However, for large values of $n$, convergence can be slow.

A faster convergent modified Maclaurin series expansion is suggested by Gill, Head-Gordon, and People\hyperref[csl:9]{(\bibinfo{author}{\bibfnamefont{P.~M.~W.} \bibnamefont{Gill}}, n.d.)} as follows:

\begin{eqnarray}
F_{n}(x) & = & \frac{1}{2} e^{-x} \frac{1}{n+1/2} \Bigg(  1 + \frac{x}{n+3/2} \Big( 1 + \frac{x^2}{n+5/2}(1+\dots) \Big) \Bigg) ,
\label{boys:func:ghp} 

\end{eqnarray}
The GHP method converges faster than Saunders' method.

For the small arguments, we can expand the Boys function as a Taylor series as follows,

\begin{eqnarray}
F_{n}(x) & = & \sum_{k=0}^{\infty} \frac{ F_{n+k}(x_0) (-\Delta x)^k}{k!} ,
\label{boys:func:taylor} 

\end{eqnarray}
where $\Delta x$ is generally chosen as 0.1 and $x_0 = x - \Delta x$. The first six terms of the Taylor series is a good approximation. Hence, we calculate $F_{n}(x_0), \dots , F_{n+k_{max}}(x_0)$ from the series expansion in Eq.\eqref{boys:func:saunders}--Eq.\eqref{boys:func:ghp} and we use these values to calculate $F_{n}(x)$. Further, we only need to compute the highest function $F_{n_{max}}(x)$ required and all other lower order functions can be calculated from the downward recursion relation

\begin{eqnarray}
F_{n}(x) & = &\frac{2x F_{n+1}(x) + exp(-x)}{2n+1} .
\label{boys:down} 

\end{eqnarray}

\paragraph{Final Strategy:}
Let us assume that we use a taylor expansion of order $k+1$ for the small argument.

\begin{enumerate}
\item Pretabulate $F_{n}(x)$ for $0<x\le30$ and $n=n_{max}+k$ with the $\Delta x = 0.1$ 
\item Obtain all grid values for $n=0-n_{max}+k-1$ from the downward recursion.
\item To compute $F_{n}(x)$ for an arbitrary $x$, find the closes point on the grid $x_g$, then use a  $k+1$-term Taylor expansion Eq.\eqref{boys:func:taylor}, where $\Delta x= x-x_g$ should be used.
\item For $x>30$ we may use the large argument.
\item For $x=30$ we may use the Eq.\eqref{boys:xzero}.

\end{enumerate}

\section*{\sffamily \Large McMurchie-Davidson (MD) Scheme}

\subsection*{Overlap Distributions from Hermite Gaussians}
Since the overlap distribution $\Omega_{ij}$ is a polynomial of degree $i+j$ in $x_p$, we can expand it as follows

\begin{eqnarray}
\Omega_{ij}(x) & = & \sum_{t=0}^{i+j} E_{t}^{ij} \Lambda_{t}(x_P) ,

\end{eqnarray}
where expansion coefficients are constant, they are independent of the electronic coordinates. The expansion coefficients satisfy the following relation

\begin{eqnarray}
E_{t}^{ij}  & = &  0, \quad t < 0 \quad or \quad t > i+j .

\end{eqnarray}

Further, the Hermite expansion coefficients obey the following recursion relations,

\begin{eqnarray}
E_{t}^{i+1,j}  & = & X_{PA} E_{t}^{ij} + \frac{1}{2p} \big( i E_{t}^{i-1,j} + j E_{t}^{i,j-1} + E_{t-1}^{ij}  \big) ,
\label{MD:E1} 

\end{eqnarray}

\begin{eqnarray}
E_{t}^{i,j+1}  & = & X_{PB} E_{t}^{ij} + \frac{1}{2p} \big( i E_{t}^{i-1,j} + j E_{t}^{i,j-1} + E_{t-1}^{ij}  \big) ,
\label{MD:E2} 

\end{eqnarray}
with

\begin{eqnarray}
E_{0}^{00}  & = & K_{AB}^{x} ,\\
E_{t}^{00}  & = & 0, \quad t > 0 .

\end{eqnarray}

\subsection*{One-Electron Integrals: Overlap, Kinetic, and Dipole Integrals}
The overlap, kinetic, and dipole integrals are defined as follows:

\begin{eqnarray}
S_{ab} & = & \langle G_{a} | G_{b} \rangle ,

\end{eqnarray}

\begin{eqnarray}
T_{ab} & = & -\frac{1}{2} \langle G_{a} | \nabla^{2} | G_{b} \rangle ,

\end{eqnarray}

\begin{eqnarray}
\mu_{ab} & = & \langle G_{a} | {\bf r}_C | G_{b} \rangle ,

\end{eqnarray}
where $S_{ab}$, $T_{ab}$, and $\mu_{ab}$ are overlap, kinetic, and dipole integrals, respectively, and ${\bf C}$ is the origin of the Cartesian dipole-moment. The evaluation of these integrals with the Hermite expansion is well-documented by Helgaker \emph{et al} .\hyperref[csl:1]{(\bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Helgaker}}, n.d.)} Hence, we will not repeat it.

\subsection*{Hermite Coulomb Integrals}
Hermite Coulomb integrals are required for the evaluation of the potential energy and electron repulsion integrals.

\begin{eqnarray}
R_{tuv}(p, R_{PC})  & = &  \Big( \dpd{}{P_x} \Big)^{t} \Big( \dpd{}{P_y} \Big)^{u} \Big( \dpd{}{P_z} \Big)^{v} F_{0}(pR_{PC}^2) ,

\end{eqnarray}
which can be obtained from the following auxiliary function

\begin{eqnarray}
R_{tuv}^{n}(p, R_{PC}) & = & (-2p)^n \Big( \dpd{}{P_x} \Big)^{t} \Big( \dpd{}{P_y} \Big)^{u} \Big( \dpd{}{P_z} \Big)^{v} F_{n}(pR_{PC}^2) .
\label{hci} 

\end{eqnarray}

Hermite coulomb integrals satisfy the following recursion relations :

\begin{eqnarray}
R_{t+1,u,v}^{n} & = &  t R_{t-1,u,v}^{n+1} + X_{PC} R_{tuv}^{n+1} ,
\label{hci:rec1} 

\end{eqnarray}

\begin{eqnarray}
R_{t,u+1,v}^{n} & = &  u R_{t,u-1,v}^{n+1} + Y_{PC} R_{tuv}^{n+1} ,
\label{hci:rec2} 

\end{eqnarray}

\begin{eqnarray}
R_{t,u,v+1}^{n} & = &  v R_{t,u,v-1}^{n+1} + Z_{PC} R_{tuv}^{n+1} ,
\label{hci:rec3} 

\end{eqnarray}
Starting from the highest

\begin{eqnarray}
R_{000}^{n_{max}}(p, R_{PC})  & = & (-2p)^{n_{max}} F_{n_{max}}(pR_{PC}^2) , 
\label{hci:max} 

\end{eqnarray}
and using the recursion relations of Eqs.\eqref{hci:rec1}--\eqref{hci:rec3} we can generate all Hermite Coulomb integrals. The maximum value of $n$ is

\begin{eqnarray}
n_{max}  & = & t + u + v .

\end{eqnarray}

\subsection*{Potential Energy Integrals}
The potential energy integral can be written as

\begin{eqnarray}
V_{ab} & = & -\sum_{C} \langle G_{a} | \frac{Z_C}{r_C} | G_{b} \rangle ,

\end{eqnarray}
which can be written explicitly as follows:

\begin{eqnarray}
V_{ab} & = & -\frac{2\pi}{p} \sum_{C}  \sum_{tuv} Z_C E_{tuv}^{ab}  R_{tuv}(p, R_{PC}) .

\end{eqnarray}
Now we may write,

\begin{eqnarray}
V_{ab} & = & \frac{2\pi}{p} \sum_{tuv} E_{tuv}^{ab}  \Big[-\sum_{C}  Z_C R_{tuv}(p, R_{PC}) \Big] ,

\end{eqnarray}
let us define

\begin{eqnarray}
R^\prime_{tuv}  & = & -\sum_{C}  Z_C R_{tuv}(p, R_{PC}) ,

\end{eqnarray}
hence, our final formula becomes:

\begin{eqnarray}
V_{ab} & = & \frac{2\pi}{p} \sum_{tuv} E_{tuv}^{ab} R^\prime_{tuv} .

\end{eqnarray}
Therefore, the contributions from all nuclei can be computed before the transformation.

\subsection*{Hermite Expansion for GTOs}
Before proceed to the metric and 3-index ERIs, let us consider the Hermite expansion for a GTO. For the $x$--component it can be written as:

\begin{eqnarray}
G_{i} & = &  \sum_{t=0}^{i} E_{t}^{i} \Lambda_{t}(x_A) ,

\end{eqnarray}

Starting from the following Using derivative relations of GTO we obtain the following recursion relations:

\begin{eqnarray}
E_{t}^{i+1}  & = & (t+1) E_{t+1}^{i}  + \frac{1}{2a} E_{t-1}^{i}  ,

\end{eqnarray}
and

\begin{eqnarray}
E_{t}^{i+1}  & = & \frac{1}{2a} \big( i E_{t}^{i-1} + E_{t-1}^{i} \big) .

\end{eqnarray}

\subsection*{Metric Integrals}
The metric integral is defined as follows,

\begin{eqnarray}
J_{PQ} & = & \int\int  \varphi_{P}({\bf r_{1}}) \frac{1}{r_{12}} \varphi_{Q}({\bf r_{2}})  \,d{\bf r_{1}} d{\bf r_{2}} ,

\end{eqnarray}
which may be written as (in terms of primitives):

\begin{eqnarray}
(a|b)  & = &  \int \int \frac{G_{a}({\bf r}_1) G_{b}({\bf r}_2)}{r_{12}}  \,d{\bf r}_1 d{\bf r}_2 ,

\end{eqnarray}
Using Hermite expansion and the Laplace transform of $r_{12}$ we obtain the following final expression:

\begin{eqnarray}
(a|b)  & = &  \frac{2 \pi^{5/2}}{ab \sqrt{a+b}} \sum_{tuv} E_{tuv}^{a} 
\sum_{t^\prime u^\prime v^\prime} (-1)^{t^\prime + u^\prime + v^\prime}  E_{t^\prime u^\prime v^\prime}^{b} R_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{AB}) ,
\label{peri:J} 

\end{eqnarray}
where

\begin{eqnarray}
\alpha & = &  \frac{ab}{a+b} .

\end{eqnarray}

\subsubsection*{The MD4 Algorithm for Metric Integrals}
Let us define

\begin{eqnarray}
F_{t^\prime u^\prime v^\prime}^{b} & = &   (-1)^{t^\prime + u^\prime + v^\prime} E_{t^\prime u^\prime v^\prime}^{b},

\end{eqnarray}
and

\begin{eqnarray}
{\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{AB})  & = &  \frac{2 \pi^{5/2}}{ab \sqrt{a+b}} R_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{AB})  .

\end{eqnarray}
Then, we may write

\begin{eqnarray}
(a|b)  & = &  \sum_{tuv} E_{tuv}^{a} \sum_{t^\prime u^\prime v^\prime} F_{t^\prime u^\prime v^\prime}^{b} {\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{AB}) .

\end{eqnarray}

For efficiency we evaluate above equation in two steps as follows

\begin{eqnarray}
(a|b)  & = &  \sum_{tuv} E_{tuv}^{a} g_{tuv}^{b}  ,

\end{eqnarray}
where

\begin{eqnarray}
g_{tuv}^{b}  & = & \sum_{t^\prime u^\prime v^\prime} F_{t^\prime u^\prime v^\prime}^{b} {\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{AB}) .

\end{eqnarray}

Further efficiency can be achieved by decomposing ERI and $g_{tuv}^{b} $ into Cartesian components. Let us start with $g_{tuv}^{b} $

\begin{eqnarray}
g_{tuv}^{b}  & = & g_{tuv}^{i^\prime k^\prime m^\prime} ,

\end{eqnarray}
and

\begin{eqnarray}
g_{tuv}^{i^\prime k^\prime m^\prime }  & = &
\sum_{t^\prime} F_{t^\prime}^{i^\prime}  \sum_{u^\prime} F_{u^\prime}^{k^\prime} 
\sum_{v^\prime} F_{v^\prime}^{m^\prime} {\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{AB}) .

\end{eqnarray}

Then, we can write

\begin{eqnarray}
(a|b)  & = &  g_{ikm}^{i^\prime k^\prime m^\prime } ,

\end{eqnarray}

\begin{eqnarray}
g_{ikm}^{i^\prime k^\prime m^\prime }  & = &  
\sum_{t} E_{t}^{i} \sum_{u} E_{u}^{k}  \sum_{v} E_{v}^{m} g_{tuv}^{i^\prime k^\prime m^\prime} .

\end{eqnarray}
Following Helgaker et al,\hyperref[csl:1]{(\bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Helgaker}}, n.d.)} we call this algorithm as the MD4 algorithm.

\subsection*{3-Index Integrals}
3-index two-electron integral (TEI) can be written as

\begin{eqnarray}
(\mu \nu | Q) & = & \int\int \chi_{\mu}({\bf r_{1}}) \chi_{\nu}({\bf r_{1}}) \frac{1}{r_{12}} \varphi_{Q}({\bf r_{2}})  \,d{\bf r_{1}} d{\bf r_{2}} ,

\end{eqnarray}

3-index TEI may be written as (in terms of primitives):

\begin{eqnarray}
(ab|C)  & = &  \int \int \frac{\Omega_{ab}({\bf r}_1) G_{C}({\bf r}_2)}{r_{12}}  \,d{\bf r}_1 d{\bf r}_2,

\end{eqnarray}
where $G_{C} $ is a primitive function of auxiliary basis set.

Applying the same machinery of the MD method we obtain:

\begin{eqnarray}
(ab|C)  & = &  \frac{2 \pi^{5/2}}{pq \sqrt{p+q}} \sum_{tuv} E_{tuv}^{ab} 
\sum_{t^\prime u^\prime v^\prime} (-1)^{t^\prime + u^\prime + v^\prime}  E_{t^\prime u^\prime v^\prime}^{C} R_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{PC}) ,
\label{peri:3index} 

\end{eqnarray}
where
$c$ is the exponent of $G_{C}$ and

\begin{eqnarray}
P_{x1} & = & \frac{ a A_{x1} + b B_{x1} }{p} ,

\end{eqnarray}

\begin{eqnarray}
p  & = &  a + b  ,

\end{eqnarray}

\begin{eqnarray}
\alpha & = &  \frac{pc}{p+c} .

\end{eqnarray}

\subsubsection*{The MD4 Algorithm for 3-Index TEI}
Let us define

\begin{eqnarray}
F_{t^\prime u^\prime v^\prime}^{C} & = &   (-1)^{t^\prime + u^\prime + v^\prime} E_{t^\prime u^\prime v^\prime}^{C} ,

\end{eqnarray}
and

\begin{eqnarray}
{\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{PC})  & = &  \frac{2 \pi^{5/2}}{pc \sqrt{p+c}} R_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{PC}) ,

\end{eqnarray}
Then, we may write

\begin{eqnarray}
(ab|C)  & = &  \sum_{tuv} E_{tuv}^{ab} \sum_{t^\prime u^\prime v^\prime} F_{t^\prime u^\prime v^\prime}^{C} {\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{PC}) .

\end{eqnarray}

For efficiency we evaluate above equation in two steps as follows

\begin{eqnarray}
(ab|C)  & = &  \sum_{tuv} E_{tuv}^{ab} g_{tuv}^{C}  ,

\end{eqnarray}
where

\begin{eqnarray}
g_{tuv}^{C}  & = & \sum_{t^\prime u^\prime v^\prime} F_{t^\prime u^\prime v^\prime}^{Q} {\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{PC}) .

\end{eqnarray}

Further efficiency can be achieved by decomposing ERI and $g_{tuv}^{C} $ into Cartesian components. Let us start with $g_{tuv}^{C} $

\begin{eqnarray}
g_{tuv}^{C}  & = & g_{tuv}^{i^\prime k^\prime m^\prime } ,

\end{eqnarray}
and

\begin{eqnarray}
g_{tuv}^{i^\prime k^\prime m^\prime }   & = &
\sum_{t^\prime} F_{t^\prime}^{i^\prime}  \sum_{u^\prime} F_{u^\prime}^{k^\prime} 
\sum_{v^\prime} F_{v^\prime}^{m^\prime} {\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{PC}) .

\end{eqnarray}

Then, we can write

\begin{eqnarray}
(ab|C)  & = &  g_{ijklmn}^{i^\prime k^\prime m^\prime } ,

\end{eqnarray}

\begin{eqnarray}
g_{ijklmn}^{i^\prime k^\prime m^\prime }   & = &  
\sum_{t} E_{t}^{ij} \sum_{u} E_{u}^{kl}  \sum_{v} E_{v}^{mn} g_{tuv}^{i^\prime k^\prime m^\prime } .

\end{eqnarray}

\section*{\sffamily \Large Obra-Saika (OS) Scheme}

\subsection*{3-Index Integrals}
Let us remember

\begin{eqnarray}
(ab|C)  & = &  \sum_{tuv} E_{tuv}^{ab} 
\sum_{t^\prime u^\prime v^\prime} (-1)^{t^\prime + u^\prime + v^\prime}  E_{t^\prime u^\prime v^\prime}^{C} {\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{PC}) ,
\label{peri:3index2} 

\end{eqnarray}
where $p=a+b$ and $c$ is the exponent of $G_{C}$ and

\begin{eqnarray}
\alpha & = &  \frac{pc}{p+c} ,

\end{eqnarray}
and

\begin{eqnarray}
{\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{PC}) & = &  \frac{2 \pi^{5/2}}{pc \sqrt{p+c}} R_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{PC}) .

\end{eqnarray}

Now, let us introduce an auxiliary integral as follows:

\begin{eqnarray}
\Theta_{i_x j_x k_x; i_y j_y k_y; i_z j_z k_z}^{N} & = &  (-2\alpha)^{-N} \sum_{tuv} E_{t}^{i_x j_x} E_{u}^{i_y j_y} E_{v}^{i_z j_z} 
\nonumber\\
& \times & 
\sum_{t^\prime u^\prime v^\prime} (-1)^{t^\prime + u^\prime + v^\prime}  E_{t^\prime}^{k_x}  E_{u^\prime}^{k_y}  E_{v^\prime}^{k_z}  
{\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}^{N}(\alpha, R_{PC}) ,
\nonumber\\
\label{3index:aux} 

\end{eqnarray}
note that

\begin{eqnarray}
\Theta_{i_x j_x k_x; i_y j_y k_y i_z j_z k_z}^{0}  & = &  g_{i_x j_x k_x; i_y j_y k_y; i_z j_z k_z} .

\end{eqnarray}

Now let us define

\begin{eqnarray}
{\cal F}_{N}(\alpha, R_{PC})  & = & \frac{2 \pi^{5/2}}{pc \sqrt{p+c}} F_{N}(\alpha R_{PC}^2) ,

\end{eqnarray}
and  introduce an auxiliary function as follows:

\begin{eqnarray}
\Theta_{ijk}^{N} & = &  (-2\alpha)^{-N} \sum_{tuv} \sum_{t^\prime u^\prime v^\prime} (-1)^{t^\prime + u^\prime + v^\prime}  E_{t}^{ij} E_{t^\prime}^{k}  E_{uv u^\prime v^\prime} 
{\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}^{N}(\alpha, R_{PC}) ,
\label{3index:aux2} 

\end{eqnarray}
where

\begin{eqnarray}
E_{uvu^\prime v^\prime} & = & E_{u}^{i_y j_y} E_{v}^{i_z j_z} E_{u^\prime}^{k_y}  E_{v^\prime}^{k_z}  ,

\end{eqnarray}
and $N \le l_1 + l_2 + l_3$ .

\subsubsection*{The OS2 Algorithm for 3-Index TEI}
In this section let us use an abstract notation for the auxiliary integral as follows:

\begin{eqnarray}
\Theta_{abC}^{N} & = &  (-2\alpha)^{-N} \sum_{tuv} E_{tuv}^{ab} 
\sum_{t^\prime u^\prime v^\prime} (-1)^{t^\prime + u^\prime + v^\prime}  E_{t^\prime u^\prime v^\prime}^{C}  
{\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}^{N}(\alpha, R_{PC}) ,
\label{3index:aux:abs} 

\end{eqnarray}
note that

\begin{eqnarray}
\Theta_{abC}^{0}  & = &  (ab|C) .

\end{eqnarray}
In this notation angular momentum components of a $G_a$ function are $i_a,j_a,k_a$ for $x,y,z,$ directions, respectively.

We may start with the following auxiliary integral:

\begin{eqnarray}
\Theta_{000}^{N} & = &  K_{AB}^{x} K_{AB}^{y} K_{AB}^{z}  {\cal F}_{N}(\alpha R_{PC}^2)  ,
\label{peri:3index:aux:abs} 

\end{eqnarray}
where $N \le l_1 + l_2 + l_3$

Then, we use Ahlrich's OS2 algorithm.\hyperref[csl:52]{(\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Ahlrichs}}, n.d.)} In the OS2 algorithm we start with Eq.\eqref{peri:3index:aux:abs} and apply vertical recurrence relation (VRR)  to increment the angular momentum of the first function on the bra side (given here for the $x$ direction)
as:

\begin{eqnarray}
\Theta_{a+1_x,0,0}^{N} & = &  X_{PA} \Theta_{a00}^{N}  \ - \ \frac{\alpha}{p} X_{PC} \Theta_{a00}^{N+1}
\ + \ \frac{i_a}{2p} \Big( \Theta_{a-1_x,0,0}^{N} -  \frac{\alpha}{p} \Theta_{a-1_x,0,0}^{N+1}\Big),
\label{3index:vrr:bra:v2} 

\end{eqnarray}
where $i_a \le l_1 + l_2 $.

When all $\Theta_{a+b,0,0}^{N}$ functions are computed, then we proceed to the ket side as the following VRR:

\begin{eqnarray}
\Theta_{a,0,C +1_x}^{N} & = &  \frac{\alpha}{c} X_{PC} \Theta_{a,0,C}^{N+1} \ + \ \frac{i_a \alpha}{2pc} \Theta_{a-1_x,0,C}^{N+1} 
\ + \ \frac{i_C}{2c} \Big( \Theta_{a,0,C-1_x}^{N} -  \frac{\alpha}{c} \Theta_{a,0,C-1_x}^{N+1}\Big) ,
\label{3index:vrr:ket:v2} 

\end{eqnarray}
if auxiliary function is a solid harmonic then,

\begin{eqnarray}
\Theta_{a,0, C +1_x}^{N} & = &  \frac{\alpha}{c} X_{PC} \Theta_{a,0,C}^{N+1} \ + \ \frac{i_a}{2(p+c)} \Theta_{a-1_x,0,C}^{N+1} ,
\label{3index:vrr:ket:sh:v2} 

\end{eqnarray}
where $i_C \le l_3$.

Now we have all $\Theta_{a+b,0,C}^{N}$ functions. Hence, we can obtain all $(a+b,0| C)$ type integrals from $\Theta_{a+b,0,C}^{0}$. The next step is \emph{horizontal recursion relation} (HRR) in the contracted basis:\hyperref[csl:8]{(\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Head-Gordon}} \bibnamefont{and}, n.d.)}

\begin{eqnarray}
(a, b +1_x | C) & = &  (a +1_x, b | C)   \ + \ R_{AB} (ab | C) .

\end{eqnarray}

\subsubsection*{The OS1 Algorithm for 3-Index TEI}
In the OS1 algorithm we start with Eq.\eqref{peri:3index:aux:abs} with $N \le l_1 + l_2 + l_3$ and apply VRR to increment the angular momentum of the first function on the bra side (given here for the $x$ direction)
as:

\begin{eqnarray}
\Theta_{a+1_x,0,0}^{N} & = &  X_{PA} \Theta_{a00}^{N}  \ - \ \frac{\alpha}{p} X_{PC} \Theta_{a00}^{N+1}
\ + \ \frac{i_a}{2p} \Big( \Theta_{a-1_x,0,0}^{N} -  \frac{\alpha}{p} \Theta_{a-1_x,0,0}^{N+1}\Big) ,
\label{3index:vrr:bra:v3} 

\end{eqnarray}
where $i_a \le l_1 + l_2 + l_3$.

When all $\Theta_{a+b+C,0,0}^{N}$ functions are computed, then we proceed to the ket side as the following the \emph{electron transfer recursion relation} ,\hyperref[csl:10]{(\bibinfo{author}{\bibfnamefont{T.~P.} \bibnamefont{Hamilton}} \bibnamefont{and}, n.d.)} ETRR:

\begin{eqnarray}
\Theta_{a,0,C+1_x}^{N} & = & -\frac{b}{c} X_{AB} \Theta_{a0C}^{N}  
\ + \ \frac{i_a}{2c} \Theta_{a-1_x,0,C}^{N} \ + \ \frac{i_C}{2c} \Theta_{a,0,C-1_x}^{N}
\ - \ \frac{p}{c} \Theta_{a+1_x,0,C}^{N} ,
\label{3index:trr} 

\end{eqnarray}
if auxiliary function is a solid harmonic then,

\begin{eqnarray}
\Theta_{a,0,C+1_x}^{N} & = & -\frac{b}{c} X_{AB} \Theta_{a0C}^{N}  
\ + \ \frac{i_a}{2c} \Theta_{a-1_x,0,C}^{N} 
\ - \ \frac{p}{c} \Theta_{a+1_x,0,C}^{N} .
\label{3index:trr2} 

\end{eqnarray}
Further, since ETRR has the same order in both side of equations, it is enough to use ($N=0$) in the ETRR step:

\begin{eqnarray}
(a,0| C+1_x)& = & -\frac{b}{c} X_{AB} (a,0| C)
\ + \ \frac{i_a}{2c} (a-1_x,0| C)
\ - \ \frac{p}{c} (a+1_x,0| C) ,
\label{3index:trr3} 

\end{eqnarray}
It is also possible to setup a reverse ETRR as follows:

\begin{eqnarray}
(a+1_x,0| C) & = & -\frac{b}{p} X_{AB} (a,0| C)
\ + \ \frac{i_a}{2p} (a-1_x,0| C)
\ - \ \frac{c}{p} (a,0| C+1_x) .
\label{3index:trr:rev} 

\end{eqnarray}
Now we have all $(a+b,0| C)$ type integrals. The next step is the HRR as in the case of OS2.

\section*{\sffamily \Large First Derivatives of Integrals}
The first derivative relation for a GTO can be written as

\begin{eqnarray}
\dpd{G_{i}}{A_x} & = &   2a G_{i+1} - i G_{i-1} ,
\label{gto:grad} 

\end{eqnarray}
To proceed higher derivative let us introduce the notation

\begin{eqnarray}
\dpd[q]{G_{i}}{A_x}  & = &  G_{i}^q ,

\end{eqnarray}
Then we have

\begin{eqnarray}
G_{i}^{q+1} & = &  2a G_{i+1}^q - i G_{i-1}^q  .

\end{eqnarray}

Further, it is also interesting to consider the first derivatives of the overlap distribution. Differentiating the overlap distribution we get the following equations:

\begin{eqnarray}
\dpd{\Omega_{ij}^{x}}{A_x} & = &   2a \Omega_{i+1,j}^{x} - i \Omega_{i-1,j}^{x} ,
\label{overlap:dist:grad:A} 

\end{eqnarray}

\begin{eqnarray}
\dpd{\Omega_{ij}^{x}}{B_x} & = &   2b \Omega_{i,j+1}^{x} - j \Omega_{i,j-1}^{x} .
\label{overlap:dist:grad:B} 

\end{eqnarray}
Further, we may expand the derivative relations using the Hermite functions as follows:\hyperref[csl:60]{(\bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Helgaker}} \bibnamefont{and}, n.d.)}

\begin{eqnarray}
\dpd{\Omega_{ij}^{x}}{A_x} & = &   \sum_{t=0}^{i+j+1} \Big( 2a E_{t}^{i+1,j} - i E_{t}^{i-1,j} \Big) \Lambda_{t} ,
\label{overlap:dist:grad:A3} 

\end{eqnarray}

\begin{eqnarray}
\dpd{\Omega_{ij}^{x}}{B_x} & = &   \sum_{t=0}^{i+j+1} \Big( 2b E_{t}^{i,j+1} - j E_{t}^{i,j-1} \Big) \Lambda_{t} .
\label{overlap:dist:grad:B3} 

\end{eqnarray}

\section*{\sffamily \Large The MD Scheme for the First Derivatives}

\subsection*{Overlap Derivatives}
At first note that we have a \emph{translational invarince} relation (TIR):

\begin{eqnarray}
\frac{\partial S_{ab} }{\partial A_x} \ + \ \frac{\partial S_{ab} }{\partial B_x} & = &  0 ,

\end{eqnarray}
where $A$ and $B$ are the the centers of the first and the second GTO, respectively.

Derivatives of overlap with respect to the $x$ coordinate of the nucleus $N$ can be written as:

\begin{eqnarray}
\frac{\partial S_{ab}}{\partial N_x} & = &  \frac{\partial S_{ab}}{\partial A_x} \delta_{AN} \ + \ \frac{\partial S_{ab}}{\partial B_x} \delta_{BN} ,

\end{eqnarray}
using TIR we may write:

\begin{eqnarray}
\frac{\partial S_{ab}}{\partial N_x} & = &  \frac{\partial S_{ab}}{\partial A_x} \ \big(\delta_{AN} - \delta_{BN}\big) .

\end{eqnarray}
Hence, we just need to differentiate with respect to one of centers only.

Using derivative properties of overlap distribution we obtain the following final formula:

\begin{eqnarray}
\frac{\partial S_{ab}}{\partial A_x} & = &  \big(2a E_{0}^{i+1,j} - i E_{0}^{i-1,j}\big) E_{0}^{kl} E_{0}^{mn} \big( \frac{\pi}{p} \big)^{3/2} ,

\end{eqnarray}

\subsection*{Dipole Derivatives}
Let us express $x$-component of $\mu_{ab}^x$ more explicitly:

\begin{eqnarray}
\mu_{ab}^x  & = &  \langle G_{a} | {\bf x}_C | G_{b} \rangle ,

\end{eqnarray}
where ${\bf C}$ is the origin of the Cartesian dipole-moment, typically zero.

For the dipole derivatives, the translational invariance relation is:

\begin{eqnarray}
\frac{\partial \mu_{ab} }{\partial A_x} \ + \ \frac{\partial \mu_{ab} }{\partial B_x} \ + \ \frac{\partial \mu_{ab} }{\partial C_x} & = &  0 ,

\end{eqnarray}
The last term is equal:

\begin{eqnarray}
\frac{\partial \mu_{ab} }{\partial C_x} & = &  -S_{ab} ,

\end{eqnarray}
hence,

\begin{eqnarray}
\frac{\partial \mu_{ab} }{\partial A_x} \ + \ \frac{\partial \mu_{ab} }{\partial B_x} \ - \ S_{ab} & = &  0 ,

\end{eqnarray}
If the dipole center is zero, then:

\begin{eqnarray}
\frac{\partial \mu_{ab} }{\partial A_x} \ + \ \frac{\partial \mu_{ab}}{\partial B_x} & = &  0 .

\end{eqnarray}

Derivatives of overlap with respect to the $x$ coordinate of the nucleus $N$ can be written as:

\begin{eqnarray}
\frac{\partial \mu_{ab}}{\partial N_x} & = &  \frac{\partial \mu_{ab} }{\partial A_x} \delta_{AN} \ + \ \frac{\partial \mu_{ab} }{\partial B_x} \delta_{BN} 
\ + \ \frac{\partial \mu_{ab} }{\partial C_x} \delta_{CN} , 

\end{eqnarray}
If the dipole origin is zero, then:

\begin{eqnarray}
\frac{\partial \mu_{ab}}{\partial N_x} & = &  \frac{\partial \mu_{ab} }{\partial A_x} \big(\delta_{AN} - \delta_{BN} \big) .

\end{eqnarray}

Using derivative relation of primitive Gaussians, we obtain derivatives for the $x$-component of the dipole integrals as follows:

\begin{eqnarray}
\frac{\partial \mu_{ab}^x }{\partial A_x} & = & \big(2a \mu_{i+1,j} - i \mu_{i-1,j}\big) S_{kl} S_{mn} ,

\end{eqnarray}

Similarly, the derivatives for the $y$-component of the dipole integrals can be written as follows:

\begin{eqnarray}
\frac{\partial \mu_{ab}^y }{\partial A_x} & = & \big(2a S_{i+1,j} - i S_{i-1,j}\big)  \mu_{kl} S_{mn} ,

\end{eqnarray}

Finally, the derivatives for the $z$-component of the dipole integrals are:

\begin{eqnarray}
\frac{\partial \mu_{ab}^z }{\partial A_x} & = & \big(2a S_{i+1,j} - i S_{i-1,j}\big) S_{kl} \mu_{mn} ,

\end{eqnarray}

\subsection*{Kinetic Energy Derivatives }
Using the TIR, derivatives of kinetic energy integrals with respect to the $x$ coordinate of the nucleus $N$ can be written as:

\begin{eqnarray}
\frac{\partial T_{ab}}{\partial N_x} & = &  \frac{\partial T_{ab}}{\partial A_x} \ \big(\delta_{AN} - \delta_{BN}\big) ,

\end{eqnarray}
Using derivative properties of overlap distribution we obtain the following final formula for the derivative with respect to $A_x$:

\begin{eqnarray}
\frac{\partial T_{ab}}{\partial A_x}  & = &  \big(2a T_{i+1,j} - i T_{i-1,j} \big)  S_{kl} S_{mn} 
\ + \ \big(2a S_{i+1,j} - i S_{i-1,j}  \big) \big(T_{kl} S_{mn} + S_{kl} T_{mn} \big) .

\end{eqnarray}
Explicit equations for overlap ($S_{ij}$) and kinetic integrals ($T_{ij}$) were given in Eq.\eqref{overlap} and Eq.\eqref{kinetic}. One may write the derivatives with respect to $A_y$ and $A_z$ in a similar way.

\subsection*{Potential Energy Derivatives}
Before starting to derivatives, let us define a more general potential function:

\begin{eqnarray}
V_{ab,C}^{qrs} & = & \langle G_{a} | \Big( \dpd{}{C_x} \Big)^{q} \Big( \dpd{}{C_y} \Big)^{r}  \Big( \dpd{}{C_z} \Big)^{s} \frac{1}{r_C} | G_{b} \rangle ,

\end{eqnarray}
It becomes:

\begin{eqnarray}
V_{ab,C}^{qrs}  & = & (-1)^{q+r+s} \ \frac{2\pi}{p} \sum_{tuv} E_{tuv}^{ab}  R_{t+q,u+r,v+s}(p, R_{PC})  ,

\end{eqnarray}

Hence, derivatives with respect to the nucleus $N$ can be written as:

\begin{eqnarray}
\frac{\partial V_{ab}}{\partial N_x} & = &  \frac{\partial V_{ab}}{\partial A_x} \ \delta_{AN} \ + \ \frac{\partial V_{ab}}{\partial B_x} \ \delta_{BN} \ - \ \bigg[ \sum_{C} Z_C V_{ab,C}^{100}  \ \delta_{CN}\bigg] ,

\end{eqnarray}
where

\begin{eqnarray}
V_{ab,C}^{100}  & = & - \ \frac{2\pi}{p} \sum_{tuv} E_{tuv}^{ab}  R_{t+1,u,v}(p, R_{PC}) ,

\end{eqnarray}

The TIR for potential energy derivatives can be written as follows:\hyperref[csl:61]{(\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Komornicki}}, n.d.)}

\begin{eqnarray}
\frac{\partial V_{ab,C}}{\partial A_x} \ + \ \frac{\partial V_{ab,C}}{\partial B_x} \ + \ V_{ab,C}^{100}  & = &  0 ,

\end{eqnarray}

Therefore, the derivative with respect to $N_x$ can be written as follows:

\begin{eqnarray}
\frac{\partial V_{ab}}{\partial N_x} & = &  \frac{\partial V_{ab}}{\partial A_x} \ \delta_{AN} \ + \ \frac{\partial V_{ab}}{\partial B_x} \ \delta_{BN} 
\ + \ Z_N \Big( \frac{\partial V_{ab,N}}{\partial A_x}  \ + \ \frac{\partial V_{ab,N}}{\partial B_x} \Big) .

\end{eqnarray}
Hence, the explicit equations become:

\begin{eqnarray}
\frac{\partial V_{ab}}{\partial A_x} & = &   \frac{2\pi}{p} \sum_{tuv} \Big( 2a E_{t}^{i+1,j} - i E_{t}^{i-1,j} \Big) E_{u}^{kl}  E_{v}^{mn} R^\prime_{tuv} ,

\end{eqnarray}
and

\begin{eqnarray}
\frac{\partial V_{ab}}{\partial B_x} & = &   \frac{2\pi}{p} \sum_{tuv} \Big( 2b E_{t}^{i,j+1} - j E_{t}^{i,j-1} \Big) E_{u}^{kl}  E_{v}^{mn} R^\prime_{tuv} .

\end{eqnarray}

\subsection*{Metric Derivatives}
At first note that we have the following TIR for metric integrals:

\begin{eqnarray}
\Big(\frac{\partial G_a }{\partial A_x} | \frac{1}{r_{12}}| G_b \Big) \ + \ \Big( G_a | \frac{1}{r_{12}} | \frac{\partial G_b }{\partial B_x} \Big) & = &  0 ,
\label{eq:metric:trans} 

\end{eqnarray}
Hence, if both orbitals have the same center, then the first derivative is zero. The derivative with respect to the center $N$ may be written as follows:

\begin{eqnarray}
\frac{\partial (G_a|G_b) }{\partial N_x} & = & \Big(\frac{\partial G_a }{\partial A_x} | G_b \Big) \ \big(\delta_{AN} - \delta_{BN} \big) .

\end{eqnarray}

Using derivative relation of primitive Gaussians, we obtain derivatives along the $x$-axis as follows:

\begin{eqnarray}
\frac{\partial (a|b) }{\partial A_x} & = &  \sum_{tuv} \big(2a E_{t}^{i+1} - i E_{t}^{i-1} \big) E_{u}^{k}  E_{v}^{m} 
\sum_{t^\prime u^\prime v^\prime}   F_{t^\prime u^\prime v^\prime}^{b} {\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{AB}) ,

\end{eqnarray}

\subsection*{3-Index Derivatives}
At first note that we have the following TIR for 3-index integrals:

\begin{eqnarray}
\frac{\partial (ab|C) }{\partial A_x} \ + \ \frac{\partial (ab|C) }{\partial B_x} \ + \ \frac{\partial (ab|C) }{\partial C_x}  & = &  0 ,
\label{eq:3index:trans} 

\end{eqnarray}
Hence, if all orbitals have the same center, then the first derivative is zero. The derivative with respect to the center $N$ may be written as follows:

\begin{eqnarray}
\frac{\partial (ab|C) }{\partial N_x} & = & \frac{\partial (ab|C) }{\partial A_x} \ \big(\delta_{AN} - \delta_{CN} \big) \ + \ \frac{\partial (ab|C) }{\partial B_x} \ \big(\delta_{BN} - \delta_{CN} \big) .

\end{eqnarray}
Therefore, we need only derivatives with respect to the centers $A$ and $B$.

Using derivative relation of primitive Gaussians, we obtain derivatives along the $x$-axis as follows:

\begin{eqnarray}
\frac{\partial (ab|C)  }{\partial A_x} & = &  
\sum_{tuv} \big(2a E_{t}^{i+1,j} - i E_{t}^{i-1,j} \big) E_{u}^{kl}  E_{v}^{mn} 
\nonumber\\
& \times & \sum_{t^\prime u^\prime v^\prime} F_{t^\prime u^\prime v^\prime}^{C} {\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{PC}) ,

\end{eqnarray}

\begin{eqnarray}
\frac{\partial (ab|C)  }{\partial B_x} & = &  
\sum_{tuv} \big(2b E_{t}^{i,j+1} - j E_{t}^{i,j-1} \big) E_{u}^{kl}  E_{v}^{mn} 
\nonumber\\
& \times & \sum_{t^\prime u^\prime v^\prime} F_{t^\prime u^\prime v^\prime}^{C} {\cal R}_{t + t^\prime, u + u^\prime, v + v^\prime}(\alpha, R_{PC}) .

\end{eqnarray}
Similar expressions can be readily written for the $y$ and $z$ components of the centers.

\section*{\sffamily \Large 3-Index Integral Derivatives in the OS Scheme}
Using the derivaive properties of GTOs we can write:

\begin{eqnarray}
\frac{\partial (ab|C) }{\partial A_\delta}  & = &  2a (a+1_\delta, b|C) \ - \ i_a  (a-1_\delta, b|C) ,
\label{os:grad:a} 

\end{eqnarray}
and

\begin{eqnarray}
\frac{\partial (ab|C) }{\partial B_\delta}  & = &  2b (a, b+1_\delta,|C) \ - \ i_b  (a, b-1_\delta|C) ,
\label{os:grad:b} 

\end{eqnarray}
where $\delta = x,y,z$. Hence, following prescription of Head-Gordon and Pople,\hyperref[csl:8]{(\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Head-Gordon}} \bibnamefont{and}, n.d.)} we can generate the all derivative equations from Eq.\eqref{os:grad:a} and \eqref{os:grad:b}.

\section*{\sffamily \Large Density-Fitting}
The AO basis integrals can be cast into the following form with the help of DF approach:

\begin{eqnarray}
(\mu \nu | \lambda \sigma)_{DF} & = & \sum_{Q}^{N_{aux}} b_{\mu \nu}^{Q} b_{\lambda \sigma}^{Q},

\end{eqnarray}
The DF tensors $b_{\mu \nu}^{Q}$ are defined as follows:

\begin{eqnarray}
b_{\mu \nu}^{Q} & = & \sum_{P}^{N_{aux}} (\mu \nu | P) [{\bf J}^{-1/2}]_{PQ},
\label{eq:df:factors} 

\end{eqnarray}
where

\begin{eqnarray}
(\mu \nu | P) & = & \int\int \chi_{\mu}({\bf r_{1}}) \chi_{\nu}({\bf r_{1}}) \frac{1}{r_{12}} \varphi_{P}({\bf r_{2}})  \,d{\bf r_{1}} d{\bf r_{2}} ,

\end{eqnarray}
and

\begin{eqnarray}
J_{PQ} & = & \int\int  \varphi_{P}({\bf r_{1}}) \frac{1}{r_{12}} \varphi_{Q}({\bf r_{2}})  \,d{\bf r_{1}} d{\bf r_{2}}, 

\end{eqnarray}
where $\chi_{\mu}({\bf r})$ and $\varphi_{P}({\bf r})$ are the primary and auxiliary basis set members, respectively.

\section*{\sffamily \Large Integral Direct Fock}
In the AO basis, the Fock matrix can be written as follows:

\begin{eqnarray}
f_{\mu\nu} & = & h_{\mu\nu} \ + \ \sum_{\lambda,\sigma} D_{\lambda \sigma} \big[ (\mu \nu | \lambda \sigma) \ - \ \frac{1}{2}(\mu \lambda | \nu \sigma) \big] ,

\end{eqnarray}
where $\bm{D}$ is the density matrix. For example, for the closed-shell case:

\begin{eqnarray}
D_{\mu\nu} & = & 2\sum_{i}^{occ} C_{\mu i} C_{\nu i} ,

\end{eqnarray}
where $C_{\mu i} $ are the MO coefficients.

Two-electron part of the Fock matrix consists of so called $J$- and $K$-terms. The $J$-term can be written as:

\begin{eqnarray}
J_{\mu\nu} & = &  \sum_{\lambda,\sigma} D_{\lambda \sigma}  (\mu \nu | \lambda \sigma) .

\end{eqnarray}

Applying the DF technique we obtain the following final formula for $J_{\mu\nu}$:

\begin{eqnarray}
J_{\mu\nu} & = &  \sum_{P}^{N_{aux}} (\mu \nu | P) D^{\prime\prime\prime}_{P} ,

\end{eqnarray}
where,

\begin{eqnarray}
D^{\prime\prime\prime}_{P}  & = &  \sum_{R}^{N_{aux}} [{\bf J}^{-1/2}]_{PR} D^{\prime\prime}_{R}, 

\end{eqnarray}

\begin{eqnarray}
D^{\prime\prime}_{R}  & = &  \sum_{Q}^{N_{aux}}  [{\bf J}^{-1/2}]_{RQ} D^{\prime}_{Q} ,

\end{eqnarray}

\begin{eqnarray}
D^{\prime}_{Q} & = &  \sum_{\lambda,\sigma} D_{\lambda \sigma}  (\lambda \sigma | Q) .

\end{eqnarray}
Hence, we can compute the 3-index integrals on the fly and contract them with the density matrix to obtain $J_{\mu\nu}$.

The $K$-term can be written as:

\begin{eqnarray}
K_{\mu\nu} & = &  \sum_{\lambda,\sigma} D_{\lambda \sigma} (\mu \lambda | \nu \sigma),

\end{eqnarray}

We can compute the $K$-term in a \emph{semi-direct} algorithm with three steps, as follows:

\begin{itemize}
\item Step-1:

\begin{eqnarray}
(Q | \nu i)  & = &  \sum_{\sigma}  C_{\sigma i} (Q | \nu \sigma)  \quad cost=O*N^2*N_{aux} ,

\end{eqnarray}

\item Step-2:

\begin{eqnarray}
C_{\nu i}^{P}  & = &  \sum_{Q}^{N_{aux}}  [{\bf J}^{-1/2}]_{PQ} (Q | \nu i) \quad cost=O*N*N_{aux}^2 ,

\end{eqnarray}

\item Step-3:

\begin{eqnarray}
K_{\mu\nu} & = &  \sum_{P}^{N_{aux}}  \sum_{i}^{occ}  C_{\mu i}^{P}  C_{\nu i}^{P}  \quad cost=O*N^2*N_{aux} .

\end{eqnarray}

\end{itemize}
The scaling of this algorithm is $O(N^4)$ and memory usage is $2*O*N*N_{aux}$.

\section*{\sffamily \Large Analytic Gradients}
For a general density-fitted wave function based method, the analytic gradient equation in the AO basis can be written as follows\hyperref[csl:42]{(\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 1241)}

\begin{eqnarray}
\dod{E}{x} &=& \sum_{\mu\nu} \gamma_{\mu \nu}  h_{\mu \nu}^{x}  \ - \ \sum_{\mu\nu} F_{\mu \nu} S_{\mu \nu}^{x}
\ + \ \sum_{Q}^{N_{aux}} \sum_{\mu \nu} \Gamma_{\mu \nu}^{Q} (Q|\mu \nu)^{x} 
\ - \  \sum_{P,Q}^{N_{aux}} \Gamma_{PQ} {\bf J}_{PQ}^{x} .
\label{dfgrad:ao} 

\end{eqnarray}
where $\gamma_{\mu \nu}$ is the one-particle density matrix (OPDM), $F_{\mu \nu}$ is the generalized Fock matrix (GFM), $\Gamma_{PQ}$ and $\Gamma_{\mu \nu}^{Q}$ are the 2- and 3-index two-particle density matrices (TPDM), $h_{\mu \nu}^{x}$, $S_{\mu \nu}^{x}$, ${\bf J}_{PQ}^{x}$, and $(Q|\mu \nu)^{x} $ are the core Hamiltonian, overlap, metric, and 3-index integral derivatives, respectively. All explicit equations for the PDMs and GFM were presented in our previous papers.\hyperref[csl:42]{(\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 1241}; \hyperref[csl:50]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}} \bibnamefont{and}, n.d.}; \hyperref[csl:51]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}} \bibnamefont{and}, n.d.}; \hyperref[csl:48]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, n.d.)}

\section*{\sffamily \Large {\sc Molint} Framework}
The {\sc Molint} framework consists of  {\sc Molint}, {\sc Tensor}, {\sc Molecule}, {\sc Options}, {\sc Master}, and {\sc Basis} libraries (Figure \ref{fig:molint:frame}). Before illustrating the usage of  the {\sc Molint} framework, let us introduce the core libraries.

\subsection*{{\sc Tensor} Library}
The {\sc Tensor} library, which is written in {\sc C++},  consists of six classes. These are {\tt Tensor1d}, {\tt Tensor2d}, and {\tt Tensor3d} for one-, two-, and three-dimensional arrays of {\tt double}, respectively, and {\tt Tensor1i}, {\tt Tensor2i}, and {\tt Tensor3i} for one-, two-, and three-dimensional arrays of {\tt integer}, respectively. The {\sc Tensor} library includes numerous functions to handle, manipulate, and operate large arrays of pointers. As an example, we illustrate some prominent features of the {\tt Tensor2d} class.

The {\tt Tensor2d} class is developed to perform various tensor operations in a convenient way. It handles 2, 3, and 4 dimensional tensors, and store them like a matrix. Further, the class take advantages of  \emph{shared pointers} . The memory allocation for a 2D tensor can be performed with the {\tt Tensor2d}  class as follows:

\begin{description}
\item  {\tt Tensor2d(std::string name,  \textcolor{mydarkgreen}{int} d1, \textcolor{mydarkgreen}{int} d2);}

\end{description}
where, {\tt name}: The lable of tensor, {\tt d1} and {\tt d2} are the first and second dimensions of the tensor. An example of memory allocation is:

\begin{description}
\item {\tt SharedTensor2d CmoA = \\ 
SharedTensor2d(\textcolor{orange}{new} Tensor2d("\textcolor{myred}{Alpha MO Coefficients}", nao, nmo));}

\end{description}
where, {\tt CmoA}  is the MO coefficient matrix, {\tt nao} and {\tt nmo} are the number of the AOs and MOs, respectively.

Similarly, 3- and 4-dimensional arrays can be allocated as follows:

\begin{description}
\item {\tt K = SharedTensor2d(\textcolor{orange}{new} Tensor2d("\textcolor{myred}{DFBCC B (Q|IA)}", nQ, naoccA, navirA));}
\item {\tt T = SharedTensor2d(\textcolor{orange}{new} Tensor2d("\textcolor{myred}{T2 <Ij|Ab>}", naoccA, naoccB, navirA, navirB));}

\end{description}

Some important functions and example usages are:

\begin{itemize}
\item \textbf{Print:} {\tt \textcolor{mydarkgreen}{void} print(type output);}\\
\textbf{Example:} {\tt T->print(output);} where {\tt output} can be type of {\tt std::string}, {\tt const char*}, and {\tt FILE*}.


\item \textbf{Set value:} {\tt \textcolor{mydarkgreen}{void} set(\textcolor{mydarkgreen}{int} i, \textcolor{mydarkgreen}{int} j, \textcolor{mydarkgreen}{double} value);}\\
\textbf{Example:} {\tt T->set(3, 5, 12.0);}

\item \textbf{Set Tensor:} {\tt \textcolor{mydarkgreen}{void} set(SharedTensor2d \&A);}\\
\textbf{Example:} {\tt T->set(A);}

\item \textbf{Add value:} {\tt \textcolor{mydarkgreen}{void} add(\textcolor{mydarkgreen}{int} i, \textcolor{mydarkgreen}{int} j, \textcolor{mydarkgreen}{double} value);}\\
\textbf{Example:} {\tt T->add(3, 5, 12.0);}

\item \textbf{Add Tensor:} {\tt \textcolor{mydarkgreen}{void} add(SharedTensor2d \&A);}\\
\textbf{Example:} {\tt T->add(A);}

\item \textbf{Scale:} {\tt \textcolor{mydarkgreen}{void} scale(\textcolor{mydarkgreen}{double} alpha);} where {\tt alpha} is the scaling value.\\
\textbf{Example:} {\tt T->scale(2.0);}

\item \textbf{AXPY:} The {\tt axpy} function is used to add/subtract a SharedTensor via {\sc Blas} routine.\\
{\tt \textcolor{mydarkgreen}{void} axpy(\textcolor{mydarkgreen}{const} SharedTensor2d \&A, \textcolor{mydarkgreen}{double} alpha);}\\
\textbf{Example:} {\tt T->axpy(2.0);}

\item \textbf{Transpose:} {\tt \textcolor{mydarkgreen}{void} trans(\textcolor{mydarkgreen}{const} SharedTensor2d \&A);}\\
\textbf{Example:} {\tt T->trans(A);}

\item \textbf{Diagonalize:} {\tt \textcolor{mydarkgreen}{void} diagonalize(\textcolor{mydarkgreen}{const} SharedTensor2d \&eigvectors, \\ 
\textcolor{mydarkgreen}{const} SharedTensor1d \&eigvalues, \textcolor{mydarkgreen}{double} cutoff, \textcolor{mydarkgreen}{bool} ascending);}\\
\textbf{Example:} {\tt FockA->diagonalize(UeigA, eigA, 1.0E-10, \textcolor{myred}{true});}

\item \textbf{Linear equation solver for ${\bf A x}  =   {\bf b} $:} \\
{\tt \textcolor{mydarkgreen}{void} cdgesv(\textcolor{mydarkgreen}{const} SharedTensor1d \&b, \textcolor{mydarkgreen}{int} errcode);} \\
\textbf{Example:} {\tt A->cdgesv(b, errcode);}

\item \textbf{Matrix multiplication of ${\bf C} =   \alpha {\bf A} {\bf B} + \beta {\bf C}$ type:} \\
{\tt \textcolor{mydarkgreen}{void} gemm(\textcolor{mydarkgreen}{bool} transa, \textcolor{mydarkgreen}{bool} transb, \textcolor{mydarkgreen}{const} SharedTensor2d \&A, \\
\textcolor{mydarkgreen}{const} SharedTensor2d \&B, \textcolor{mydarkgreen}{double} alpha, \textcolor{mydarkgreen}{double} beta);}\\
\textbf{Example:} {\tt C->gemm(\textcolor{myred}{false}, \textcolor{myred}{false}, A, B, 1.0, 0.0);}

\item \textbf{General contraction of $C_{m,n} =   \alpha \sum_{k} A_{m,k} B_{k,n} + \beta C_{m,n}$ type:} \\
{\tt \textcolor{mydarkgreen}{void} contract(\textcolor{mydarkgreen}{bool} transa, \textcolor{mydarkgreen}{bool} transb, \textcolor{mydarkgreen}{int} m, \textcolor{mydarkgreen}{int} n, \textcolor{mydarkgreen}{int} k,\\ \textcolor{mydarkgreen}{const} SharedTensor2d \&A, \textcolor{mydarkgreen}{const} SharedTensor2d \&B, \textcolor{mydarkgreen}{double} alpha, \textcolor{mydarkgreen}{double} beta);}\\
\textbf{Example:} Let us try to evaluate the contraction of $x_{ij}^{ab} =  \sum_{e} t_{ij}^{ae} F_{be}$, where $i,j$ active occupied and $a,b,e$ are active virtual indices:\\
{\tt X->contract(\textcolor{myred}{false}, \textcolor{myred}{true}, naoccA*navirA*naoccA, navirA, navirA, T, F, 1.0, 0.0);}

\item \textbf{Transformation of ${\bf C} =   {\bf L}^T {\bf A} {\bf L} $ type:} \\
{\tt \textcolor{mydarkgreen}{void} transform(\textcolor{mydarkgreen}{const} SharedTensor2d \&A, \textcolor{mydarkgreen}{const} SharedTensor2d \&L);}\\
\textbf{Example:} {\tt C->transform(A, L);}

\item \textbf{Dot product of $value  = \sum_{mnk\dots}  A_{mnk\dots} B_{mnk\dots}$ type:} \\
{\tt \textcolor{mydarkgreen}{double} vector\_dot(\textcolor{mydarkgreen}{const} SharedTensor2d \&rhs);}\\
\textbf{Example:} {\tt \textcolor{mydarkgreen}{double} value = A->vector\_dot(B);}

\item \textbf{Orthogonalization via modified Gram-Schmid algorithm:} \\
{\tt \textcolor{mydarkgreen}{void} mgs();}\\
\textbf{Example:} {\tt A->mgs();}


\end{itemize}

\subsection*{{\sc Options} Library}
The {\sc Options} library ({\tt liboptions}) read the input file, default is {\tt input.inp}, which includes geometry and the MFW options. The MFW accepts the following user options:

\begin{itemize}
\item \textbf{\tt nocom} : Do not shift the geometry to the center-of-mass.{\tt Type}: Boolean, {\tt Default}: false.

\item \textbf{\tt puream} : Do use pure angular momentum functions?{\tt Type}: Boolean, {\tt Default}: true.

\item \textbf{\tt print\_level} : This options is a flag to control printing extra information.{\tt Type}: integer, {\tt Default}: 0.



\item \textbf{\tt geom\_units} : Unit of the geometry. {\tt Type}: string, {\tt Possible values}: ang (default), bohr.

\item \textbf{\tt integral\_cutoff} : Integral cutoff value. {\tt Type}: double, {\tt Default}: 1.0E-10.

\item \textbf{\tt boys\_cutoff} : The cutoff value used in the evaluation of the Boys function. {\tt Type}: double, {\tt Default}: 1.0E-15.

\item \textbf{\tt boys\_zero} : The computational zero used in the evaluation of the Boys function. {\tt Type}: double, {\tt Default}: 1.0E-6.

\item \textbf{\tt delta\_x} : The step size in the Boys grid. {\tt Type}: double, {\tt Default}: 0.1.

\item \textbf{\tt x\_max} : The maximum value of $x$ for the Boys grid. {\tt Type}: int, {\tt Default}: 30.

\item \textbf{\tt boys\_maxiter} : The maximum number of iterations in computing values of the Boys grid. {\tt Type}: int, {\tt Default}: 1000.

\item \textbf{\tt taylor\_order} : The order of Taylor expansion in computing the Boys function. {\tt Type}: int, {\tt Default}: 6.

\item \textbf{\tt boys\_method} : Method to compute the Boys grid. {\tt Type}: string, {\tt Possible values}: ghp (default), saunders.


\end{itemize}

An input example for the {\tt liboptions} is given as follows:

\begin{description}
\item

{\tt \quad \quad \$set\_ints \\
\quad nocom true \\
\quad integral\_cutoff 1.0e-10 \\
\quad boys\_cutoff 1.0e-15 \\
\quad delta\_x 0.1 \\
\quad x\_max 30 \\
\quad taylor\_order 6 \\
\quad print\_level 0 \\
\quad boys\_method ghp \\
\$end\_ints}\\

\end{description}

\subsection*{{\sc Molecule} Library}
The {\sc Molecule} library ({\tt libmolecule}) read geometry, charge, and multiplicity information, and pass all necessary information to {\sc Libmolint}, which includes molecular integrals. An example input is given as follows:

\begin{description}
\item
{\tt ~~~~\$set\_molecule \\
0 1 \\
O          0.000000000000      0.000000000000    $-$0.065775570538 \\
H          0.000000000000    $-$0.759061990794     0.521953018295 \\
H          0.000000000000     0.759061990794     0.521953018295 \\
\$end\_molecule} \\

\end{description}
The current version of the {\sc Molecule} library accepts Cartesian coordinates only. The {\tt libmolecule} is directly called by the {\tt libmolint}, hence; there is no need to call it when interfacing the MFW

\subsection*{{\sc Master} Library}
The {\sc Master} library ({\tt libmaster}) is built for the purpose of the communication within the MFW.  In the MFW, users directly call {\tt libmolint} only, and the constructor function creates the {\tt libmaster} object, let us denote it {\tt master}.  Then, {\tt master} gets all necessary options using the {\tt liboptions} and store them. Similarly, the {\tt master} get all necessary molecular information, such as coordinates, inter atomic distances, nuclear repulsion energy and gradients, and store them. Users can access the {\tt master} object through the {\tt libmolint}, and can get all information stored in it. Let us denote the object of the {\tt libmolint} as {\tt ints}. Then, users can access the {\tt master} as follows:

\begin{description}
\item {\tt SharedMaster master = ints->master();}

\end{description}
The object {\tt master} is reset in the destructor function of the {\tt libmolint}.

\subsection*{{\sc Basis} Library}
The {\sc Basis} library ({\tt libbasis}) is a built-in library, which reads in all required basis set(s) information for a given molecule. The {\tt libbasis} is directly called by in the {\tt libmolint}, hence; there is no need to call it when interfacing the MFW. Further,  {\tt libbasis} includes common basis functions in a  {\tt share} directory in the {\sc Gaussian} format. One can readily add a new basis set as a  {\tt .gbs} file. The address of the {\tt share} directory is defined as an environmental variable {\tt MFWDATADIR}. Hence, users should export it, for example:

\begin{description}
\item {\tt export MFWDATADIR=/Users/ugur/lib/molintfw/share}

\end{description}

\subsection*{{\sc Molint} Library}
The usage and features of the core libraries were explained in the previous sections. In here, we focus on the {\tt Molint} library ({\tt libmolint}), which is capable of computing overlap, kinetic, potential, dipole, metric, and 3-index integrals, as well as their first derivatives (Figure \ref{fig:molint:features}).

Now, let us demonstrate how to use the MFW. The constructor function for the {\tt libmolint} is described as follows:

\begin{description}
\item  {\tt Integrals(std::string file\_name, std::string basis, std::string aux)}

\end{description}
where {\tt file\_name} is the name of input file, {\tt basis} and {\tt aux} are the names of the primary and auxiliary basis sets, for which there should be {\tt basis.gbs} and {\tt aux.gbs} files, respectively. At first we need to construct the {\tt libmolint} object as follows:

\begin{description}
\item {\tt SharedIntegrals ints = \\ 
SharedIntegrals(new Integrals("input.inp", "cc-pvdz", "cc-pvdz-ri"));}

\end{description}
Basic information such as the number of atoms, atomic orbitals, auxiliary basis functions, and shells can be obtained from the object as follows:

\begin{description}
\item {\tt natom = ints->natom();} 
\item {\tt nao = ints->nao();} 
\item {\tt naux = ints->naux();} 
\item {\tt nshell = ints->nshell();} 

\end{description}

\subsubsection*{Common Integrals}
The overlap integrals can be obtained form the object as follows:

\begin{description}
\item {\tt Sao->copy(ints->compute\_overlap\_ints());} 

\end{description}
Of course, the {\tt Sao} array should be formed previously, such as follows:

\begin{description}
\item {\tt auto Sao = SharedTensor2d(new Tensor2d("AO-basis Overlap Ints", nao, nao));} 

\end{description}

Similarly, the kinetic and potential integrals can be obtained as follows:

\begin{description}
\item {\tt Tao->copy(ints->compute\_kinetic\_ints());} 

\end{description}

\begin{description}
\item {\tt Vao->copy(ints->compute\_potential\_ints());} 

\end{description}
where {\tt Tao} and {\tt Vao} are two-dimensional tensors (Tensor2d).

The dipole integrals can be obtained as follows:

\begin{description}
\item {\tt Mdm = ints->compute\_dipole\_ints();} 

\end{description}
where {\tt Mdm} is a three-dimensional tensor (Tensor3d), which can be formed as follows:

\begin{description}
\item {\tt  auto Mdm = SharedTensor3d(new Tensor3d("Mu (x,y,z|mu,nu)", 3, nao, nao));} 

\end{description}

The metric integrals can be computed as follows:

\begin{description}
\item {\tt  ints->compute\_metric\_ints();} 
\item {\tt  Jmetric->copy(ints->metric());}

\end{description}
where {\tt Jmetric} is a two-dimensional tensor of size {\tt naux $\times$ naux}. One may obtain $[{\bf J}^{-1/2}]_{PQ}$ directly instead of metric integrals as follows:

\begin{description}
\item {\tt  Jmhalf->copy(ints->metric\_mhalf());}

\end{description}

The DF factors can be computed in two steps. In the first step one can compute three-index integrals as follows:

\begin{description}
\item {\tt  ints->compute\_three\_index\_ints();} 

\end{description}
and in the second step the DF factors can be computed, according to Eq.\eqref{eq:df:factors}, as follows:

\begin{description}
\item {\tt ints->compute\_df\_factors();} 

\end{description}
Then, the DF factors can be obtained as follows:

\begin{description}
\item {\tt Bdf->copy(ints->df\_factors());} 

\end{description}
where {\tt Bdf} is a two-dimensional tensor with row and column dimensions of {\tt naux} and  {\tt nao*nao}, respectively. Similarly, the three-index integrals can be obtained as follows:

\begin{description}
\item {\tt  Qmn->copy(ints->three\_index\_tei());} 

\end{description}
However, in order to proceed in this way, one should compute metric integrals before the DF factors. Further, one may combine the computation of metric and DF factors, as follows:

\begin{description}
\item {\tt ints->compute\_df\_integrals();} 

\end{description}

\subsubsection*{Density-Fitted Fock Matrix}
So far, we have assume that there is enough memory to keep 3-index integrals in the core memory. However, for larger computations one may prefer an integral direct approach, where $J$- and $K$-terms of the Fock matrix can be computed. In our implementation we have a full direct approach for the $J$-term and a semi-direct approach for the $K$-term. The memory requirement of the semi-direct $K$ (SDK) approach is $2*O*N*N_{aux}$, which is significantly low (by $\frac{N}{2O}$) compare with $N^2*N_{aux}$ for the incore case. In our $JK$ algorithm we build the both terms simultaneously. For example, for the restricted HF (RHF) case:

\begin{description}
\item {\tt  ints->compute\_dffock\_JK(noccA, FaoA, CmoA, DaoA, Jmhalf);} 

\end{description}
where {\tt noccA} is the number of $\alpha$ spin occupied orbitals, {\tt FaoA} is the AO basis Fock matrix, {\tt CmoA} is the MO coefficients matrix,  {\tt DaoA} is the density matrix, and {\tt Jmhalf} is $\bm{J}^{-1/2}$.

For the unrestricted HF (UHF) case we need to call {\tt libmolint} twice:

\begin{description}
\item {\tt  ints->compute\_dffock\_JK(noccA, FaoA, CmoA, Dtot, Jmhalf);} 
\item {\tt  ints->compute\_dffock\_JK(noccB, FaoB, CmoB, Dtot, Jmhalf);} 

\end{description}
where {\tt Dtot} is the total density matrix.

\subsubsection*{Integral Derivatives}
The core-Hamiltonian term of analytic gradients in Eq.\eqref{dfgrad:ao} $\big(\sum_{\mu\nu} \gamma_{\mu \nu}  h_{\mu \nu}^{x} \big)$ can be computed as follows:

\begin{description}
\item {\tt  EgradH->copy(ints->compute\_kinetic\_grad(G1ao));} 
\item {\tt  EgradH->add(ints->compute\_potential\_grad(G1ao));} 

\end{description}
where {\tt G1ao} is the AO basis OPDM ($\gamma_{\mu \nu}$) and {\tt EgradH} is a two-dimensional tensor, which can be allocated as follows:

\begin{description}
\item {\tt  SharedTensor2d EgradH = \\ SharedTensor2d(new Tensor2d("-Core Hamiltonian Gradient:", natom, 3));} 

\end{description}
where {\tt natom} is the total number of atoms.

The overlap gradient $\big(-\sum_{\mu\nu} F_{\mu \nu} S_{\mu \nu}^{x}\big)$ can be computed as follows:

\begin{description}
\item {\tt  EgradS->copy(ints->compute\_overlap\_grad(GFao));} 
\item {\tt  EgradS->scale(-1.0);} 

\end{description}
where {\tt GFao} is the AO basis GFM ($F_{\mu \nu}$) and {\tt EgradS} is a two-dimensional tensor, which can be allocated as follows:

\begin{description}
\item {\tt  SharedTensor2d EgradS = \\ SharedTensor2d(new Tensor2d("-Overlap Gradient:", natom, 3));} 

\end{description}

The two-electron part of analytic gradients should be computed in two steps since we employ different auxiliary basis sets for the reference and correlation energies. The details of this separation procedure were given in our previous studies.\hyperref[csl:42]{(\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 1241}; \hyperref[csl:47]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 4389)} The reference and separable parts (RefSep) of 2- and 3-index PDMs can be merged before the contraction. The RefSep contribution for the metric gradients $\big( -\sum_{P,Q}^{JKFIT} \Gamma_{PQ}^{(RefSep)} {\bf J}_{PQ}^{x} \big)$ can be computed as follows:
\textbf{\begin{description}
\item {\tt EgradM->copy(ints->compute\_metric\_grad(Gaux\_ref));} 
\item {\tt  EgradM->scale(-1.0);} 
\end{description}} 
where {\tt Gaux\_ref} is the RefSep part of the 2-index TPDM ($\Gamma_{PQ}^{(RefSep)} $) and {\tt EgradM} is a two-dimensional tensor.

The RefSep contribution for the 3-index gradient $\big(\sum_{Q}^{JKFIT} \sum_{\mu \nu} \Gamma_{\mu \nu}^{Q} (Q|\mu \nu)^{x} \big)$ can be computed as follows:
\textbf{\begin{description}
\item {\tt Egrad3I->copy(ints->compute\_3index\_grad(gQao\_ref));} 
\end{description}} 
where {\tt gQao\_ref} is the RefSep part of the 3-index TPDM ($\Gamma_{\mu \nu}^{Q(RefSep)}$) and {\tt Egrad3I} is a two-dimensional tensor.

The correlation parts of the metric and 3-index gradients can be computed similarly. One needs to call the same function but with the correlation parts of 2- and 3-index TPDMs. Finally all gradient components are merged in to the final gradient.

\section*{\sffamily \Large Verification of the {\sc Molint} Framework}
The MFW code has been verified with respect to the {\sc Psi4}\hyperref[csl:57]{(\bibinfo{author}{\bibfnamefont{D.~G.~A.} \bibnamefont{Smith}}, n.d.)} package for the computation of energies, dipole moments, and analytic gradients for various density-fitted methods, such as the Hartree-Fock (DF-HF), second-order perturbation theory (DF-MP2), coupled-cluster singles and doubles (DF-CCSD), and the CCSD with perturbative triples [DF-CCSD(T)]. All comparisons were performed with my {\sc Dfocc} module,\hyperref[csl:62]{(\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 2371}; \hyperref[csl:42]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 1241}; \hyperref[csl:47]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 4389}; \hyperref[csl:36]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 1179}; \hyperref[csl:37]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, n.d.}; \hyperref[csl:50]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}} \bibnamefont{and}, n.d.}; \hyperref[csl:63]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, 1441}; \hyperref[csl:51]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}} \bibnamefont{and}, n.d.}; \hyperref[csl:48]{\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bozkaya}}, n.d.)} which is available in the {\sc Psi4}\hyperref[csl:57]{(\bibinfo{author}{\bibfnamefont{D.~G.~A.} \bibnamefont{Smith}}, n.d.)} package, as well as it has an external version.

\section*{\sffamily \Large Illustrative Applications}
Results from the DF-MP2 method were obtained for a set of alkanes for the assessment of computational cost for single-point energy and analytic gradient computations. For the alkanes set, Dunning's correlation-consistent polarized valence triple-, quadruple-, and quintuple-$\zeta$ basis sets (cc-pVTZ, cc-pVQZ, and cc-pV5Z) were employed with the frozen core approximation.\hyperref[csl:64]{(\bibinfo{author}{\bibfnamefont{T.~H.} \bibnamefont{Dunning}}, n.d.}; \hyperref[csl:65]{\bibinfo{author}{\bibfnamefont{D.~E.} \bibnamefont{Woon}} \bibnamefont{and}, n.d.)} For the cc-pVXZ (X=T, Q, 5) primary basis sets, cc-pVXZ-JKFIT\hyperref[csl:26]{(\bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Weigend}}, n.d.)} and cc-pVXZ-RI\hyperref[csl:66]{(\bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Weigend}}, n.d.)} auxiliary basis sets were employed for reference and correlation energies, respectively.  The DF-MP2 computations were performed with the external version of the {\sc Dfocc} module.

We consider a set of alkanes (C$_n$H$_{2n+2}$, $n=1-10$) to assess the efficiency of the MFW.  The total computational (wall) time for single-point frozen-core energy and analytic gradient computations using the DF-MP2 method are presented in Tables \ref{tab:mp2:tz:e}--\ref{tab:mp2:5z:g}.  These computations were performed on an Intel(R) Xeon(R) CPU E5-2620 v4 $@$ 2.10 GHz computer ($\sim$ 120 GB memory). In the MFW, 3-index integrals are parallelized using {\sc OpenMP} API in the {\sc Fortran} side, while all other integrals use only 1 core. The computations of overlap, kinetic, dipole, and potential integrals are already very fast; hence, there is no need to consider parallelism for them. However, the computational time for metric integrals also includes the cost of diagonalization and $\bm{J}^{-1/2}$ formation. Hence, our metric code is partially parallelized, through basic linear algebra subprogram routines (BLAS and LAPACK).

Tables \ref{tab:mp2:tz:e}--\ref{tab:mp2:5z:e} reports computational (wall) time for single-point frozen-core DF-MP2/cc-pVXZ (X=T,Q,5) energy computations with 1 and 8 cores. For the largest system considered, \ce{C7H16} with cc-pV5Z/cc-pV5Z-RI, there are 1517 primary and  2807 auxiliary basis functions. For this system, the computation of metric integrals and the formation of $\bm{J}^{-1/2}$ is performed just in 15 and 8 seconds, with 1 and 8 cores, respectively. Hence, the computations of metric and one-electron integrals are already very fast, and we need to focus on the cost of 3-index integrals, which are the time determining step in molecular integrals. The number of primary and auxiliary basis functions ($N_{bf}/N_{aux}$) in the largest systems considered are 608/1470 (cc-pVTZ/cc-pVTZ-RI), 1210/2530 (cc-pVQZ/cc-pVQZ-RI), and 1517/2807 (cc-pV5Z/cc-pV5Z-RI). For these systems wall time for serial (1 core) and parallel (8 cores) computations are 92.84 and 27.13 (cc-pVTZ), 630.15 and 166.49 (cc-pVQZ), and 1523.50 and 322.42 (cc-pV5Z) seconds.  For the largest system considered, 1517 primary and  2807 auxiliary basis functions in the cc-pV5Z/cc-pV5Z-RI pair, the 3-index integrals are computed just in 5.4 minutes, with 8 cores, on our moderately fast server. Furthermore, for this system, our parallel code is 4.7-fold faster than the serial version. These results demonstrate that the MFW is satisfactorily efficient for the computation of 3-index integrals.

The MFW is capable of computing analytic gradients providing the AO basis PDMs and GFM. Tables \ref{tab:mp2:tz:g}--\ref{tab:mp2:5z:g} reports computational (wall) time for single-point frozen-core DF-MP2/cc-pVXZ (X=T,Q,5) analytic gradient computations with 2 cores. Our first derivative code has not been parallelized yet. However, in the case of 2 or more cores, we compute derivatives for two different centers simultaneously. For the largest system considered, \ce{C5H12} with cc-pV5Z/cc-pV5Z-RI, there are 1115 primary and 2057 auxiliary basis functions. For this system, the computation of the metric term of the analytic gradient is performed just in 8 seconds. Hence, the computations of the metric and one-electron terms are already very fast, and we need to focus on the cost of the 3-index term, which is the time determining step in analytic gradient evaluation. The number of primary and auxiliary basis functions ($N_{bf}/N_{aux}$) in the \ce{C5H12} molecule are 318/765 (cc-pVTZ/cc-pVTZ-RI), 635/1320 (cc-pVQZ/cc-pVQZ-RI), and 1115/2057 (cc-pV5Z/cc-pV5Z-RI). For these systems wall time, with 2 cores, are 50.18 (cc-pVTZ), 358.90 (cc-pVQZ), and 2285.34 (cc-pV5Z) seconds.  These results demonstrate that the MFW is reasonably efficient for the evaluation of the analytic gradients for the density-fitted methods.

\section*{\sffamily \Large Future Directions}
The first version of MFW, {\sc Molint 1.0}, includes the essential integrals and their first derivatives, as well as integral direct/semi-direct Fock and analytic gradients. In future releases, we intend to include more challenging features such as the second derivatives, R12 and/or F12 integrals, and gauge including atomic orbital (GIAO) integrals. Furthermore, a parallel API framework for the graphical process unit (GPU) is in our ``to do" list. However, all these developments require a significant effort; hence, we believe that the MFW will cover these features in the next years.

\section*{\sffamily \Large Conclusions}
In this research, an application programming interface (API) framework, denoted {\sc Molint} (MFW),  for the computation of molecular integrals and their first derivatives, over contracted Gaussian functions, for the density-fitted methods has been reported. The MFW includes overlap, dipole, kinetic, potential, metric, and 3-index integrals as well as their first derivatives. Furthermore, the MFW provides a smooth approach to build the Fock matrix and evaluate analytic gradients. The MFW is a {\sc C++/Fortran} hybrid code, which can take advantage of shared-memory parallel programming techniques. Moreover, the MFW is free software, which can be obtained through an e-mail request from the developer ({\tt http://www.bozkayalab.hacettepe.edu.tr/en/menu/molint-47}).  The illustrative applications have demonstrated that the MFW is an efficient and user-friendly API for the computation of molecular integrals and their first derivatives. Consequently, considering both the computational efficiency and an easy-to-use structure of the MFW, we conclude that the MFW emerges as a very useful tool for computational quantum chemistry.

\section*{Acknowledgements}
I thank Yavuz Alag\"{o}z and Bet\"{u}l Ermi\c{s} for their assistance to run test computations.

((Additional Supporting Information may be found in the online version of this article.))

\clearpage

\clearpage

\begin{sidewaystable}\selectlanguage{english}
\begin{table}[H] %TZ
\caption{{The number of primary and auxiliary basis functions ($N_{bf}$  and $N_{aux}$) in the cc-pVTZ and cc-pVTZ-RI basis sets, wall time (in seconds) for the metric and 3-index integrals for DF-MP2 energy computations. All computations were performed on a single node (up to 8 cores) Intel(R) Xeon(R) CPU E5-2620 v4 $@$ 2.10 GHz computer (memory $\sim$ 120 GB).}} 
\begin{tabular}{l cccccc}
\hline\hline 
\textbf{Molecule} & $\bm{N_{bf}}$ & $\bm{N_{aux}}$ & \textbf{Metric (1 core)} & {\textbf{Metric (8 cores)}} & \textbf{3-Index (1 core)} & {\textbf{3-Index (8 cores)}} \\
\hline
\ce{CH4}    & 86                   & 201              & 0.02                     & 0.03                      & 0.26                      & 0.12                       \\
\ce{C2H6}   & 144                  & 342              & 0.07                     & 0.08                      & 1.24                      & 0.41                       \\
\ce{C3H8}   & 202                  & 483              & 0.15                     & 0.13                      & 3.44                      & 1.00                       \\
\ce{C4H9}   & 260                  & 624              & 0.26                     & 0.21                      & 7.32                      & 2.16                       \\
\ce{C5H10}  & 318                  & 765              & 0.44                     & 0.32                      & 13.36                     & 3.98                       \\
\ce{C6H14}  & 376                  & 906              & 0.64                     & 0.45                      & 22.05                     & 6.51                       \\
\ce{C7H16}  & 434                  & 1047             & 0.91                     & 0.60                      & 33.73                     & 9.65                       \\
\ce{C8H18}  & 492                  & 1188             & 1.24                     & 0.78                      & 48.98                     & 14.06                      \\
\ce{C9H20}  & 550                  & 1329             & 1.67                     & 0.99                      & 68.43                     & 19.60                      \\
\ce{C10H22} & 608                  & 1470             & 2.15                     & 1.22                      & 92.84                     & 27.13                     \\
\hline\hline
\end{tabular}
\label{tab:mp2:tz:e}

\end{table}

\end{sidewaystable}

\clearpage

\begin{sidewaystable}\selectlanguage{english}
\begin{table}[H] %QZ
\caption{{The number of primary and auxiliary basis functions ($N_{bf}$  and $N_{aux}$) in the cc-pVQZ and cc-pVQZ-RI basis sets, wall time (in seconds) for the metric and 3-index integrals for DF-MP2 energy computations. All computations were performed on a single node (up to 8 cores) Intel(R) Xeon(R) CPU E5-2620 v4 $@$ 2.10 GHz computer (memory $\sim$ 120 GB).}} 
\begin{tabular}{l cccccc}
\hline\hline 
\textbf{Molecule} & $\bm{N_{bf}}$ & $\bm{N_{aux}}$ & \textbf{Metric (1 core)} & {\textbf{Metric (8 cores)}} & \textbf{3-Index (1 core)} & {\textbf{3-Index (8 cores)}} \\
\hline
\ce{CH4}    & 175                  & 352              & 0.00                     & 0.11                      & 0.03                      & 0.50                       \\
\ce{C2H6}   & 290                  & 594              & 0.28                     & 0.23                      & 8.67                      & 2.33                       \\
\ce{C3H8}   & 405                  & 836              & 0.61                     & 0.45                      & 23.81                     & 6.11                       \\
\ce{C4H9}   & 520                  & 1078             & 1.09                     & 0.75                      & 50.71                     & 12.60                      \\
\ce{C5H10}  & 635                  & 1320             & 2.48                     & 2.13                      & 92.45                     & 22.99                      \\
\ce{C6H14}  & 750                  & 1562             & 2.68                     & 1.62                      & 152.45                    & 39.18                      \\
\ce{C7H16}  & 865                  & 1804             & 3.90                     & 2.22                      & 233.55                    & 60.32                      \\
\ce{C8H18}  & 980                  & 2046             & 5.53                     & 2.98                      & 339.10                    & 88.72                      \\
\ce{C9H20}  & 1095                 & 2288             & 7.49                     & 3.76                      & 467.37                    & 119.32                     \\
\ce{C10H22} & 1210                 & 2530             & 10.09                    & 4.96                      & 630.15                    & 166.49                    \\
\hline\hline
\end{tabular}
\label{tab:mp2:qz:e}

\end{table}

\end{sidewaystable}

\clearpage

\begin{sidewaystable}\selectlanguage{english}
\begin{table}[H] %5Z
\caption{{The number of primary and auxiliary basis functions ($N_{bf}$  and $N_{aux}$) in the cc-pV5Z and cc-pV5Z-RI basis sets, wall time (in seconds) for the metric and 3-index integrals for DF-MP2 energy computations. All computations were performed on a single node (up to 8 cores) Intel(R) Xeon(R) CPU E5-2620 v4 $@$ 2.10 GHz computer (memory $\sim$ 120 GB).}} 
\begin{tabular}{l cccccc}
\hline\hline 
\textbf{Molecule} & $\bm{N_{bf}}$ & $\bm{N_{aux}}$ & \textbf{Metric (1 core)} & {\textbf{Metric (8 cores)}} & \textbf{3-Index (1 core)} & {\textbf{3-Index (8 cores)}} \\
\hline
\ce{CH4}    & 311                  & 557              & 0.47                     & 0.31                                             & 12.34                     & 2.90                                              \\
\ce{C2H6}   & 512                  & 932              & 0.97                     & 0.70                                             & 57.32                     & 12.80                                             \\
\ce{C3H8}   & 713                  & 1307             & 2.03                     & 1.45                                             & 157.29                    & 38.98                                             \\
\ce{C4H9}   & 914                  & 1682             & 3.63                     & 2.30                                             & 329.34                    & 75.08                                             \\
\ce{C5H10}  & 1115                 & 2057             & 6.19                     & 3.56                                             & 601.72                    & 132.32                                            \\
\ce{C6H14}  & 1316                 & 2432             & 9.98                     & 5.33                                             & 991.67                    & 211.05                                            \\
\ce{C7H16}  & 1517                 & 2807             & 15.25                    & 8.05                                             & 1523.50                   & 322.42                                           \\
\hline\hline
\end{tabular}
\label{tab:mp2:5z:e}

\end{table}

\end{sidewaystable}

\clearpage\selectlanguage{english}
\begin{table}[hbt!] %TZ-grad
\caption{{The number of primary and auxiliary basis functions ($N_{bf}$  and $N_{aux}$) in the cc-pVTZ and cc-pVTZ-RI basis sets, wall time (in seconds) for the metric and 3-index terms in DF-MP2 analytic gradient computations. All computations were performed on a single node (2 cores) Intel(R) Xeon(R) CPU E5-2620 v4 $@$ 2.10 GHz computer (memory $\sim$ 120 GB).}} 
\begin{tabular}{l cccc}
\hline\hline 
\textbf{Molecule} & $\bm{N_{bf}}$ & $\bm{N_{aux}}$ & \textbf{Metric} & \textbf{3-Index}  \\
\hline
\ce{CH4}    & 86                   & 201              & 0.05                   & 0.89                    \\
\ce{C2H6}   & 144                  & 342              & 0.17                   & 4.55                    \\
\ce{C3H8}   & 202                  & 483              & 0.31                   & 12.73                   \\
\ce{C4H9}   & 260                  & 624              & 0.49                   & 27.68                   \\
\ce{C5H10}  & 318                  & 765              & 0.73                   & 50.18                  \\
\hline\hline
\end{tabular}
\label{tab:mp2:tz:g}

\end{table}

\clearpage\selectlanguage{english}
\begin{table}[hbt!] %QZ-grad
\caption{{The number of primary and auxiliary basis functions ($N_{bf}$  and $N_{aux}$) in the cc-pVQZ and cc-pVQZ-RI basis sets, wall time (in seconds) for the metric and 3-index terms in DF-MP2 analytic gradient computations. All computations were performed on a single node (2 cores) Intel(R) Xeon(R) CPU E5-2620 v4 $@$ 2.10 GHz computer (memory $\sim$ 120 GB).}} 
\begin{tabular}{l cccc}
\hline\hline 
\textbf{Molecule} & $\bm{N_{bf}}$ & $\bm{N_{aux}}$ & \textbf{Metric} & \textbf{3-Index}  \\
\hline
\ce{CH4}    & 175                  & 352              & 0.19                   & 8.59                    \\
\ce{C2H6}   & 290                  & 594              & 0.52                   & 32.78                   \\
\ce{C3H8}   & 405                  & 836              & 1.03                   & 92.24                   \\
\ce{C4H9}   & 520                  & 1078             & 1.71                   & 196.61                  \\
\ce{C5H10}  & 635                  & 1320             & 2.61                   & 358.90                 \\
\hline\hline
\end{tabular}
\label{tab:mp2:qz:g}

\end{table}

\clearpage\selectlanguage{english}
\begin{table}[hbt!] %5Z-grad
\caption{{The number of primary and auxiliary basis functions ($N_{bf}$  and $N_{aux}$) in the cc-pV5Z and cc-pV5Z-RI basis sets, wall time (in seconds) for the metric and 3-index terms in DF-MP2 analytic gradient computations. All computations were performed on a single node (2 cores) Intel(R) Xeon(R) CPU E5-2620 v4 $@$ 2.10 GHz computer (memory $\sim$ 120 GB).}} 
\begin{tabular}{l cccc}
\hline\hline 
\textbf{Molecule} & $\bm{N_{bf}}$ & $\bm{N_{aux}}$ & \textbf{Metric} & \textbf{3-Index}  \\
\hline
\ce{CH4}    & 311                  & 557              & 0.54                   & 45.85                   \\
\ce{C2H6}   & 512                  & 932              & 1.57                   & 220.37                  \\
\ce{C3H8}   & 713                  & 1307             & 3.21                   & 607.48                  \\
\ce{C4H9}   & 914                  & 1682             & 5.31                   & 1257.25                 \\
\ce{C5H10}  & 1115                 & 2057             & 8.03                   & 2285.34                \\
\hline\hline
\end{tabular}
\label{tab:mp2:5z:g}

\end{table}

\clearpage\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-molint-frame/default-figure}
\caption{{\label{fig:molint:frame} Components of the {\sc Molint} framework.%
}}
\end{center}
\end{figure}

\clearpage\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-molint-features/default-figure}
\caption{{\label{fig:molint:features} Capabilities of the {\sc Molint} framework.%
}}
\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\section*{References}\sloppy
\phantomsection
\label{csl:1}

\phantomsection
\label{csl:2}

\phantomsection
\label{csl:3}

\phantomsection
\label{csl:4}

\phantomsection
\label{csl:5}

\phantomsection
\label{csl:6}

\phantomsection
\label{csl:7}

\phantomsection
\label{csl:8}

\phantomsection
\label{csl:9}

\phantomsection
\label{csl:10}

\phantomsection
\label{csl:11}

\phantomsection
\label{csl:12}

\phantomsection
\label{csl:13}(5198).

\phantomsection
\label{csl:14}

\phantomsection
\label{csl:15}

\phantomsection
\label{csl:16}

\phantomsection
\label{csl:17}

\phantomsection
\label{csl:18}

\phantomsection
\label{csl:19}

\phantomsection
\label{csl:20}

\phantomsection
\label{csl:21}

\phantomsection
\label{csl:22}

\phantomsection
\label{csl:23}

\phantomsection
\label{csl:24}

\phantomsection
\label{csl:25}

\phantomsection
\label{csl:26}

\phantomsection
\label{csl:27}

\phantomsection
\label{csl:28}

\phantomsection
\label{csl:29}

\phantomsection
\label{csl:30}

\phantomsection
\label{csl:31}

\phantomsection
\label{csl:32}

\phantomsection
\label{csl:33}

\phantomsection
\label{csl:34}

\phantomsection
\label{csl:35}

\phantomsection
\label{csl:36}(1179).

\phantomsection
\label{csl:37}

\phantomsection
\label{csl:38}(4415).

\phantomsection
\label{csl:39}

\phantomsection
\label{csl:40}

\phantomsection
\label{csl:41}

\phantomsection
\label{csl:42}(1241).

\phantomsection
\label{csl:43}

\phantomsection
\label{csl:44}

\phantomsection
\label{csl:45}

\phantomsection
\label{csl:46}

\phantomsection
\label{csl:47}(4389).

\phantomsection
\label{csl:48}

\phantomsection
\label{csl:49}

\phantomsection
\label{csl:50}

\phantomsection
\label{csl:51}

\phantomsection
\label{csl:52}

\phantomsection
\label{csl:53}

\phantomsection
\label{csl:54}

\phantomsection
\label{csl:55}(2020).

\phantomsection
\label{csl:56}

\phantomsection
\label{csl:57}

\phantomsection
\label{csl:58}

\phantomsection
\label{csl:59}

\phantomsection
\label{csl:60}

\phantomsection
\label{csl:61}

\phantomsection
\label{csl:62}(2371).

\phantomsection
\label{csl:63}(1441).

\phantomsection
\label{csl:64}

\phantomsection
\label{csl:65}

\phantomsection
\label{csl:66}

\end{CJK}\end{document}

