\documentclass[aip,jcp,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,longtable,wrapfig,bbold,siunitx,xspace}
\usepackage[version=4]{mhchem}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}

\usepackage{hyperref}
\hypersetup{
    colorlinks,
    linkcolor={red!50!black},
    citecolor={red!70!black},
    urlcolor={red!80!black}
}

\newcommand{\mc}{\multicolumn}
\newcommand{\fnm}{\footnotemark}
\newcommand{\fnt}{\footnotetext}
\newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}

\newcommand{\titou}[1]{\textcolor{red}{#1}}
\newcommand{\toto}[1]{\textcolor{blue}{#1}}
\newcommand{\joonho}[1]{\textcolor{purple}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{r\ed}{\sout{#1}}}
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}
\newcommand{\AS}[1]{\toto{(\underline{\bf AS}: #1)}}
\newcommand{\JL}[1]{\joonho{(\underline{\bf JL}: #1)}}

\usepackage{listings}
\definecolor{codegreen}{rgb}{0.58,0.4,0.2}
\definecolor{codegray}{rgb}{0.5,0.5,0.5}
\definecolor{codepurple}{rgb}{0.25,0.35,0.55}
\definecolor{codeblue}{rgb}{0.30,0.60,0.8}
\definecolor{backcolour}{rgb}{0.98,0.98,0.98}
\definecolor{mygray}{rgb}{0.5,0.5,0.5}

\definecolor{sqred}{rgb}{0.85,0.1,0.1}
\definecolor{sqgreen}{rgb}{0.25,0.65,0.15}
\definecolor{sqorange}{rgb}{0.90,0.50,0.15}
\definecolor{sqblue}{rgb}{0.10,0.3,0.60}

\lstdefinestyle{mystyle}{
    backgroundcolor=\color{backcolour},
    commentstyle=\color{codegreen},
    keywordstyle=\color{codeblue},
    numberstyle=\tiny\color{codegray},
    stringstyle=\color{codepurple},
    basicstyle=\ttfamily\footnotesize,
    breakatwhitespace=false,
    breaklines=true,
    captionpos=b,
    keepspaces=true,
    numbers=left,
    numbersep=5pt,
    numberstyle=\ttfamily\tiny\color{mygray},
    showspaces=false,
    showstringspaces=false,
    showtabs=false,
    tabsize=2
  }

  \newcolumntype{d}{D{.}{.}{-1}}

\lstset{style=mystyle}

% addresses
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
\newcommand{\Vienna}{Vienna}

\begin{document}	

\title{Stochastically accelerated perturbative triples in coupled cluster calculations}

% Alphabetic order
\author{Yann \surname{Damour}}
	\affiliation{\LCPQ}
\author{Alejandro \surname{Gallo}}
	\affiliation{\Vienna}
\author{Andreas \surname{Irmler}}
	\affiliation{\Vienna}
\author{Andreas \surname{Gr\"uneis}}
	\affiliation{\Vienna}
\author{Anthony \surname{Scemama}}
	\email{scemama@irsamc.ups-tlse.fr}
	\affiliation{\LCPQ}
	
\begin{abstract}
We introduce a novel algorithm that leverages stochastic sampling
techniques to approximate perturbative triples correction in the
coupled-cluster (CC) framework.
By combining elements of randomness and determinism, our algorithm
achieves a favorable balance between accuracy and computational cost.
The main advantage of this algorithm is that it allows for
calculations to be stopped at any time, providing an unbiased
estimate, with a statistical error that goes to zero as the exact
calculation is approached.
We provide evidence that our semi-stochastic algorithm achieves
substantial computational savings compared to traditional
deterministic methods.
Specifically, we demonstrate that a precision of 0.5 milliHartree can
be attained with only 10\% of the computational effort required by the
full calculation.
This work opens up new avenues for efficient and accurate
computations, enabling investigations of complex molecular systems
that were previously computationally prohibitive.
\bigskip
\begin{center}
%	\boxed{\includegraphics[width=0.5\linewidth]{TOC}}
\end{center}
\bigskip
\end{abstract}

\maketitle

%=================================================================%
\section{Introduction}
\label{sec:introduction}

Coupled cluster (CC) theory is a powerful quantum mechanical approach widely used in computational chemistry and physics to describe the electronic structure of atoms, molecules, and materials.
It offers a systematic and rigorous framework for accurate predictions of molecular properties and reactions by accounting for electron correlation effects beyond the mean-field approximation.
Among the various variants of the CC method, the CCSD(T) method stands as the gold standard of quantum chemistry.
CCSD(T), which incorporates singles, doubles, and a perturbative treatment of triples, has demonstrated exceptional accuracy and reliability, making it one of the preferred choices for benchmark calculations and highly accurate predictions.
The CCSD(T) method has found successful applications in a diverse range of areas, including spectroscopy,\cite{villa_2011,watson_2016,vilarrubias_2020} reaction kinetics,\cite{dontgen_2015,castaneda_2012} and materials design,\cite{zhang_2019} and has played a pivotal role in advancing our understanding of complex chemical phenomena.

In the context of CC theory, perturbative triples represent an important contribution to the accuracy of electronic structure calculations.\cite{stanton_1997}
However, the computational cost associated with their calculation can be prohibitively high, especially for large molecular systems.
The CCSD(T) method, which includes the perturbative treatment of triples, is known to have a computational scaling of $\order{N^7}$, where $N$ represents the system size.
This scaling can rapidly become impractical, posing significant challenges in terms of computational resources and time requirements.

To address this computational bottleneck, our goal is to develop a novel semi-stochastic algorithm that brings back the computational time to a level comparable to that of the CCSD method, which has a scaling of $\order{N^6}$, while ensuring well-controlled approximations.
Our algorithm strikes a balance between computational efficiency and
accuracy, making calculations for larger basis sets more feasible without compromising precision.
By incorporating stochastic sampling techniques, our approach provides an alternative avenue for approximating perturbative triples, relieving the computational burden inherent in traditional deterministic methods. This not only reduces the computational time to a more favorable level but also preserves the parallelism capabilities of CC calculations, ensuring efficient utilization of computational resources.

In the following sections of this paper, we will provide a brief introduction to the computation of perturbative triples in coupled cluster theory. We will explain the principles of our semi-stochastic algorithm, outlining its key features and advantages. Additionally, we will present implementation details, discussing the technical aspects and considerations involved in the algorithm's practical realization. To demonstrate the effectiveness and applicability of our approach, we finally present illustrative examples that showcase the algorithm's performance and compare it with the conventional algorithm.

%=================================================================%
\section{Theoretical Background}
\label{sec:theory}

%a. Brief overview of coupled cluster theory
%b. Perturbative triples in coupled cluster theory
%c. Challenges and limitations of traditional approaches


The perturbative triples correction,
\begin{equation}
E_{(T)}  =  \sum_{ijkabc} E_{ijk}^{abc} 
\end{equation}
is a sum of $N=N_o^3 \times N_v^3$ terms,
\begin{equation}
 E_{ijk}^{abc}  =  \frac{(4 W_{ijk}^{abc} +
              W_{ijk}^{bca} + W_{ijk}^{cab})
              (V_{ijk}^{abc} - V_{ijk}^{cba})}{\epsilon_i + \epsilon_j + \epsilon_k -
\epsilon_a - \epsilon_b - \epsilon_c} 
\end{equation}
which depend on the canonical orbital energies $\epsilon$, and the tensors $W$ and $V$.
The indices $i,j,k$ and $a,b,c$ denote respectively occupied and virtual orbitals.

The bottleneck is the computation of $W$, which requires $\order{N_o^3 \times
N_v^4}$ operations. However, most of the operations involved in the computation of $W$
can be recast into matrix multiplications, which are among the most efficient
operations than can be executed on modern CPUs and
accelerators.\cite{ma_2011,haidar_2015,dinapoli_2014,springer_2018}

%=================================================================%
\section{Semi-Stochastic Algorithm}
\label{sec:algorithm}

%a. Explanation of the algorithm's main principles and methodology
% - Rewriting with nV in outer-most loops
% - Implies summing all permutations of (a,b,c) to avoid extra dgemms
%  - Comparison in number of Gflops to theory
%b. Description of the stochastic sampling procedure
%    - Memoization
%    - Importance function
%    - Semi-stochastic algorithm
%c. Discussion of the algorithm's advantages and potential trade-offs
%d. Detailed pseudocode or algorithmic steps, if applicable

\subsection{Stochastic formulation}

\subsection{Test code}
\label{subsec:test_code}
% Include the test code here, if applicable.


%=================================================================%
\section{Implementation Details}
\label{sec:implementation}

%a. Description of the computational framework and software used
%b. Discussion of any specific optimizations or parallelization techniques employed
% - Explain that form_w and form_v can be entirely produced by dgemm
% - Show implementation with OpenMP tasks
%c. Any relevant technical considerations or limitations
% - Limitations: memory because in-core algorithm.

%=================================================================%
\section{Examples of applications}

%a. Presentation of benchmark systems and datasets used for evaluation
% - Benzene TZ
% - Streptocyanine QZ: Small molecule in a large basis set
% - Caffeine def2-svp: Large molecule in a small basis set
% - Vibrational frequency of F2/cc-pvqz
%b. Discussion of the obtained results, comparing against other methods
% - Measure flops and compare to the peak
%c. Analysis of the algorithm's accuracy, efficiency, and scalability
%d. Discussion of any observed limitations or challenges

\subsection{Vibrational frequency of \ce{F2}}

In this example, we compute the vibrational frequency of \ce{F2} by
computing the potential energy curve, and fitting it with a Morse
potential
\begin{equation}
  E(r) = D_e \left( 1 - e^{-a (r - r_e)} \right)^2 + E(r_e)
\end{equation}
where $E(r)$ is the energy at distance $r$, $D_e$ is the well depth,
$r_e$ is the equilibrium bond distance, and $a$ is a parameter
controlling the width of the potential well.
The vibrational frequency $\nu$ is calculated as
\begin{equation}
  \nu = \frac{1}{2 \pi c} \sqrt{\frac{2D_e a^2}{\mu}
\end{equation}
where $\mu$ is the mass of the Fluorine atom, and $c$ is the speed of
light in cm/s.

% CCSD
%a               = 2.2936           +/- 0.006318     (0.2755%)
%De              = 0.125888         +/- 0.0005213    (0.4141%)
%re              = 1.3893           +/- 0.0003428    (0.02468%)
%E0              = -199.338         +/- 6.422e-05    (3.222e-05%)

% CCSD(T) exact
%a               = 2.65592          +/- 0.0403       (1.518%)
%De              = 0.0718253        +/- 0.001879     (2.617%)
%re              = 1.4105           +/- 0.00215      (0.1524%)
%E0              = -199.358         +/- 0.0003179    (0.0001595%)

\section{Conclusion}
\label{sec:conclusion}

%a. Summary of the algorithm and its advantages
%b. Recapitulation of the key findings and contributions
%c. Final remarks and encouragement for further research

%=================================================================%

%%%%%%%%%%%%%%%%%%%%%%%%
\acknowledgements{
This work was supported by the European Centre of Excellence in Exascale Computing TREX --- Targeting Real Chemical Accuracy at the Exascale. 
This project has received funding from the European Union's Horizon 2020 — Research and Innovation program --- under grant agreement No.~952165.
This work was performed using HPC resourced from CALMIP (Toulouse) under allocations p18005 and p22001.}
%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Data availability statement}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The data that supports the findings of this study are available within the article and its supplementary material.

%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{stochastic_triples}
%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}
