Chimie Théorique » scripts_chimie4psmn

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% This is a (brief) model paper using the achemso class

%% The document class accepts keyval options, which should include

%% the target journal and optionally the manuscript type.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\documentclass[journal=jacsat,manuscript=article,draft]{achemso}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Place any additional packages needed here.  Only include packages

%% which are essential, to avoid problems later.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\usepackage{chemformula} % Formula subscripts using \ch{}

\usepackage[T1]{fontenc} % Use modern font encodings

\usepackage{amsmath} % collection de symboles mathématiques

\usepackage{amssymb} % collection de symboles mathématiques

\usepackage{mathtools} % \mathclap for long subscipts under sums

\newcommand{\lgi}[1]{\limits_{\mathclap{\substack{#1}}}}

\DeclarePairedDelimiter\bra{\langle}{\rvert}

\DeclarePairedDelimiter\ket{\lvert}{\rangle}

\DeclarePairedDelimiterX\braket[2]{\langle}{\rangle}{#1 \delimsize\vert #2}

\usepackage{dsfont} % mathds for unitary matrix

% define "struts", as suggested by Claudio Beccari in

%    a piece in TeX and TUG News, Vol. 2, 1993.

\newcommand\Tstrut{\rule{0pt}{2.6ex}}         % = `top' strut

\newcommand\Bstrut{\rule[-0.9ex]{0pt}{0pt}}   % = `bottom' strut

\usepackage{subcaption} % subfigure

\usepackage[utf8]{inputenc}       % utilisation directe des caractères accentués sur pc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% If issues arise when submitting your manuscript, you may want to

%% un-comment the next line.  This provides information on the

%% version of every file you have used.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%\listfiles

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Place any additional macros here.  Please use \newcommand* where

%% possible, and avoid layout-changing macros (which are not used

%% when typesetting).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand*\mycommand[1]{\texttt{\emph{#1}}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Meta-data block

%% ---------------

%% Each author should be given as a separate \author command.

%%

%% Corresponding authors should have an e-mail given after the author

%% name as an \email command. Phone and fax numbers can be given

%% using \phone and \fax, respectively; this information is optional.

%%

%% The affiliation of authors is given after the authors; each

%% \affiliation command applies to all preceding authors not already

%% assigned an affiliation.

%%

%% The affiliation takes an option argument for the short name.  This

%% will typically be something like "University of Somewhere".

%%

%% The \altaffiliation macro should be used for new address, etc.

%% On the other hand, \alsoaffiliation is used on a per author basis

%% when authors are associated with multiple institutions.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\author{Stephan Steinmann}

\email{stephan.steinmann@ens-lyon.fr}

\author{Ruben Staub}

\email{ruben.staub@ens-lyon.fr}

%\phone{+123 (0)123 4445556}

%\fax{+123 (0)123 4445557}

\affiliation[ENS de Lyon]

{Laboratoire de Chimie, ENS de Lyon, Lyon}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% The document title should be given as usual. Some journals require

%% a running title from the author: this should be supplied as an

%% optional argument to \title.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title[Extension of ALMO formalism to metallic systems in CP2K]

  {Extension of ALMO formalism to metallic systems in CP2K\footnote{A footnote for the title}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Some journals require a list of abbreviations or keywords to be

%% supplied. These should be set up here, and will be printed after

%% the title and author information, if needed.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\abbreviations{IR,NMR,UV}

\keywords{American Chemical Society, \LaTeX}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% The manuscript does not need to include \maketitle, which is

%% executed automatically.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% The "tocentry" environment can be used to create an entry for the

%% graphical table of contents. It is given here as some journals

%% require that it is printed as part of the abstract page. It will

%% be automatically moved as appropriate.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{tocentry}

Some journals require a graphical entry for the Table of Contents.

This should be laid out ``print ready'' so that the sizing of the

text is correct.

Inside the \texttt{tocentry} environment, the font used is Helvetica

8\,pt, as required by \emph{Journal of the American Chemical

Society}.

The surrounding frame is 9\,cm by 3.5\,cm, which is the maximum

permitted for  \emph{Journal of the American Chemical Society}

graphical table of content entries. The box will not resize if the

content is too big: instead it will overflow the edge of the box.

This box and the associated title will always be printed on a

separate page at the end of the document.

\end{tocentry}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% The abstract environment will automatically gobble the contents

%% if an abstract is not used by the target journal.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstract}

  In the DFT formalism, charge transfer is usually poorly described, since DFT methods have a tendency to overestimate electron delocalization.

This issue can be solved with the use of localized orbitals within the ALMO (also called BLW) formalism, which forces electrons to remain localized on predefined fragments.

This method has already been implemented in the CP2K modeling tool. However, it's current development and implementation is not compatible with partially occupied orbitals (i.e. electronic smearing). Therefore, metallic systems cannot be currently handled with this formalism.

Here we show that the ALMO formalism can be practically extended to metallic systems, under a reasonable approximation.

We found that common simplifications used in mixed-state theory cannot be applied within the ALMO formalism, which makes the exact approach unpractical, as it is a combinatorial problem with exponential complexity. However, under a basic mean field approximation, we were able to unify mixed-state theory with ALMO formalism into a simple, practical formulation, that have been implemented in CP2K as S-ALMO.

This formulation is of critical importance for electro-catalysis applications, where we have both metals and electrolytes whose charge must remain localized, for example. Another application of this work is to allow charge transfer analysis for metallic systems.

We believe that this tool will provide a deep understanding of charge transfer impact on catalytic processes, especially for alloys where much is yet to be discovered. This could even represent a first step in the rational design of alloys catalysts, a major challenge for green chemistry development.

\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Start the main part of the manuscript here.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

\subsection{Context}

\paragraph{}The Absolutely Localized Molecular Orbital (ALMO) formalism was developed with the aim to force electrons to remain localized on predefined fragments. This idea has been proposed by multiple authors (Stoll\cite{stoll_use_1980}, Nagata\cite{nagata_perturbation_2004}, Gianinetti\cite{gianinetti_modification_1996}) and is known under various names in the literature (BLW, ELMO\cite{couty_extremely_1997}\cite{genoni_novel_2005}, \ldots).

\paragraph{}ALMO formalism has been implemented\cite{khaliullin_efficient_2013}\cite{khaliullin_efficient_2006} in multiple codes, including the open source simulation software CP2K\cite{noauthor_cp2k_nodate}.

The most important applications are:

\begin{itemize}

\item Charge transfer analysis\cite{steinmann_norbornene_2010}\cite{mo_charge_2004}\cite{fornili_determination_2003}\cite{khaliullin_analysis_2008} (by providing currently one of the most rigorous way to switch on/off specific charge transfers), and more generally energy decomposition analysis\cite{mo_energy_2000}\cite{khaliullin_unravelling_2007}. Alternatives to this formalism can be found with Constrained DFT\cite{kaduk_constrained_2012} (CDFT). However, ALMO formalism provides generally with a more convenient framework to work with, and lower computational cost for large systems.

\item Computational speed-up when dealing with large systems, leading to potential linear scaling.

\item Design of transferable densities/orbitals (for hybrid QM/MM calculations\cite{fornili_suitability_2006}, and X-Ray structure interpretation\cite{genoni_x-ray_2013}).

\item Reduced Basis Set Superposition Error (BSSE, due to electron delocalization in finite basis computations, leading to energy mismatch).

\end{itemize}

\paragraph{}In the DFT formalism (see Annexes section \ref{annexes}), electron delocalization has a tendency to be over-estimated\cite{steinmann_why_2012}, leading to artificial fractional charges for ions at large distances\cite{cohen_challenges_2012}. For example, a partial unphysical charge transfer from PF$_6^-$ to Bu$_4$N$^+$ could occur in solution. This is particularly problematic for electrochemical simulation where many ions are present. ALMO, where such electron delocalization can be suppressed thus provides a valuable methodology in this context.

\subsection{Problematic}

\paragraph{}Unfortunately, in its current development, ALMO formalism is incompatible with metallic systems (and more generally, any system whose ground-state cannot be described by a single quantum state).

To be precise, mixed-state theory is required to properly describe metallic systems. However, ALMO has not yet been developed within the mixed-state theory.

\paragraph{}This is particularly problematic for heterogeneous catalysis, and especially electrocatalysis (where both the electrolytes and the metallic part are present).

\paragraph{}Therefore, we are interested here in unifying ALMO formalism and mixed-state theory.

\subsection{Mixed-state theory}

\paragraph{}Mixed-state theory is required to investigate the ground-state of metallic systems, and more generally any system with high degeneracy near the Fermi level. Indeed, the ground-state of such a system cannot be described by a single quantum state anymore, and a statistical ensemble of accessible quantum states must be considered: a mixed-state.

\subsubsection{Mixed-state algebra}

\paragraph{}Let $\hat{A}$ be an observable, and $\Omega$ is the set of all possible quantum states of a given system, the expectation value $A^{avg}$ associated with $\hat{A}$ on this system is:

\begin{equation}

A^{avg} = \sum\lgi{\mathcal{S}\in\Omega}\bra{\Psi_\mathcal{S}}\hat{A}\ket{\Psi_\mathcal{S}}p_\mathcal{S}

\end{equation}

where $\Psi_\mathcal{S}$ is the wavefunction associated with the quantum state $\mathcal{S}$. $p_\mathcal{S}$ is the probability that the real system is in the quantum state $\mathcal{S}$, based on the energy of $\mathcal{S}$.

\paragraph{}For the ALMO formalism, we are mainly interested in the reformulation of the $1$-electron density $\hat{\rho}$. In the Hartree-Fock approximation, and in the case of orthonormal molecular orbitals, the $1$-electron density $\hat{\rho}$ can be rewritten:

\begin{equation}

\hat{\rho} = \sum\lgi{\text{occ }\psi_i\in\Psi}\ket{\psi_i}\!\bra{\psi_i}

\end{equation}

\paragraph{}The mixed-state $1$-electron density operator $\hat{\rho}^{avg}$ can then be written:

\begin{equation}

\hat{\rho}^{avg} = \sum\lgi{\mathcal{S}\in\Omega}p_\mathcal{S}\hat{\rho}_{\mathcal{S}} = \sum\lgi{\mathcal{S}\in\Omega}p_\mathcal{S}\sum\lgi{\psi_i\in\Psi_\mathcal{S}}\ket{\psi_i}\!\bra{\psi_i}

\end{equation}

\paragraph{}Since the same global set of orbitals is used to construct all the quantum states of the system (see section \ref{mo_state}), this sum can be rewritten in terms of orbitals instead of quantum-states. In other words, instead of considering the population of the quantum-states, one can consider the population of the orbitals:

\begin{equation}

\hat{\rho}^{avg} = \sum\lgi{i}\ket{\psi_i}\!\bra{\psi_i}p_i

\end{equation}

where the sum is over all orbitals (solutions of the Fock equation), and $p_i$ is the probability for the orbital $\psi_i$ to be occupied (based on the energy of the orbital).

\subsubsection{Mathematical trick: rescaled orbitals}

\paragraph{}The mixed-state $1$-electron density can be rewritten with modified orbitals:

\begin{equation}

\hat{\rho}^{avg} = \sum\lgi{i}\ket{\psi_i}\!\bra{\psi_i}p_i = \sum\lgi{i}\ket{\sqrt{p_i}\psi_i}\!\bra{\sqrt{p_i}\psi_i} = \sum\lgi{i}\ket{\psi_i^\prime}\!\bra{\psi_i^\prime}

\end{equation}

where $\psi_i^\prime = \sqrt{p_i}\psi_i$ is called a rescaled orbital.

\paragraph{}Therefore, with these rescaled orbitals, the mixed-state $1$-electron density can be computed as usual\footnote{This trick even works for virtual orbital dependent approaches such as MP2 and RPA\cite{yang_extension_2013}.}. In other words, to compute the mixed-state $1$-electron density, one can simply rescale each orbital $\psi_i$ by $\sqrt{p_i}$ to obtain rescaled orbitals $\psi_i^\prime$, and then compute the density as usual (as with a single quantum-state) by using the rescaled orbitals in the equations instead of the original ones. Instead of changing the equations, the orbitals are changed. This rescaling trick saves therefore a lot of computational time, since we do not need to consider separately every single quantum-state possible anymore (all states are basically considered at once, since we use orbitals that are common to all those states).

\subsection{ALMO formalism presentation}

\paragraph{}The ALMO (Absolutely Localized Molecular Orbitals) formalism has several names in the literature (BLW (Block Localized Wavefunctions), ELMO (Extremely Localized Molecular Orbitals), \ldots), all based on the same idea: rewrite the Hatree-Fock equations to allow the use of localized molecular orbitals.

\subsubsection{Block localization}

\paragraph{}In order to define localized orbitals, one must first define blocks within which orbitals are to be localized.

\paragraph{Definition:}A block $B_i$ (or domain, fragment, \ldots) is defined as an exclusive set of basis set functions $\{\phi_a^i\}_a$, such that each basis set function is associated with exactly one block. In other words, blocks can be seen as partitions of the set of basis set functions.

It is common that a block represents a meaningful fragment, as a molecule. In that case, the set of basic set functions $\{\phi_a^i\}_a$ associated with the block $B_i$ is the union of all basis set functions $\{\phi_k^j\}_k$ used to describe each atom $j$ of $B_i$.

\begin{equation}

\{\phi_{a}^i\}_a = \bigcup\lgi{\text{atom }j \in B_i} \ \{\phi_k^j\}_k

\end{equation}

It is important to note that blocks are well defined only in the case of finite basis set\footnote{In practice, even with large basis sets, ALMO is still applicable.}, and this partition requires localized basis functions (i.e. plane waves should not be used).

\paragraph{Definition:}A block-localized molecular orbital is a linear combination of basis set functions associated with the same block. Such orbital will be simply referred as localized orbital in this work. In other words, a localized orbital is an orbital whose basis set function expansion is restrained to a given block.

Therefore, an occupied orbital $\psi_{a}^i$ localized on block $B_i$ is written:

\begin{equation}

\psi_{a}^i = \sum\lgi{k}T_{ik,ia}\phi_{k}^i

\end{equation}

where $T$ is the orbital coefficient matrix, and $\{\phi_{k}^i\}_k$ is the set of basis set functions associated with block $B_i$.

\paragraph{}As a consequence, the orbital coefficient matrix $T$ has a block-diagonal structure (i.e. $\forall a,k\quad i \neq j \Rightarrow T_{ik,ja}=0$):

\begin{equation}

T = \left(\begin{matrix}

T^1 & 0 & \cdots & 0 \\

0 & T^2 & \cdots & 0 \\

\vdots & \vdots & \ddots & \vdots \\

0 & 0 & \cdots & T^n \\

\end{matrix}\right)

\end{equation}

where $T^i$ is the orbital coefficient matrix restricted to the block $B_i$.

%T_{11,11} & T_{11,12} & \cdots & T_{11,1n_1} & 0 & \cdots & 0 \\

%T_{12,11} & T_{12,12} & \cdots & T_{12,1n_1} & 0 & \cdots & 0 \\

%\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots \\

%T_{1m_1,11} & T_{1m_1,12} & \cdots & T_{1m_1,1n_1} & 0 & \cdots & %0 \\

\paragraph{Note:}Thanks to this sparsity, matrix multiplications involving $T$ are particularly less costly in computation time. This contributes to the global efficiency of the ALMO formalism when it comes to large systems.

\subsubsection{Non-orthogonal orbitals}

\paragraph{}Most complications arising from the use of localized orbitals are due to the fact that localized orbitals are inherently non-orthogonal.

\begin{equation}

\braket{\psi_i}{\psi_j} = \sigma_{ij} \neq 0

\end{equation}

where $\sigma$ is called the overlap matrix.

Indeed, due to degree of freedom considerations, locality and orthogonality constraints cannot be both satisfied in the general case. However, one can require that orthogonality is preserved within each block, without any loss of generality (i.e. $\sigma$ has identity-like diagonal blocks). We will work within this additional constraint in this study.

\subsubsection{Updated formulae}

\paragraph{}We will not enter here in the details of the derivations (that can be found in the original paper from Stoll\cite{stoll_use_1980}), but instead provide with the main ideas and results.

\paragraph{}In Stoll's algorithm, the localized occupied orbitals are orthonormalized\footnote{Orthonormalized orbitals are no longer localized.} using the Löwdin orthonormalization:

\begin{equation}

\ket{\psi_i^L} = \sum\lgi{\text{occ }j} \ket{\psi_j}\sigma^{-\frac{1}{2}}_{ij},\quad T^L = T\sigma^{-\frac{1}{2}}

\end{equation}

where $\sigma^{-\frac{1}{2}}$ is the square root inverse of the overlap matrix. In terms of orbital coefficient matrices, it is equivalent to define a new coefficient matrix $T^L$ for the orthonormal orbitals $\left\{\ket{\psi_i^L}\right\}_i$. Therefore, the usual machinery (e.g. gradient, forces computation, \ldots) can be used with these orthonormal orbitals.

\paragraph{}Since square root computations of matrices are costly, in practice we work with reciprocal occupied orbitals $\ket{\tilde{\psi_i}}$, such that $\braket{\tilde{\psi_i}}{\psi_j} = \delta_{ij}$:

\begin{equation}

\ket{\tilde{\psi_i}} = \sum\lgi{\text{occ }j} \ket{\psi_j}\sigma^{-1}_{ij},\quad \tilde{T} = T\sigma^{-1}

\end{equation}

where $\tilde{T}$ is the coefficient matrix of reciprocal orbitals. 

\paragraph{}Once the orbitals are orthonormalized, the $1$-electron density is computed as usual:

\begin{equation}

\hat{\rho} = \sum\lgi{\text{occ }i}\ket{\psi_i^L}\!\bra{\psi_i^L} = \sum\lgi{\text{occ }i}\ket{\tilde{\psi_i}}\!\bra{\psi_i} = \sum\lgi{\text{occ }i}\sum\lgi{\mu,\nu}\tilde{T}_{\mu i}\ket{\phi_\mu}\!\bra{\phi_\nu}T_{\nu i} = \sum\lgi{\mu,\nu}\sum\lgi{\text{occ }i}T_{\nu i}\tilde{T}^\dagger_{i\mu}\ket{\phi_\mu}\!\bra{\phi_\nu} = \sum\lgi{\mu,\nu}R_{\mu\nu}\ket{\phi_\mu}\!\bra{\phi_\nu}

\end{equation}

\begin{equation}

R = T\tilde{T}^\dagger = T(T\sigma^{-1})^\dagger = T\sigma^{-1}T^\dagger

\end{equation}

where $R$ is the $1$-electron density matrix.

\paragraph{}However, the use of reciprocal orbitals changes the Fock equations. This leads to the definition of projected Fock operators $\{\hat{F}^j\}_j$ defined on each block $B_j$:

\begin{equation}

\hat{F}^j\psi_i^j = \epsilon_i^j\psi_i^j,\quad \hat{F}^j = \left(\hat{\mathds{1}}-\hat{\rho}+\hat{\rho}^{j\dagger}\right)\hat{F}\left(\hat{\mathds{1}}-\hat{\rho}+\hat{\rho}^{j}\right)

\end{equation}

\begin{equation}

\hat{\rho}^{j} = \sum\lgi{i \in B_j}\ket{\tilde{\psi}_i^j}\!\bra{\psi_i^j} = \sum\lgi{\mu,\nu}R_{\mu\nu}^j\ket{\phi_\mu}\!\bra{\phi_\nu}

\end{equation}

where $\mathds{1}$ is the identity matrix, $\hat{F}$ is the usual Fock operator, $\hat{F}^j$ is the projected Fock operator on block $B_j$, $\hat{\rho}$ is the usual density operator, $\hat{\rho}^j$ is the fragment density operator defined by using only localized orbitals of block $B_j$, and $R^j$ is the associated fragment density matrix.

\subsubsection{Updated SCF loop}

\paragraph{}Instead of solving a single Fock eigenvalue equation for the whole system, one solves $n$ independent projected Fock equations ($n$ is the number of blocks in the system) in the ALMO formalism. The projected Fock equation is solved independently for each block, leading to a new SCF scheme:

The density $\hat{\rho}$ is used to construct the usual Fock operator $\hat{F}$ ; this Fock operator is projected using the fragment densities ; projected Fock operators $\hat{F}^j$ are diagonalized, providing with new orbitals $T$ ; a new density $\hat{\rho}$, and new fragment densities $\hat{\rho}^j$ are computed from the new orbitals ; and so forth until convergence.% (see figure ???).

%TODO include updated SCF loop

\subsubsection{Discussion}

\paragraph{}The main differences introduced by the ALMO formalism can be summed up in the table \ref{ALMO_comp}:

\begin{table}[h!]

\begin{center}

\begin{tabular}{|c|c|c|}

\hline

operator & orthogonal MO & ALMO\\

\hline

\begin{tabular}{@{}c@{}}one-electron \\ density\end{tabular} & $\hat{\rho} = \sum\lgi{\text{occ }i} \ket{\psi_i}\!\bra{\psi_i}$ & $\hat{\rho} = \displaystyle\sum\lgi{\text{block }j} \hat{\rho}^j = \displaystyle\sum\lgi{\text{block }j \\\text{occ }i \in B_j} \ket{\psi_{i}^{j}}\sigma^{-1}\!\bra{\psi_{i}^{j}}$\Tstrut\\

\hline

Fock & $\hat{F}$ & $\hat{F}^j = (\hat{\mathds{1}}-\hat{\rho}+\hat{\rho}^{j\dagger})\hat{F}(\hat{\mathds{1}}-\hat{\rho}+\hat{\rho}^j)$\Tstrut\\

\hline

\end{tabular}

\end{center}

%\vspace{-0.6cm}

\caption{\small Synoptic table comparing formulae with and without ALMO formalism. \label{ALMO_comp}}

\end{table}

\paragraph{}It is important to notice that density operators have a crucial role in the ALMO formalism. Indeed, each block "sees" the other blocks through these density operators $\hat{\rho}^j$, allowing to treat independently each block (for a given set of density operators).

Consequently, in order to adapt the ALMO formalism to mixed-state theory, the major and only challenge is to adapt the density matrices $\{R^j\}_j$.

%Besides, in the equation, those density matrices are constructed from only two elements: the localized orbitals and the inverse overlap matrix.

\paragraph{}Each projected Fock equation is solved independently for each block, diagonalizing each projected Fock operator separately, instead of diagonalizing the Fock operator for the whole system. Consequently the resolution of the Fock equation in the ALMO formalism can be achieved through linear scaling when dealing with large systems (assuming blocks of similar sizes). 

Therefore, the computational bottleneck becomes the reconstruction of the global Fock operator $\hat{F}$ for large system. This step is particularly fast with the CP2K modelling tool\footnote{Thanks to the use of plane waves as auxiliary basis, with the mixed Gaussian and plane waves approach (GPW) of CP2K\cite{vandevondele_quickstep:_2005}.}, explaining the choice of CP2K for the implementation of ALMO.

However, the use of non-orthogonal orbitals induces more costly equations, and the need to compute at each SCF step $\sigma$, and its inverse $\sigma^{-1}$.

Nevertheless, this added complexity is partially compensated by the sparsity of $T$, which allow faster matrix multiplications by block.

\section{Rigorous approach}

\label{rigorous_section}

\subsection{General formula}

\paragraph{}Let us apply\footnote{Indeed, the general mixed-state formula is general enough to be compatible with ALMO formalism, since it basically manipulates independent instances, each one representing a given quantum state. Therefore, as long as a quantum state is defined, the general principles of mixed-state theory can be applied.} the general mixed-state formula for the computation of the mixed-state density matrix $^\delta\!R$ in the ALMO formalism:

\begin{equation}

^\delta\!R = \sum\lgi{\mathcal{S}\in\Omega}p_\mathcal{S}R_\mathcal{S} = \sum\lgi{\mathcal{S}\in\Omega}p_\mathcal{S}T_\mathcal{S}\sigma^{-1}_\mathcal{S}T_\mathcal{S}^\dagger = \sum\lgi{\mathcal{S}\in\Omega}p_\mathcal{S}T_\mathcal{S}(T_\mathcal{S}^\dagger S T_\mathcal{S})^{-1}T_\mathcal{S}^\dagger

\label{general_ALMO_formula}

\end{equation}

where $p_\mathcal{S}$ is the probability that the real system is in the quantum state $\mathcal{S}$, $R_\mathcal{S}$ is the density matrix of this quantum state (since $\mathcal{S}$ is a single quantum state, the usual ALMO formalism can be applied), $T_\mathcal{S}$ is the orbital coefficient matrix associated with the wavefunction $\Psi_\mathcal{S}$ (i.e. $T_\mathcal{S} = \left(\begin{matrix}\ket{\psi_1^\mathcal{S}}\!\!&\!\!\cdots\!\!&\!\!\ket{\psi_n^\mathcal{S}}\end{matrix}\right)$, where $\psi_i^\mathcal{S}$ is one of the (doubly occupied) orbitals selected to construct $\Psi_\mathcal{S}$), $\sigma^{-1}_\mathcal{S}$ is the overlap matrix of the quantum state $\mathcal{S}$, and $S$ is the basis set function overlap matrix common for all quantum states.

\paragraph{}This formula is particularly unpractical, since it even doesn't provide with a single general $^\delta T$ matrix, and a general $^\delta\sigma$ overlap matrix. It is however possible to slightly rewrite the equations, in order to obtain a more elegant formulation, of the form:

\begin{equation}

^\delta\!R = \ \!^\delta T\ \!^\delta\sigma^{-1}\ \!^\delta T^\dagger

\end{equation}

\subsection{Elegant reformulation}

\paragraph{}Let us consider a general orbital coefficient matrix $^\delta T$ containing all the localized (doubly occupied) orbitals used to construct every $T_\mathcal{S}$. In other words, $^\delta T = \left(\begin{matrix}\ket{\psi_1}\!\!&\!\!\cdots\!\!&\!\!\ket{\psi_{n+k}}\end{matrix}\right)$, where $\psi_i$ is such that a quantum state $\mathcal{S}\in \Omega$ satisfying $\psi_i\in \Psi_\mathcal{S}$ exists. Therefore, one can construct any $T_\mathcal{S}$ from $^\delta T$:

\begin{equation}

T_\mathcal{S} = \ \!^\delta T\Delta_\mathcal{S}= \ \!^\delta T\left(\begin{matrix}

\delta_1 & 0 & \cdots & 0\\

0 & \delta_2 & \cdots & 0\\

\vdots & \vdots & \ddots & \vdots\\

0 & 0 & \cdots & \delta_{n+k}

\end{matrix}\right),\quad \delta_i = \begin{cases}

1 & \text{if } \psi_i\in\Psi_\mathcal{S},\\

0 & \text{otherwise}.\end{cases}

\label{T_proj}

\end{equation}

where $\Delta_\mathcal{S}$ can be seen as a rescaling matrix.

\paragraph{}By injecting equation (\ref{T_proj}) in equation (\ref{general_ALMO_formula}), we obtain the following more elegant reformulation:

\begin{equation}

^\delta\!R = \sum\lgi{\mathcal{S}\in\Omega}p_\mathcal{S}R_\mathcal{S} = \sum\lgi{\mathcal{S}\in\Omega}p_\mathcal{S}T_\mathcal{S}\sigma^{-1}_\mathcal{S}T_\mathcal{S}^\dagger = \sum\lgi{\mathcal{S}\in\Omega}p_\mathcal{S}\ \!^\delta T\Delta_\mathcal{S}\sigma^{-1}_\mathcal{S}\Delta_\mathcal{S}\ \!^\delta T^\dagger = \ \!^\delta T\left(\sum\lgi{\mathcal{S}\in\Omega}p_\mathcal{S}\Delta_\mathcal{S}\sigma^{-1}_\mathcal{S}\Delta_\mathcal{S}\right)\!^\delta T^\dagger = \ \!^\delta T\ \!^\delta\sigma^{-1}\ \!^\delta T^\dagger

\end{equation}

where $\ \!^\delta\sigma^{-1} = \sum\lgi{\mathcal{S}\in\Omega}p_\mathcal{S}\Delta_\mathcal{S}\sigma^{-1}_\mathcal{S}\Delta_\mathcal{S}$. In this formulation, it is clear that all the complexity of the equations is still here, but simply hidden in the inverse overlap matrix $\ \!^\delta\sigma^{-1}$.

\subsection{Computational considerations}

\paragraph{}This exact method is the rigorous unification of mixed-state theory and ALMO formalism. Let us evaluate the computational cost of this approach.

\paragraph{Combinatorial explosion:}The number of quantum states to consider is theoretically infinite, since one must in theory include an infinite number of additional orbitals. However, it is common to include a small amount\footnote{In this study, the number of added orbitals $k$ is arbitrarily taken as the number of atoms in the system.} $k$ of additional orbitals to the initial $n$ orbitals (since the occupation number decreases rapidly as the orbital energy increases). Combinatorial considerations show that the number $|\Omega|$ of quantum states in $\Omega$ is:

\begin{equation}

|\Omega| = \left(\begin{matrix}n+k\\n\end{matrix}\right) = C^n_{n+k} = \frac{(n+k)!}{n!\,k!} \approx_{k=n} \frac{2^{2n}}{\sqrt{\pi n}}

\end{equation}

Therefore, in order to compute $\ \!^\delta\sigma^{-1}$, one has to compute $|\Omega|$ terms, which is exponential in $n$ in the worst case (when $k=n$). For each term, one must construct $\Delta_\mathcal{S}$, compute the energy of $\mathcal{S}$ and determine $p_\mathcal{S}$, perform $6$ sparse matrix multiplications, and $1$ matrix inversion.

Clearly, the most costly operation is the matrix inversion, so let us just focus on that. In the current implementation of ALMO, this matrix inversion is performed by the Hotelling-Bodewig\cite{hotelling_analysis_1933} algorithm\footnote{The Hotelling-Bodewig algorithm\cite{hotelling_analysis_1933}\cite{haghani_new_2014} is based on the recurrence relation: $A_{n+1} = A_n(2\times\mathds{1}-BA_n)$, where $B^{-1} = A$.}, an iterative procedure that can use an initial guess of the inverse to converge faster. Currently, the initial guess for the inverse overlap matrix is the $\sigma^{-1}$ from the previous SCF step. Therefore, in order to avoid longer inversions, one must store $|\Omega|$ overlap matrices (one for each quantum state).

As a conclusion, the time complexity (normalized cost in term of execution time) of this method is already prohibitive for most simulations. Besides, to avoid an even bigger time complexity, the space complexity (normalized cost in term of RAM consumption) of this method must be exponential in the worst cases. Therefore, this method cannot be used in real applications for practical reasons.

\paragraph{Parallelization:}Even if the computational cost of $\ \!^\delta\sigma^{-1}$ is \textit{a priori} prohibitive, each term can be computed independently from the others. It is therefore easy to design a parallel version of this method, allowing to distribute the workload over a large number of processing units (through GPU computing for example).

However, for most real applications, the computational cost of this method would still be prohibitive.

\section{Occupation-state dependency}

\subsection{Problematic}

\paragraph{}The rigorous approach formulation is unsuitable for practical applications. However, the same formula is used in usual mixed-state theory. Only that in the case of pure mixed-state theory, a mathematical trick was developed in order to avoid the combinatorial exploration.

\paragraph{}Basically, the rescaling trick exploits the fact that each orbital's impact on the electronic density is common to all quantum states constructed with this occupied orbital. Therefore, the density corresponding to this orbital can be computed only once, and weighted by the probability that the real system is in a quantum state containing this orbital (or equivalently, the probability that this orbital is occupied in the real system). In other words, the orbitals and their impact on the density are basically reused from one state to the other, avoiding the combinatorial approach.

\paragraph{}We would like to apply such a simplification to our rigorous approach in order to make it suitable for practical applications. For such a simplification to be applicable on our formula, we need to be able to reuse the orthonormalized orbitals computed from one quantum state to the other.

\paragraph{}Unfortunately, it can be mathematically proven that there is a strong occupation-state dependency of the orthonormalized orbitals in the ALMO formalism.

\paragraph{}To be more precise, let us consider a quantum state $\mathcal{S}$, with associated wave-function $\Psi_\mathcal{S}$, and two arbitrary orbitals $\psi_1$ and $\psi_2$ such that $\psi_1\in\Psi_\mathcal{S}$ and $\psi_2\notin\Psi_\mathcal{S}$.

The orthonormalized orbital $\psi_1^L$ in the quantum state $\mathcal{S}$ does not change by adding the orbital $\psi_2$ to $\Psi_\mathcal{S}$, if and only if $\psi_1$ and $\psi_2$ does not overlap with each other, and there exists no $\psi_i$ in that quantum state $\mathcal{S}$ such that both $\psi_1$ and $\psi_2$ overlap with $\psi_i$:

\begin{equation}

\ket{\psi_1^L} = \ket{{\psi_1^L}^\prime} \Leftrightarrow \braket{\psi_1}{\psi_2} = 0\ \wedge\ \left(\nexists \psi_i\in\Psi_\mathcal{S},\ \braket{\psi_1}{\psi_i} \neq 0\ \wedge\ \braket{\psi_2}{\psi_i} \neq 0\right)

\end{equation}

\paragraph{}Therefore, the occupation-state dependency of the orthonormalized orbitals is always present, except when the added or removed orbitals does not overlap with the rest of the system.

\subsection{Chemical interpretation}

\paragraph{}The average contribution of an orbital to the density cannot be evaluated, since its orthonormalized orbital varies from one quantum state to the other. Let us try to understand these mathematical results based on chemical intuition.

\paragraph{Dependancy between blocks:}Let us imagine a $2$-blocks system. The orthonormalized orbitals of the first block $B_1$ are $\{\psi_i^L\}_i$, and the density contribution of $B_1$ is $\hat{\rho_1}$.

The second block $B_2$ can be seen here as an environment\footnote{An environment can be simply seen as a modification of the potential $v(\mathbf{r})$ on which the Hamiltonian of the system $\hat{H}$ depends.} for the first block $B_1$ (unless $B_2$ and $B_1$ are very far apart, i.e. non overlapping).

Consequently, a modification of the block $B_2$ (like a modification of the occupation-state of $B_2$) can be seen as a modification of the environment for $B_1$. But the orthonormalized orbitals (and their corresponding density) depends on the environment, since the environment has a direct impact on the Hamiltonian of the system.

Therefore, it becomes clear now how the density (and the orthonormal orbitals) of a block can be affected by the occupation-state of other blocks.

\paragraph{Dependancy within blocks:}However, the mathematical derivations provide with an even more powerful statement: Even within a block, the occupation-state of an orbital $\psi_1$ can affect the orthonormalization of $\psi_2$, through the intermediate of other blocks.

Indeed, the modification of a block $B_1$, can be seen as a modification of environment for the neighboring blocks $\{\mathrm{neigh}(B_1)\}$ of $B_1$. Therefore, the blocks $\{\mathrm{neigh}(B_1)\}$ are changed, impacting their neighbors $\{\mathrm{neigh}(\{\mathrm{neigh}(B_1)\})\}$, and so forth, including $B_1$ itself.

\label{subway_analogy}

\paragraph{}As a conclusion, since the orbitals overlap (i.e. interact) one with an other, when an orbital is modified the whole system has to re-adapt.% One can picture an expressive analogy with people in the subway: the position of a person depends on the disposition and presence of other people, just as the spatial probability distribution of an (orthonormalized) orbital depends on the presence of other orbitals.

\subsection{Conclusion}

\paragraph{}The orthonormalized orbitals (and their impact on the density) cannot be reused from one quantum state to the other, since they depend on the occupation-state of the other orbitals (through the $\sigma^{-1}$ term). In other words, the localized orbitals are treated (in order to evaluate their impact on the density) differently in each quantum state, since the inverse overlap matrix $\sigma^{-1}_\mathcal{S}$ depends on the quantum state $\mathcal{S}$.

This occupation-state dependency of the ALMO formalism prevents therefore the use of the rescaling trick to avoid a combinatorial approach (as achieved in pure mixed-state theory).

\paragraph{}In fact, the rescaling trick is only applicable when a non interacting sets of orbitals (i.e. non interacting system) is considered. But in such a case, it is always faster to perform an independent simulation of this isolated sub-system.

Therefore, one cannot usually apply this mathematical trick, and even when it is possible to apply this trick, it means that one should not be running such a simulation in the first place!

\section{Mean-field approximation}

\paragraph{}Since the rigorous exact approach cannot be simplified by the usual reductions, we propose here a mean-field approximation in order to avoid the combinatorial enumeration.

\subsection{Chemical intuition}

\paragraph{}For each quantum-state $\mathcal{S}$, a set of orbitals is selected, generating interactions and overlap between the orbitals. These interactions are consequently specific to the quantum-state $\mathcal{S}$ (these interactions are summarized in the overlap matrix $\sigma_\mathcal{S}$). Therefore, for each quantum state $\mathcal{S}$ the system must re-adapt its orbitals because of these specific interactions $\sigma_\mathcal{S}$. The combinatorial problem is then coming from the fact that each quantum-state $\mathcal{S}$ has its own interactions $\sigma_\mathcal{S}$.

\paragraph{}However, one could consider average interactions to apply to all quantum-states independently. In other words, one could simply formulate a mean-field approximation: instead of treating each quantum-state with its specific interactions, we compute an average interaction field (over all quantum-states) that we apply to every quantum-state.

\subsection{Mathematical derivation}

\paragraph{}This mean-field approximation is about using average interactions between orbitals. These interactions (i.e. overlap) are represented by the orbital overlap matrix. Therefore, averaging the interactions between orbitals is equivalent to obtaining an average overlap matrix.

\paragraph{}In order to average the interaction between $\psi_i$ and $\psi_j$, we decided\footnote{Although more sophisticated weights could also be considered, this rescaling was chosen due to its similarities with the rescaling trick.} to simply rescale their overlap by $\sqrt{p_i}\sqrt{p_j}$. This rescaling provides with an overlap matrix $^\epsilon\sigma$:

\begin{equation}

^\epsilon\sigma_{ij}=\begin{cases}

\braket{\psi_i}{\psi_i} & \text{if } i=j,\\

\braket{\psi_i}{\psi_j}\sqrt{p_i}\sqrt{p_j} & \text{otherwise}.\end{cases}

\end{equation}

\paragraph{}This mean-field orbital overlap matrix $^\epsilon\sigma$ can then be applied in every quantum state. Since the same $^\epsilon\sigma$ is used for all quantum state (i.e. $\forall \mathcal{S},\ ^\epsilon\sigma_\mathcal{S} =\ \!^\epsilon\sigma$), then the rescaling trick can be used with this approximation. Therefore, the density matrix $^\epsilon\!R$ can be written:

\begin{equation}

^\epsilon\!R = T^\prime\ \!^\epsilon\sigma^{-1}T^{\prime\dagger}

\end{equation}

where $T^\prime$ is the usual rescaled orbital coefficient matrix (i.e. the orbital coefficient matrix corresponding to the rescaled orbitals $\psi_i^\prime = \psi_i\sqrt{p_i}$).

%\paragraph{}Mathematically, this approximation is therefore equivalent to the following approximation:

%\begin{equation}

%\sum\lgi{\mathcal{S}\in\Omega}p_\mathcal{S}\sigma_\mathcal{S}^{-1} \approx \left(\sum\lgi{\mathcal{S}\in\Omega}(p_\mathcal{S}\sigma_\mathcal{S}+(1-p_\mathcal{S})\mathds{1})\right)^{-1}

%\end{equation}

%which is therefore exact when $\Omega$ contains only a single non-negligible quantum state (i.e. infinitesimal smearing), or if $\sigma_\mathcal{S} = \mathds{1}$ for each quantum-state $\mathcal{S}$ (i.e. $1$-block systems).

\paragraph{}It is important to note here that under this approximation, the global form of the ALMO equations is not changed. Simply, the following changes are made:

\begin{itemize}

\item $\sigma$ is replaced by $^\epsilon\sigma$,

\item and $T$ is replaced by $T^\prime$, the rescaled orbital coefficient matrix.

\end{itemize}

Indeed, updating the fragment density matrices $\{R^j\}_j$ is enough to adapt the whole ALMO formalism, and those matrices depend only on $\sigma$ and $T$. Therefore, updating $\sigma$ and $T$ is enough to adapt the whole ALMO formalism.

\subsection{New mathematical trick derivation}

\paragraph{}In pure mixed-state theory, a mathematical trick (i.e. the rescaling trick) was developed so that instead of changing the equations, one could simply change the orbitals, and perform the computations as in the usual formalism.

\paragraph{}In analogy with this rescaling trick, it is possible to develop a modification of the orbitals (i.e. a mathematical trick), so that with these new orbitals, the usual ALMO formalism can be applied to compute our mean-field approximation of mixed-state ALMO (called S-ALMO for smearing ALMO).

\paragraph{}The new orbitals $\psi_i^s$ to consider in this mathematical trick have rescaled interactions (as in the rescaling trick), except with themselves. This is the crucial difference between the rescaled orbitals $\psi_i^\prime$ and these new orbitals $\psi_i^s$:

\begin{equation}

\forall \mu\ \braket{\psi_i^s}{\phi_\mu} = \sqrt{p_i}\braket{\psi_i}{\phi_\mu}\quad;\quad\forall j\neq i\ \braket{\psi_i^s}{\psi_j^s} = \sqrt{p_i}\sqrt{p_j}\braket{\psi_i}{\psi_j}\quad;\ \text{but }\braket{\psi_i^s}{\psi_i^s}=\braket{\psi_i}{\psi_i}

\end{equation}

Therefore, our mean-field approximation of mixed-state ALMO (S-ALMO) can be performed by a usual ALMO computation, with the orbitals $\psi_i$ replaced by these new orbitals $\psi_i^s$.

\subsection{Chemical interpretation: selfish orbitals}

\paragraph{}These new orbitals $\psi_i^s$ define a new concept: selfish orbitals. These orbitals have rescaled interactions with everything except themselves, hence the name.%: selfish orbitals.

Unlike rescaled orbitals, selfish orbitals cannot be considered shrunk. They just somehow interact less with their environment. Hence the name.

\paragraph{}The concept of selfish orbitals is a more subtle concept than rescaled orbitals, since selfish orbitals can replace rescaled orbitals for the usual rescaling trick (but the converse is not true).

Besides, the normalization is preserved with selfish orbitals, unlike with rescaled orbitals.

\paragraph{}In the usual mixed-state theory, working with non normalized orbitals is not a major problem in itself. But in the ALMO formalism, the orbitals are forced to be normalized (through the $\sigma$ inversion), causing the rescaling trick to be inapplicable.

%\section{Computational considerations}

\section{Testing and results}

\paragraph{}S-ALMO was successfully tested on physically meaningful and reasonable simulations. Even if those simulations were performed solely for test purposes, the main results are summarized here.

\subsection{Implementation}

\paragraph{}Our newly developed mean-field approximation of mixed-state ALMO (called S-ALMO) has been successfully implemented in the CP2K simulation tool software, in its version 5.0 (Development Version)

\paragraph{}Thanks to the selfish orbital trick, S-ALMO has been integrated to the pre-existing FORTRAN subroutines developed for usual ALMO formalism\cite{khaliullin_efficient_2013}\cite{khaliullin_efficient_2006}. This integration has been thought to have the minimum impact on the algorithmic flow of ALMO computations.

\paragraph{}Indeed, the modifications can be summarized by the addition of the following steps:

\begin{itemize}

\item During the diagonalization of the Fock operators, the orbital energies are retrieved $\{\epsilon_i\}_i$.

\item These orbital energies $\{\epsilon_i\}_i$ are converted into probabilities $\{p_i\}_i$ for each block independently, using a Fermi-Dirac smearing model.

\item The orbitals $\{\psi_i\}_i$ obtained during the diagonalization of the Fock operators are rescaled (before computing the new density matrices) $\psi_i^\prime = \psi_i\sqrt{p_i}$.

\item The overlap matrix $\sigma$ is computed using the rescaled orbitals, then the diagonal is set to $1$ (to obtain the effect of selfish orbitals).

\item The rest of the algorithm and equations (inversion of $\sigma$, computation of new Fock operators, convergence check, \ldots) are left unchanged with these new $T$ and $\sigma$.

\end{itemize}

\paragraph{}Finally, the computational cost of this extension is negligible thanks to the use of selfish orbitals. And the additional storage used is simply limited to an array containing the orbital energies.

\subsection{Quality evaluation}

\paragraph{}S-ALMO is, before anything, an approximation. Even if a mean-field approximation is common and often reasonable (in fact, the DFT SCF formulation itself is a mean-field approximation), it is crucial to design methods to estimate the error introduced by this approximation.

\paragraph{}As shown before, S-ALMO is exact in the case of infinitesimal smearing, and in the case of $1$-block systems (i.e. orthogonal orbitals). Thus, one needs to focus on multiple block systems with important smearing. Heterogeneous catalysis simulations could therefore provide with interesting tests.

Unfortunately, ALMO simulations of metallic systems was not possible before the development of S-ALMO. It is therefore particularly problematic to find a reference against which our results could be compared.

\subsubsection{Density recovery}

\paragraph{}Even if no proper references are available, a quality estimator can be conceived: the number of electrons recovered. Indeed, since we are constructing the density of a N-electrons system, then this density must account for exactly N electrons (i.e. the N electrons have to be recovered in the density).

\paragraph{}The number of electrons recovered $N_{calc}$ can be computed by integrating the density over the whole space. In matrix form, this calculation becomes:

\begin{equation}

N_{calc} = \sum\lgi{\mu,\nu}R_{\mu\nu}S_{\mu\nu} = \sum\lgi{\mu,\nu}R_{\mu\nu}S_{\nu\mu} = \sum\lgi{\mu}(RS)_{\mu\mu} = \mathrm{Tr}(RS)

\end{equation}

\paragraph{}Since this estimator needs to be computed only once, at the end of a simulation, it can then be computed systematically with negligible additional cost.

%\paragraph{}Therefore, this estimator can be easily computed, with a quite low computational cost (since this estimator needs to be computed only once, at the end of a simulation). This estimator can then be computed systematically with negligible additional cost.

\paragraph{Counterexample evaluation:} As an example, the number of electron pairs recovered from the 2-electron pairs counterexample studied section \ref{counterexample} is $N_{calc} \approx 1.83 \neq 2$. Therefore, $8.33\%$ of the total density was lost.

%However, the non-reliability of this estimator can be mathematically proven.

\paragraph{}However, this estimator provides with a test that is only a necessary, not sufficient condition. Thus, a total recovery of the density would not necessarily mean that the S-ALMO result is exact\ldots

Hence, a critical assessment of the physical meaning of the results is required.

%\subsubsection{Infinitesimal smearing limit prediction}

%\paragraph{}The energy $E_{avg}$ of a mixed-state is simply the average energy of all quantum states, according to mixed-state theory (i.e. $E_{avg} = \sum\lgi{\mathcal{S}\in\Omega}p_\mathcal{S}E_\mathcal{S}$)

%Can we really predict $E_S$? Because we can indeed predict easily the energy of each block, but the total energy? I doubt so...

\subsection{Charge transfer analysis}

\paragraph{}The newly developed S-ALMO was applied, for test purposes, to charge transfer analysis of adsorption at a metallic surface. Therefore, the charge transfer analyzed here are between the metallic surface and the molecule adsorbed.

Several systems were studied, but only the two extreme systems will be presented here.

\paragraph{}For the charge transfer energy to be computed, one must compare the energy obtained when charge transfer is allowed (usual SCF DFT single energy-point calculation), and forbidden (S-ALMO single energy-point calculation).

\paragraph{}Consequently, each system will clearly be composed of $2$ blocks: a metallic block containing the metal surface, and an organic block containing the adsorbed molecule. It is important to note that the definition of such blocks is meaningful since only adsorption is considered. Therefore, no real chemical bond must be cut (in such a case, the distribution of the bonding electron pair would be highly problematic\ldots).

\paragraph{}However, the ALMO formalism does not only suppress charge transfer between block, but it also eliminates the Basis Set Superposition Error\footnote{This error happens (in finite basis set computations) when a fragment uses the basis set functions of an other fragment for its description (which is by definition forbidden in the ALMO formalism). Basically, it leads to some regions being better described than other, leading to an energy mismatch responsible for the BSSE.} (BSSE) between both blocks.

Therefore, S-ALMO simulations should in fact be compared to BSSE-free simulations.

\subsubsection{Results}

\paragraph{Description of the simulations:}The adsorption of water and carbon monoxide at a platinum surface (Pt(111) surface, on top position, as displayed figure \ref{systems}) was studied.

The adsorbed structures were optimized with VASP, a plane-wave based electronic structure code.

For both S-ALMO and BSSE-free SCF DFT simulations: molecular orbitals were represented by a double-$\zeta$ Gaussian basis set with one set of polarization functions (i.e. DZVP). A cutoff of $400$ Ry was used to describe the electron density. The exchange-correlation (XC) energy was approximated with the PBE functional. The Brillouin zone was sampled at the $\Gamma$-point. GTH pseudo-potentials were used to describe the interactions between the valence electrons and the ionic cores, and the electronic smearing was approximated by a Fermi-Dirac distribution at $300$ K.

%The parameters used for the S-ALMO, and BSSE-free SCF DFT simulations are: PBE functional with VdW correction, DZVP MolOpt basis set, GTH-PBE pseudopotential, EPS-SCF $1\times10^{-7}$

%Description of the formula used

\paragraph{}An energy decomposition analysis (of the interactions between the metallic surface and the adsorbed molecule) was performed, separating the total interaction energy $E_{tot}$ into $3$ specific contributing terms: charge transfer ($E_{CT}$), polarization ($E_{pol}$), and finally, electrostatic and exchange repulsion ($E_{elec\text{-}xc}$)\cite{mo_energy_2000}:

% = E_{SCF} - E_{S\text{-}ALMO} - E_{BSSE}

\begin{equation}

E_{tot} = E_{CT} + E_{pol} + E_{elec\text{-}xc}

\end{equation}

\paragraph{}These terms can be expressed as:

\begin{subequations}

\begin{align}

E_{CT} & = E_{BSSE\text{-}free} - E_{S\text{-}ALMO}\\

E_{pol} & = E_{S\text{-}ALMO} - E_{frozen}\\

E_{elec\text{-}xc} & = E_{frozen} - E_{sum}

\end{align}

\end{subequations}

where $E_{BSSE\text{-}free}$ denotes the energy computed by BSSE-corrected SCF DFT calculation. $E_{S\text{-}ALMO}$ denotes the energy computed with the S-ALMO formalism. Let $B_1^{only}$ (resp. $B_2^{only}$) be a system composed of the block $B_1$ (resp. $B_2$) alone. $E_{sum} = E_{B_1^{only}}+E_{B_2^{only}}$ denotes the sum of the energy of $B_1^{only}$ and $B_2^{only}$. Finally, $E_{frozen}$ is the energy of a system where both blocks are brought together, while keeping their densities (i.e. $\hat{\rho}_{frozen} = \hat{\rho}_{B_1^{only}}+\hat{\rho}_{B_2^{only}}$).

\paragraph{}The results are summarized in the following table:

\begin{table}[h!]

\begin{center}

\begin{tabular}{|c|c|c|}

\hline

\begin{tabular}{@{}c@{}}energy contribution \\ (kcal/mol)\end{tabular} & CO@Pt(111) & H$_2$O@Pt(111)\\

\hline

$E_{CT}$ & -121.9 & -31.4\\

\hline

$E_{pol}$ & -71.4 & -11.1\\

\hline

$E_{elec\text{-}xc}$ & 151.7 & 33.4\\

\hline

$E_{tot}$ & -41.6 & -9.0\\

\hline

\hline

density recovery & \begin{tabular}{@{}c@{}}$2.6*10^{-4}\%$ \\ of density lost\end{tabular} & \begin{tabular}{@{}c@{}}$2.6*10^{-5}\%$ \\ of density lost\end{tabular}\\

\hline

\end{tabular}

\end{center}

\caption{\small Energy decomposition analysis and comparison of the adsorption of carbon monoxide and water at a platinum (Pt(111)) surface.}

\end{table}

\begin{figure}[h!]

\centering

\begin{subfigure}[b]{0.45\linewidth}

\includegraphics[width=\linewidth]{H2O_top_Pt.png}

\vspace{-0.2cm}

\caption{\centering\small{H2O@Pt: E$_{\mathrm{CT}} = -31.4$ kcal/mol. $2.6*10^{-4}\%$ of density lost.}}

\label{H2O_Pt}

\end{subfigure}

\hfill

\begin{subfigure}[b]{0.45\linewidth}

\centering

\includegraphics[width=\linewidth]{CO_top_Pt.png}

\vspace{-0.2cm}

\caption{\centering\small{CO@Pt: E$_{\mathrm{CT}} = -121.9$ kcal/mol. $2.6*10^{-5}\%$ of density lost.}}

\label{CO_Pt}

\end{subfigure}

\caption{Representation of the two systems described in this study, whose energy decomposition analysis was performed using S-ALMO and BSSE-free simulations.\label{systems}}

\end{figure}

\subsubsection{Discussion}

\paragraph{}CO is generally a much stronger ligand than H$_2$O. This is partially due to the electron donation (i.e. forward donation) and $\pi$-backdonation synergistic effects, that are particularly important in the case of carbon monoxide. These effects are basically based on charge transfer mechanisms.% forward donation

\paragraph{}Therefore, one can expect the energy contribution of charge transfer to be much higher in the adsorption of CO, than for H$_2$O. This expectation is indeed coherent with the empirical results, based on S-ALMO simulations.

\paragraph{}However, a lot more results can be extracted from this simple energy decomposition analysis (e.g. the relative contribution of charge transfer to adsorption is relatively lower with CO than H$_2$O, \ldots), that might shed new light on adsorption mechanisms, for example.

%An other important point is the polarization energy:

\paragraph{}Finally, S-ALMO allows us to assess which fraction of the energy and which adsorption behavior could be obtained in the best possible force field that excludes charge-transfer and thus covalent bond formation\cite{raimondi_ab_2001}.

\section{Conclusion and perspectives}

\subsection{Conclusion}

\paragraph{}The extension of ALMO formalism to metallic systems (and more generally, any system that requires electronic smearing to be properly described) was performed in this work. This extension has been possible through the unification of mixed-state theory and ALMO formalism.

\paragraph{}Unfortunately, the rigorous unification is highly unpractical, as its computational cost is prohibitive due to a combinatorial enumeration.

Besides, usual simplifications (like the rescaling trick) cannot be applied because of an occupation-state dependency introduced by the ALMO formalism.

\paragraph{}However, a reasonable mean-field approximation is proposed by averaging the overlap (i.e. interactions) between orbitals.

This mean-field approximation of mixed-state ALMO is called S-ALMO (for smearing ALMO), and is exact in the case of infinitesimal smearing (i.e. pure ALMO formalism) and/or $1$-block systems (i.e. pure mixed-state theory).

\paragraph{}In order to implement our newly developed S-ALMO theory  without changing the usual ALMO equations, a new mathematical trick was developed in analogy with the rescaling trick.

This new mathematical trick is based on a new concept: selfish orbitals, whose interactions are rescaled (like rescaled orbitals), except with themselves.

\paragraph{}Thanks to the use of selfish orbitals, S-ALMO has successfully been implemented in the open source simulation software CP2K, as a patch for the pre-existing ALMO subroutines.

%Besides, this new mathematical trick allows the use of S-ALMO formalism by adding only a negligible computational cost, compared with usual ALMO formalism.

S-ALMO adds only a negligible computational cost, compared with the usual ALMO formalism, making it applicable to large and diverse systems.

\paragraph{}An estimator was designed to help the user assess the quality of the approximation. Finally, S-ALMO has been successfully tested on charge transfer analyses, showing promising results.

\subsection{Perspectives}

\paragraph{}S-ALMO is a tool with numerous potential applications in various fields of chemistry. A few of them are mentioned here:

\paragraph{}First of all, the computational cost of S-ALMO is reduced compared to traditional DFT SCF calculations. This could, for example, make the use of explicit solvent far more affordable for simulations including a metallic system in solution (e.g. heterogeneous catalysis).

\paragraph{}In S-ALMO, the charges can be forced to remain localized on predefined domains. Therefore, with S-ALMO, electrocatalytic simulations would suffer far less from the DFT tendency to overestimate electron delocalization. And more generally, S-ALMO will provide a computational benefit for any system where both electrolytes and a metallic part are present, and at the same time improving the physical soundness of the description of the system.

\paragraph{}Finally, S-ALMO provides a proper and convenient (i.e. easy to use) charge transfer analysis tool, and therefore an energy decomposition analysis tool.

The application of this energy decomposition analysis tool could shine a new light on the impact of charge transfer/polarization in catalytic processes (or virtually any system containing a metallic part).

This is especially pertinent for unraveling the origin of catalytic properties of alloys, where much is yet to be discovered. S-ALMO could therefore even provide a first step toward the rational design of alloy catalysts, a major challenge for green chemistry development.

\subsection{References}

The class makes various changes to the way that references are

handled.  The class loads \textsf{natbib}, and also the

appropriate bibliography style.  References can be made using

the normal method; the citation should be placed before any

punctuation, as the class will move it if using a superscript

citation style

\cite{Mena2000,Abernethy2003,Friedman-Hill2003,EuropeanCommission2008}.

The use of \textsf{natbib} allows the use of the various citation

commands of that package: \citeauthor{Abernethy2003} have shown

something, in \citeyear{Cotton1999}, or as given by

Ref.~\citenum{Mena2000}.  Long lists of authors will be

automatically truncated in most article formats, but not in

supplementary information or reviews \cite{Pople2003}. If you

encounter problems with the citation macros, please check that

your copy of \textsf{natbib} is up to date. The demonstration

database file \texttt{achemso-demo.bib} shows how to complete

entries correctly. Notice that ``\latin{et al.}'' is auto-formatted

using the \texttt{\textbackslash latin} command.

Multiple citations to be combined into a list can be given as

a single citation.  This uses the \textsf{mciteplus} package

\cite{Johnson1972,*Arduengo1992,*Eisenstein2005,*Arduengo1994}.

Citations other than the first of the list should be indicated

with a star. If the \textsf{mciteplus} package is not installed,

the standard bibliography tools will still work but starred

references will be ignored. Individual references can be referred

to using \texttt{\textbackslash mciteSubRef}:

``ref.~\mciteSubRef{Eisenstein2005}''.

The class also handles notes to be added to the bibliography.  These

should be given in place in the document \bibnote{This is a note.

The text will be moved the the references section.  The title of the

section will change to ``Notes and References''.}.  As with

citations, the text should be placed before punctuation.  A note is

also generated if a citation has an optional note.  This assumes that

the whole work has already been cited: odd numbering will result if

this is not the case \cite[p.~1]{Cotton1999}.

\subsection{Floats}

New float types are automatically set up by the class file.  The

means graphics are included as follows (Scheme~\ref{sch:example}).  As

illustrated, the float is ``here'' if possible.

\begin{scheme}

  Your scheme graphic would go here: \texttt{.eps} format\\

  for \LaTeX\, or \texttt{.pdf} (or \texttt{.png}) for pdf\LaTeX\\

  \textsc{ChemDraw} files are best saved as \texttt{.eps} files:\\

  these can be scaled without loss of quality, and can be\\

  converted to \texttt{.pdf} files easily using \texttt{eps2pdf}.\\

  %\includegraphics{graphic}

  \caption{An example scheme}

  \label{sch:example}

\end{scheme}

\begin{figure}

  As well as the standard float types \texttt{table}\\

  and \texttt{figure}, the class also recognises\\

  \texttt{scheme}, \texttt{chart} and \texttt{graph}.

  \caption{An example figure}

  \label{fgr:example}

\end{figure}

Charts, figures and schemes do not necessarily have to be labelled or

captioned.  However, tables should always have a title. It is

possible to include a number and label for a graphic without any

title, using an empty argument to the \texttt{\textbackslash caption}

macro.

The use of the different floating environments is not required, but

it is intended to make document preparation easier for authors. In

general, you should place your graphics where they make logical

sense; the production process will move them if needed.

\subsection{Math(s)}

The \textsf{achemso} class does not load any particular additional

support for mathematics.  If packages such as \textsf{amsmath} are

required, they should be loaded in the preamble.  However,

the basic \LaTeX\ math(s) input should work correctly without

this.  Some inline material \( y = mx + c \) or $ 1 + 1 = 2 $

followed by some display. \[ A = \pi r^2 \]

It is possible to label equations in the usual way (Eq.~\ref{eqn:example}).

\begin{equation}

  \frac{\mathrm{d}}{\mathrm{d}x} \, r^2 = 2r \label{eqn:example}

\end{equation}

This can also be used to have equations containing graphical

content. To align the equation number with the middle of the graphic,

rather than the bottom, a minipage may be used.

\begin{equation}

  \begin{minipage}[c]{0.80\linewidth}

    \centering

    As illustrated here, the width of \\

    the minipage needs to allow some  \\

    space for the number to fit in to.

    %\includegraphics{graphic}

  \end{minipage}

  \label{eqn:graphic}

\end{equation}

\section{Experimental}

The usual experimental details should appear here.  This could

include a table, which can be referenced as Table~\ref{tbl:example}.

Notice that the caption is positioned at the top of the table.

\begin{table}

  \caption{An example table}

  \label{tbl:example}

  \begin{tabular}{ll}

    \hline

    Header one  & Header two  \\

    \hline

    Entry one   & Entry two   \\

    Entry three & Entry four  \\

    Entry five  & Entry five  \\

    Entry seven & Entry eight \\

    \hline

  \end{tabular}

\end{table}

Adding notes to tables can be complicated.  Perhaps the easiest

method is to generate these using the basic

\texttt{\textbackslash textsuperscript} and

\texttt{\textbackslash emph} macros, as illustrated (Table~\ref{tbl:notes}).

\begin{table}

  \caption{A table with notes}

  \label{tbl:notes}

  \begin{tabular}{ll}

    \hline

    Header one                            & Header two \\

    \hline

    Entry one\textsuperscript{\emph{a}}   & Entry two  \\

    Entry three\textsuperscript{\emph{b}} & Entry four \\

    \hline

  \end{tabular}

  \textsuperscript{\emph{a}} Some text;

  \textsuperscript{\emph{b}} Some more text.

\end{table}

The example file also loads the optional \textsf{mhchem} package, so

that formulas are easy to input: \texttt{\textbackslash ce\{H2SO4\}}

gives \ch{H2SO4}.  See the use in the bibliography file (when using

titles in the references section).

The use of new commands should be limited to simple things which will

not interfere with the production process.  For example,

\texttt{\textbackslash mycommand} has been defined in this example,

to give italic, mono-spaced text: \mycommand{some text}.

\section{Extra information when writing JACS Communications}

When producing communications for \emph{J.~Am.\ Chem.\ Soc.}, the

class will automatically lay the text out in the style of the

journal. This gives a guide to the length of text that can be

accommodated in such a publication. There are some points to bear in

mind when preparing a JACS Communication in this way.  The layout

produced here is a \emph{model} for the published result, and the

outcome should be taken as a \emph{guide} to the final length. The

spacing and sizing of graphical content is an area where there is

some flexibility in the process.  You should not worry about the

space before and after graphics, which is set to give a guide to the

published size. This is very dependant on the final published layout.

You should be able to use the same source to produce a JACS

Communication and a normal article.  For example, this demonstration

file will work with both \texttt{type=article} and

\texttt{type=communication}. Sections and any abstract are

automatically ignored, although you will get warnings to this effect.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% The "Acknowledgement" section can be given in all manuscript

%% classes.  This should be given within the "acknowledgement"

%% environment, which will make the correct section or running title.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{acknowledgement}

Please use ``The authors thank \ldots'' rather than ``The

authors would like to thank \ldots''.

The author thanks Mats Dahlgren for version one of \textsf{achemso},

and Donald Arseneau for the code taken from \textsf{cite} to move

citations after punctuation. Many users have provided feedback on the

class, which is reflected in all of the different demonstrations

shown in this document.

\end{acknowledgement}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% The same is true for Supporting Information, which should use the

%% suppinfo environment.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{suppinfo}

A listing of the contents of each file supplied as Supporting Information

should be included. For instructions on what should be included in the

Supporting Information as well as how to prepare this material for

publications, refer to the journal's Instructions for Authors.

The following files are available free of charge.

\begin{itemize}

  \item Filename: brief description

  \item Filename: brief description

\end{itemize}

\end{suppinfo}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% The appropriate \bibliography command should be placed here.

%% Notice that the class file automatically sets \bibliographystyle

%% and also names the section correctly.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibliography{achemso-demo}

\end{document}