# Value Explicit Pretraining for Goal-Based Transfer Learning

Kiran Lekkala, Henghui Bao, Sumedh Anand Sontakke, Laurent Itti
Department of Computer Science
University of Southern California
Los Angeles, CA 90089, USA
<EMAIL_ADDRESS>

###### Abstract

We propose a method that allows for learning task-agnostic representations
based on value function estimates from a sequence of observations where the
last frame corresponds to a goal. These representations would learn to relate
states across different tasks, based on the temporal distance to the goal
state, irrespective of the appearance changes and dynamics. This method could
be used to transfer learnt policies/skills to unseen related tasks.

##  1 Introduction

Contrastive Learning has gained popularity in RL for unsupervised pretraining
DBLP:conf/icra/SermanetLCHJSLB18 (7). Recent works like
DBLP:conf/nips/EysenbachZLS22 (1) proposed an alternative view of RL in the
means of Contrastive Learning. Value Implicit Pretraining (VIP)
DBLP:conf/iclr/MaSJBK023 (3) proposes a method for learning an embedding space
that could be used for Zero-shot reward specification.

We discuss a method that learns consistent representations across a set of
different, but related tasks, using a sequence of observations. Although tasks
have differences in dynamics and appearance changes, they share a similar
objective based on reaching the goal frame. We compute the value estimates for
every state in a sequence using sum of discounted rewards to the goal frame.
States with similar value estimates are embedded closer to each other,
irrespective of the task they are sampled from.

Let
$\mathcal{T}_{train}=\\{\mathcal{T}_{1},\mathcal{T}_{2},...\mathcal{T}_{m}\\}$,
$\mathcal{T}_{test}=\\{\mathcal{T}_{1},\mathcal{T}_{2},...\mathcal{T}_{n}\\}$
be the set of train and test tasks, respectively. All the tasks in
$\mathcal{T}_{train}$ and $\mathcal{T}_{test}$ share a similar objective
function. During training, we assume access to the offline sequences of
observations for all the tasks
$\mathcal{D}=\\{\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3},...\mathcal{D}_{m+n}\\}$
for each of the tasks. Our model consists of an Encoder $f$ that learns to
encode images/observations into an embedding, which is taken as an input by
the policy $\pi$ to perform a task. Note that the encoder $f$ is trained only
using $\mathcal{D}$, that is a sequence of observations and rewards. The
objective of our method is to use the sequence of observations $\mathcal{D}$
to learn an encoder $f$ such that so that the policy could be zero-shot
transferred to unseen test tasks $\mathcal{T}_{test}$.

##  2 Overview of Methods for Pretraining

In this paper, we not only propose a novel idea that enables transfer across
different related games, but also describe the existing contrastive learning-
based pretraining methods. Following are the proposed methods in a
chronological order. All the following methods uses the problem description
mentioned in the previous section.

  1. 1.

Time Contrastive Learning (TCN): TCN involves sampling the positive within a
certain margin distance $d_{margin}$ from the anchor and a negative anywhere
from the positive to the end of the trajectory
DBLP:conf/icra/SermanetLCHJSLB18 (7). If the anchor is sampled at time instant
$t_{a}$, positive is sampled at $t_{p}$ and the negative sampled at $t_{n}$,
then $|t_{n}-t_{a}|>|t_{p}-t_{a}|$. We then use the standard triplet loss for
optimization.

  2. 2.

Positives from State Occupancy Measure (SOM): To learn good representations
for pretraining, we can embed 2 representations of states sampled within a
state occupancy measure DBLP:conf/nips/EysenbachZLS22 (1). State occupancy
measure from a specific instant is sampling within the following distribution:
DBLP:journals/corr/abs-2304-00046 (4):

$\rho(s|s_{t})=(1-\gamma)\sum_{\Delta t=1}^{H}\gamma^{\Delta t-1}P[S_{t+\Delta
t}=s|s_{t};\pi]$ (1)

  3. 3.

Value Implicit Pretraining (VIP): Unlike TCN, VIP DBLP:conf/iclr/MaSJBK023 (3)
learns temporally smooth embeddings by learning by embedding the goal
(positive) and the start (anchor) image close and the middle images
(negatives) as negatives within a sampled sub-trajectory. By learning this
objective, using multiple sampled sub-trajectories, the encoder learns
continuous embeddings across sequential frames in a trajectory. We then use
the InfoMAX loss DBLP:journals/corr/abs-1807-03748 (5) to optimize the encoder
parameters as shown below:

$\min_{\phi}\mathbb{E}_{x_{\text{pos}}}\left[-\log\frac{S_{\phi}\left(x,x_{\mathrm{pos}}\right)}{\mathbb{E}_{x_{\text{neg}}}S_{\phi}\left(x,x_{\text{neg}}\right)}\right]$
(2)

Figure 1: Using the pretrained encoder trained using our method that uses the
observation sequences from both the train games and the unseen games, we can
reuse the policy learnt on the train game, for an unseen game.

All the above methods perform pretraining as an optimization problem on a
large set of trajectories. Our method, not only learns an encoder that
performs the downstream task effectively but also embeds similar states of
related tasks close to each other. This would enable the learnt skills and
policies from one task to be reused in another.

##  3 Proposed Method

We learn a representation function (Encoder) that learns embeddings based on
the value estimates. The value estimates are the returns computed by the
Bellman equation, which is the discounted sum of rewards until the goal-frame.
We use a set of Shoot’em up Atari games for our experiments. For pretraining
the encoder, we use a sequence of observations obtained from the D4RL Datasets
DBLP:journals/corr/abs-2004-07219 (2). We convert these large sequences of raw
trajectories into smaller goal-based trajectories by considering an instant
where the agent gets a reward to be a goal. The goal states across all the
games are unified, as the objective of all the Shoot’em up games is to shoot
the enemy objects. We then compute value estimates of each frame in a sequence
using the Bellman equation. We then use contrastive learning to learn task-
agnostic, value-based representations that would allow transferring the policy
to an unseen related game. Following is the outline of our training procedure.

We first sample a game $x$ and then sample a trajectory within that game. We
randomly then sample two timesteps $i^{x}$ and $j^{x}$ such that
$(j^{x}-i^{x})<d_{thresh}$ and $\|v_{i^{x}}-v_{j^{x}}\|$, where $v_{i^{x}}$
and $v_{j^{x}}$ are the value estimates at the $i^{x}$ and $j^{x}$ frames.
$d_{thresh}$ and $v_{thresh}$ are the distance and the value thresholds,
respectively, that are fixed during pretraining. The value threshold
influences sampling by limiting $\|v_{i^{x}}-v_{j^{x}}\|$, since some games
have more denser reward signals. The distance threshold $d_{thresh}$ not only
allows for sampling the positive at $j^{x}$, but also limits where the
positive could be sampled from. This limit ensures that the negative, sampled
at $k_{x}$ is farther away from the anchor than the positive sample. We sample
the negative at timestep $k_{x}$, in the same trajectory, such that
$|i^{x}-j^{x}|<|i^{x}-k^{x}|$, where $i^{x}$, $j^{x}$ and $k^{x}$ are the
indices for anchor, positive and negative, respectively.

From the above process, we use $v_{i^{x}}$ and $|v_{j^{x}}-v_{i^{x}}|$ to
sample anchor, positive and negative observations for another game $y$. We
find a sample anchor at $i^{y}$ having value $v_{i^{y}}$ closest to
$v_{i^{x}}$. The positive sample for this specific anchor at $i^{y}$ is
obtained by sampling the closest element to value
$v_{i^{y}}+|v_{j^{x}}-v_{i^{x}}|$, within the distance threshold $d_{thresh}$.
We combine the anchors, positives and negatives for games $x$ and $y$ and use
Triplet loss function to embed the positives and the anchor close to each
other and the negatives far away from the anchor.

$\mathcal{L}(o_{an},o_{po},o_{ne})=\frac{1}{B}\sum_{i=0}^{i=B}\mathbf{max}\left[\mathbf{0},\;||f(o_{an})-f(o_{po})||_{2}^{2}-||f(o_{an})-f(o_{ne})||_{2}^{2}+\epsilon\right]$
(3)

In the above equation $f(o_{an})$, $f(o_{po})$ and $f(o_{ne})$ represent the
embeddings that are obtained after passing observation $o_{an}$, $o_{po}$ and
$o_{ne}$ (anchors, positives and negatives) and are of shape [$B$, $2$, $D$],
where $B$ is the batch size and $D$ is the embedding size. Note that the
middle index with value 2 is because of sampling a pair of games $(x,y)$ every
iteration, and is folded into the batch size during training.

Figure 2: We describe our pretraining strategy. We divide raw trajectories
into goal-based sub-trajectories, where the goal frame in a sequence of
observations corresponds to a common objective across all the games. Computed
value estimates for each frame in the sequence are shown above the frame. We
then use a contrastive learning-based pretraining method that learns task-
agnostic representations based on the proximity to the goal.

##  4 Results

We perform experiments by selecting 2 games DemonAttack and SpaceInvaders both
from Atari. Both of them belong to the Shoot’em up games and share the same
action space. The 5D action space consists of NoOp, Fire, Move Left, Move
Right, Move Left and Fire, and Move Right and Fire. Although the objective in
both the games is to shoot enemy objects, there are subtle difference across
games. For example, the agent has access to shields in SpaceInvaders that
allow the agent to stay under to prevent getting shot.

We first pretrain the encoder $f$ using the method described in the previous
section and visually shown in Figure 2. This is achieved by using a sequence
of unlabelled trajectories from both the games. Once we obtain the pretrained
encoder, we using an Online RL algorithm, in our case PPO
DBLP:journals/corr/SchulmanWDRK17 (6), to train a policy. The results using
different methods are shown in Figure 3.

We used the following hyperparameters for the 3 baselines and our proposed
method. For all the methods, we used a CNN architecture of 4 blocks of 16, 32,
256 and 512 channels with avg pooling and a train batch size of 32. For TCN we
sample the positive from a maximum margin of 4 timesteps. For VIP, we set the
minimum length of the subsequence and the maximum length of the sampled
subsequence to be 10 and 50 timesteps. Sequences could range from 5 timesteps
to 600 timesteps. For SOM, we set the $\gamma$ parameter to be $0.1$, which
means that the positive was sampled close to 1-2 timesteps from the anchor. As
for our method VEP (Value Explicit pretraining), we set the $d_{thresh}$ and
$v_{thresh}$ to be $0.2*$episode length and $0.01$ respectively.

##  5 Discussion

Transferring policies to novel but related tasks is an important problem that
needs to be addressed. We formulated a method to learn representations of
states from different tasks solely based on the temporal distance to the goal
frame. This way, the skills learnt from the train tasks could be transferred
to unseen related tasks. Apart from performing more comprehensive evaluations
for Atari, we are also experimenting with other similar settings, like that of
visual navigation.

Figure 3: Performance of different pretraining methods on the respective games
as mentioned above.

## References

  * (1) Benjamin Eysenbach, Tianjun Zhang, Sergey Levine and Ruslan Salakhutdinov  “Contrastive Learning as Goal-Conditioned Reinforcement Learning”  In _NeurIPS_ , 2022  URL: http://papers.nips.cc/paper\\_files/paper/2022/hash/e7663e974c4ee7a2b475a4775201ce1f-Abstract-Conference.html
  * (2) Justin Fu et al.  “D4RL: Datasets for Deep Data-Driven Reinforcement Learning”  In _CoRR_ abs/2004.07219, 2020  arXiv: https://arxiv.org/abs/2004.07219
  * (3) Yecheng Jason Ma et al.  “VIP: Towards Universal Visual Reward and Representation via Value-Implicit Pre-Training”  In _The Eleventh International Conference on Learning Representations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023_ OpenReview.net, 2023  URL: https://openreview.net/pdf?id=YJ7o2wetJ2
  * (4) Bogdan Mazoure et al.  “Accelerating exploration and representation learning with offline pre-training”  In _CoRR_ abs/2304.00046, 2023  DOI: 10.48550/arXiv.2304.00046
  * (5) Aäron Oord, Yazhe Li and Oriol Vinyals  “Representation Learning with Contrastive Predictive Coding”  In _CoRR_ abs/1807.03748, 2018  arXiv: http://arxiv.org/abs/1807.03748
  * (6) John Schulman et al.  “Proximal Policy Optimization Algorithms”  In _CoRR_ abs/1707.06347, 2017  arXiv: http://arxiv.org/abs/1707.06347
  * (7) Pierre Sermanet et al.  “Time-Contrastive Networks: Self-Supervised Learning from Video”  In _2018 IEEE International Conference on Robotics and Automation, ICRA 2018, Brisbane, Australia, May 21-25, 2018_ IEEE, 2018, pp. 1134–1141  DOI: 10.1109/ICRA.2018.8462891


# LibRPA: A Software Package for Low-scaling First-principles Calculations of
Random Phase Approximation Electron Correlation Energy Based on Numerical
Atomic Orbitals

Rong Shi  Min-Ye Zhang<EMAIL_ADDRESS>Peize Lin  Lixin He
<EMAIL_ADDRESS>Xinguo Ren<EMAIL_ADDRESS>
###### Abstract

LibRPA is a software package designed for efficient calculations of random
phase approximation (RPA) electron correlation energies from first principles
using numerical atomic orbital (NAOs). Leveraging a localized resolution of
identity (LRI) technique, LibRPA achieves $O(N^{2})$ or better scaling
behavior, making it suitable for large-scale calculation of periodic systems.
Implemented in C++ and Python with MPI/OpenMP parallelism, LibRPA integrates
seamlessly with NAO-based density functional theory (DFT) packages through
flexible file-based and API-based interfaces. In this work, we present the
theoretical framework, algorithm, software architecture, and installation and
usage guide of LibRPA. Performance benchmarks, including the parallel
efficiency with respect to the computational resources and the adsorption
energy calculations for H2O molecules on graphene, demonstrate its nearly
ideal scalability and numerical reliability. LibRPA offers a useful tool for
RPA-based calculations for large-scale extended systems.

## Program summary

  * 1.

Program title: LibRPA

  * 2.

Developer’s repository link: https://github.com/Srlive1201/LibRPA

  * 3.

Licensing provisions: LGPL

  * 4.

Programming language: C++, Python

  * 5.

Nature of problem: Calculating RPA electron correlation energies is
computationally expensive, typically scaling as $O(N^{4})$ with system size,
hindering its application to large-scale materials science problems.

  * 6.

Solution method: LibRPA utilizes the Localized Resolution of Identity (LRI)
technique, reducing computational scaling to $O(N^{2})$ or better. Implemented
in C++ and Python with MPI/OpenMP parallelization, it integrates with NAO-
based DFT packages, facilitating efficient and accurate RPA calculations for
large-scale periodic systems.

###### keywords:

First-principles calculations , Random-phase approximation , Numerical atomic
orbitals , Resolution of Identity , Low-scaling algorithm

††journal: Computer Physics Communications

[inst1] organization=CAS Key Laboratory of Quantum Information, University of
Science and Technology of China, city=Hefei, postcode=230026, state=Anhui,
country=China

[inst2] organization=Institute of Physics, Chinese Academy of Sciences,
city=Beijing, postcode=100190, country=China

[inst4] organization=The NOMAD Laboratory at the FHI Molecular Physics
Department of the Max-Planck-Gesellschaft, city=Berlin-Dahlem, postcode=14195,
state=, country=Germany

[inst3] organization=Songshan Lake Materials Laboratory, city=Dongguan,
postcode=523808, state=Guangdong, country=China

##  1 Introduction

Random phase approximation (RPA) [1, 2, 3] is a powerful method for
calculating the ground-state energy of interacting many-electron systems from
first principles [4, 5, 6, 7]. Formulated within the framework of adiabatic-
connection fluctuation-dissipation theorem (ACFDT) [8, 9], RPA can be viewed
as a fifth-rung functional according to the Jacob’s ladder [10] that
classifies different density functional approximations [11, 12]. Compared to
lower-rung functionals, RPA captures seamlessly nonlocal electron correlation
effects and provides a balanced description of different chemical bonding
interactions, rendering it particularly effective in describing delicate
energy differences. Benchmark calculations in the past have shown the
remarkable performance of RPA in describing properties such as surface
adsorption energies [13, 14], reaction barrier heights [15, 16], and the
energy hierarchy among different polymorphs of crystals [17, 18, 19, 20, 21]
or among different conformers of clusters [22, 23]. Furthermore, RPA is also a
key component of more sophisticated methodologies such as renormalized second-
order perturbation theory [16] or double hybrid functionals [24]. Despite its
attractive features, the application of RPA in materials science has been
limited by its high computational cost [25, 13, 7, 26, 27, 28, 29, 30], which
formally scales as $O(N^{4})$[31, 32, 33, 34, 35, 26] with respect to system
size $N$.

To reduce the computational cost of RPA for large systems, considerable recent
efforts have been devoting to developing low-scaling algorithms [36, 37, 38,
39, 40, 41, 42, 43, 44]. These include the space-time method that reduces the
scaling to $O(N^{3})$ [45, 46, 47]and the use of localized basis functions
combined with the resolution of identity (RI) approximation [34, 35, 26]. Most
recently, within the numerical atomic orbital (NAO) basis set framework,
paired with the localized resolution of identity (LRI) technique [48, 49, 50,
51, 52, 53, 54], we have developed an algorithm of formally $O(N^{2})$ scaling
that works both for finite systems and periodic systems with finite
${\mathbf{k}}$-point sampling [55]. Further exploiting the decay properties of
the Green’s functions in real space, one can achieve asymptotically linear
scaling for system sizes of approximately 1000 atoms.[55].

Building on these developments, here we introduce a software package called
LibRPA, designed for performing efficient and accurate first-principles RPA
correlation energy calculations using localized atomic orbitals (AOs). By
leveraging the LRI technique, LibRPA achieves $O(N^{2})$ or better scaling
behavior for calculating the key quantity in the RPA formalism, i.e., the non-
interacting Kohn-Sham (KS) density response function. The software is written
in C++ and Python, and it is massively parallelized by means of the MPI/OpenMP
hybrid scheme. LibRPA can be integrated with AO-based density functional
theory (DFT) packages (provided that the LRI infrastructure is available) via
flexible interfaces, making it a versatile tool for computational materials
research.

The paper is organized as follows. Section 2 presents the theoretical
framework behind the LibRPA implementation, which is followed by an exposition
of the software architecture and the implementation details in Sec. 3. Then,
in Sec. 4 we describe how the LibRPA package is installed and used in
practical calculations. The performance of LibRPA in parallel scalability is
reported in Sec. 5. A case study of the adsorption of a H2O molecule on
graphene using LibRPA is presented in Sec. 6. Finally, we summarize this work
in Sec. 7.

##  2 Theoretical framework

In this section, we briefly recapitulate the key concept and the main
equations behind the LibRPA implementation. A full exposition of these
equations can be found in Ref. [55].

### The Localized Resolution of Identity (LRI)

The localized resolution of identity (LRI) approximation, also known as pair
atomic density fitting (PADF) [48, 49] or concentric atomic density fitting
(CADF) [50, 51] in the literature, is the key numerical technique that enables
low-scaling algorithms in computational chemistry. It has been exploited
effectively in various contexts, including periodic hybrid functional
calculations [52, 53, 54, 56], $G_{0}W_{0}$ computations [57], as well as in
calculating the RPA correlation energies [58] and forces [23]. Benchmark
studies have demonstrated its sufficient accuracy for Hartree-Fock exchange
evaluations within the NAO framework for both molecular [59] and periodic
systems [52, 53, 54, 56]. While the LRI errors are more pronounced in
correlated methods involving unoccupied states, high-quality auxiliary basis
sets can mitigate these errors, making LRI applicable to correlated methods
like RPA and $GW$ for molecules [59, 23] and extended systems [60, 57].
Recently, a computational scheme by Spadetto et al. further refines PADF-RPA
for large molecules by selectively projecting out the problematic high-lying
states that cannot be well described by the local RI scheme [58].

In the context of using NAOs as the basis set, the LRI approach essentially
consists in approximating the product of two NAOs by a linear combination of
auxiliary basis functions (ABFs). This is mathematically expressed as

$\displaystyle\varphi_{i}({\mathbf{r}}-{\mathbf{R}}_{1}-{\bm{\tau}}_{I})\varphi_{k}({\mathbf{r}}-{\mathbf{R}}_{2}-{\bm{\tau}}_{K})$
$\displaystyle\approx$ $\displaystyle\sum_{\mu\in
I}C_{i({\mathbf{R}}_{1}),k({\mathbf{R}}_{2})}^{\mu({\mathbf{R}}_{1})}P_{\mu}({\mathbf{r}}-{\mathbf{R}}_{1}-{\bm{\tau}}_{I})+\sum_{\mu\in
K}C_{i({\mathbf{R}}_{1}),k({\mathbf{R}}_{2})}^{\mu({\mathbf{R}}_{2})}P_{\mu}({\mathbf{r}}-{\mathbf{R}}_{2}-{\bm{\tau}}_{K})$
$\displaystyle=$ $\displaystyle\sum_{\mu\in I}C_{i({\bf
0}),k({\mathbf{R}}_{2}-{\mathbf{R}}_{1})}^{\mu({\bf
0})}P_{\mu}({\mathbf{r}}-{\mathbf{R}}_{1}-{\bm{\tau}}_{I})+\sum_{\mu\in
K}C_{i({\mathbf{R}}_{1}-{\mathbf{R}}_{2}),k({\bf 0})}^{\mu({\bf
0})}P_{\mu}({\mathbf{r}}-{\mathbf{R}}_{2}-{\bm{\tau}}_{K})\,,$ (1)

where $\varphi_{i}({\mathbf{r}})$ and $\varphi_{k}({\mathbf{r}})$ denote the
NAOs centered on atoms $I$ and $K$, located at positions $\bm{\tau}_{I}$ and
$\bm{\tau}_{K}$ within the unit cell, respectively. The ABFs
$P_{\mu}({\mathbf{r}})$ are centered on either atom $I$ or $K$, and
$C_{i(\mathbf{R}_{1}),k(\mathbf{R}_{2})}^{\mu(\mathbf{R}_{1})}$ represent the
expansion coefficients of the product of NAOs in terms of the ABFs and the
lattice vector in parentheses indicates the unit cell to which the basis
function belongs. Such a notational system facilitates the usage of the LRI
formalism to periodic systems with finite ${\mathbf{k}}$-point sampling.

This approximation allows for a significant reduction in computational
complexity, especially when dealing with the evaluation of two-electron
integrals in large systems. By leveraging the locality and compactness of NAOs
and the associated sparsity of the expansion coefficients $C$, the LRI
technique enables efficient and scalable quantum-mechanical calculations.

### Green’s Function

The Green’s function is a central quantity in quantum-mechanical many-body
physics, describing the propagation behavior of quasiparticles and quasiholes
in an interacting system. For RPA calculations, the quantity of interest is
the non-interacting Green’s function $G^{0}$ of the KS system. Within the NAO
basis-set framework, the real-space imaginary-time expression for the non-
interacting Green’s function is given by

$G^{0}({\mathbf{r}},{\mathbf{r}^{\prime}},\mathrm{i}\tau)=\sum_{i,j}\sum_{{\mathbf{R}}_{1},{\mathbf{R}}_{2}}\varphi_{i}({\mathbf{r}}-{\mathbf{R}}_{1}-{\bm{\tau}}_{I})G^{0}_{i,j}({\mathbf{R}}_{2}-{\mathbf{R}}_{1},\mathrm{i}\tau)\varphi_{j}({\mathbf{r}^{\prime}}-{\mathbf{R}}_{2}-{\bm{\tau}}_{J})\,,$
(2)

where $G^{0}_{i,j}(\mathbf{R},i\tau)$ represents the matrix form of the
imaginary-time Green’s function expanded in terms of NAO basis,

$\displaystyle
G^{0}_{i,j}(\textbf{R},\mathrm{i}\tau)=\left\\{\begin{array}[]{ll}\displaystyle-\frac{\mathrm{i}}{N_{\mathbf{k}}}\sum_{n,{\mathbf{k}}}f_{n\textbf{k}}c_{i,n}(\textbf{k})c_{j,n}^{\ast}(\textbf{k})\mathrm{e}^{-\mathrm{i}\textbf{k}\cdot\textbf{R}}\mathrm{e}^{-(\epsilon_{n,\textbf{k}}-\mu)\tau}&\tau\leq
0\,,\\\
\displaystyle\frac{\mathrm{i}}{N_{\mathbf{k}}}\sum_{n,\textbf{k}}(1-f_{n\textbf{k}})c_{i,n}(\textbf{k})c_{j,n}^{\ast}(\textbf{k})\mathrm{e}^{-\mathrm{i}\textbf{k}\cdot\textbf{R}}\mathrm{e}^{-(\epsilon_{n,{\mathbf{k}}}-\mu)\tau}&\tau>0\,.\end{array}\right.$
(5)

Here, $c_{i,n}(\mathbf{k})$ are the coefficients of the KS eigenvectors,
$f_{n\mathbf{k}}$ the occupation numbers, $\epsilon_{n,\mathbf{k}}$ the KS
orbital energies, and $\mu$ the chemical potential. The real-space Green’s
function will be used to construct the non-interacting response function
matrix, as described below.

### Non-Interacting Response Function ($\chi^{0}$)

The response function $\chi^{0}$ is an essential quantity for computing the
RPA correlation energy and various electronic response properties. In real-
space and imaginary-time domain, $\chi^{0}$ can be expressed in terms of the
Green’s function as [45, 46, 47]

$\chi^{0}(\mathbf{r},\mathbf{r^{\prime}},\mathrm{i}\tau)=-\mathrm{i}G^{0}(\mathbf{r},\mathbf{r^{\prime}},i\tau)G^{0}(\mathbf{r^{\prime}},\mathbf{r},-\mathrm{i}\tau)\,.$
(6)

By employing the expression for the Green‘s function in the basis of NAOs as
given in Eq. (2), $\chi^{0}$ can be expanded in terms of NAOs ,

$\displaystyle\chi^{0}(\textbf{r},\textbf{r}^{\prime},\mathrm{i}\tau)=-\mathrm{i}\sum_{i,j,k,l}\sum_{{\mathbf{R}}_{1},{\mathbf{R}}_{2},{\mathbf{R}}_{3},{\mathbf{R}}_{4}}$
$\displaystyle\varphi_{i}({\mathbf{r}}-{\mathbf{R}}_{1}-{\bm{\tau}}_{I})\varphi_{k}({\mathbf{r}}-{\mathbf{R}}_{3}-{\bm{\tau}}_{K})G^{0}_{i,j}({\mathbf{R}}_{2}-{\mathbf{R}}_{1},\mathrm{i}\tau)$
$\displaystyle
G^{0}_{l,k}({\mathbf{R}}_{3}-{\mathbf{R}}_{4},-\mathrm{i}\tau)\varphi_{j}({\mathbf{r}^{\prime}}-{\mathbf{R}}_{2}-{\bm{\tau}}_{J})\varphi_{l}({\mathbf{r}^{\prime}}-{\mathbf{R}}_{4}-{\bm{\tau}}_{L})\,.$
(7)

Further applying the LRI approximation leads to a more compact expression of
$\chi^{0}$ in terms of ABFs,

$\chi^{0}({\mathbf{r}},{\mathbf{r}^{\prime}},\mathrm{i}\tau)=\sum_{\mu\in\mathcal{U},\nu\in\mathcal{V}}\sum_{\textbf{R}_{1},\textbf{R}_{2}}P_{\mu}({\mathbf{r}}-\textbf{R}_{1}-{\bm{\tau}}_{\mathcal{U}})\chi^{0}_{\mu,\nu}(\textbf{R}_{2}-\textbf{R}_{1},\mathrm{i}\tau)P_{\nu}({\mathbf{r}^{\prime}}-\textbf{R}_{2}-{\bm{\tau}}_{\mathcal{V}})\,.$
(8)

By combining Eq. (7), Eq. (8) and Eq. (1), one finally arrives at the matrix
form of $\chi^{0}$ in the real-space imaginary-time domain,

$\displaystyle\chi_{\mu,\nu}^{0}(\textbf{R},\mathrm{i}\tau)$
$\displaystyle=-\mathrm{i}\left[\sum_{i\in\mathcal{U}}\sum_{k,\textbf{R}_{1}}C_{i(\textbf{0}),k(\textbf{R}_{1})}^{\mu(\textbf{0})}\left(\sum_{j\in\mathcal{V}}G_{i,j}(\textbf{R},\mathrm{i}\tau)\sum_{l,\textbf{R}_{2}}C_{j(\textbf{0}),l(\textbf{R}_{2}-\textbf{R})}^{\nu(\textbf{0})}G_{l,k}(\textbf{R}_{1}-\textbf{R}_{2},-\mathrm{i}\tau)\right.\right.$
$\displaystyle+\sum_{j\in\mathcal{V}}G^{*}_{i,j}(\textbf{R},-\mathrm{i}\tau)\sum_{l,\textbf{R}_{2}}C_{j(\textbf{0}),l(\textbf{R}_{2}-\textbf{R})}^{\nu(\textbf{0})}G^{*}_{l,k}(\textbf{R}_{1}-\textbf{R}_{2},\mathrm{i}\tau)$
$\displaystyle+\sum_{j\in\mathcal{V}}G_{j,k}(\textbf{R}_{1}-\textbf{R},-\mathrm{i}\tau)\sum_{l,\textbf{R}_{2}}C_{j(\textbf{0}),l(\textbf{R}_{2}-\textbf{R})}^{\nu(\textbf{0})}G_{i,l}(\textbf{R}_{2},\mathrm{i}\tau)$
$\displaystyle\left.\left.+\sum_{j\in\mathcal{V}}G^{*}_{j,k}(\textbf{R}_{1}-\textbf{R},\mathrm{i}\tau)\sum_{l,\textbf{R}_{2}}C_{j(\textbf{0}),l(\textbf{R}_{2}-\textbf{R})}^{\nu(\textbf{0})}G^{*}_{i,l}(\textbf{R}_{2},-\mathrm{i}\tau)\right)\right]\,.$
(9)

The derivation of the above equation takes advantage of the symmetry
attributes of the expansion coefficients along with index permutation
techniques. For an in-depth explanation of these derivations, please refer to
Ref. [55].

The evaluation of $\chi^{0}_{\mu,\nu}(\textbf{R},\mathrm{i}\tau)$ represents
the computational bottleneck in the conventional algorithm. The essential
point of our work is to reduce the computational scaling of evaluating
$\chi^{0}_{\mu,\nu}(\textbf{R},\mathrm{i}\tau)$, which is crucial for the
overall efficiency of LibRPA. To this end, we have designed an algorithm based
on atomic pair assignments which is illustrated in Algorithm 1.

1 forall  _ $<\textbf{R},\tau>$_ do

2 forall  _ $<\mathcal{U}(\textbf{0}),\mathcal{V}(\textbf{R})>$_ do

3 forall  _ $L(\textbf{R}_{2})\in\mathcal{N}[\mathcal{V}(\textbf{R})]$ _ do

4 Calculate $X^{\nu}_{i,j}(\textbf{R},\mathrm{i}\tau)=\displaystyle\sum_{l\in
L,\textbf{R}_{2}}C_{j(\textbf{0}),l(\textbf{R}_{2}-\textbf{R})}^{\nu(\textbf{0})}G_{i,l}(\textbf{R}_{2},\mathrm{i}\tau)$;

5

6 forall  _ $K(\textbf{R}_{1})\in\mathcal{N}[\mathcal{U}(\textbf{0})]$ _ do

7 forall  _ $L(\textbf{R}_{2})\in\mathcal{N}[\mathcal{V}(\textbf{R})]$ _ do

8 Calculate
$N^{\nu}_{j,k}(\textbf{R}_{1},\textbf{R},\mathrm{i}\tau)=\displaystyle\sum_{l\in
L,\textbf{R}_{2}}C_{j(\textbf{0}),l(\textbf{R}_{2}-\textbf{R})}^{\nu(\textbf{0})}G_{l,k}(\textbf{R}_{1}-\textbf{R}_{2},-\mathrm{i}\tau)$;

9

10 Calculate
$M^{\nu}_{i,k}(\textbf{R}_{1},\textbf{R},\mathrm{i}\tau)=\displaystyle\sum_{j\in\mathcal{V}}G_{i,j}(\textbf{R},\mathrm{i}\tau)N^{\nu}_{j,k}(\textbf{R}_{1},\textbf{R},\mathrm{i}\tau)$;

11 Calculate
$Z^{\nu}_{i,k}(\textbf{R}_{1},\textbf{R},\mathrm{i}\tau)=\displaystyle\sum_{j\in\mathcal{V}}G_{j,k}(\textbf{R}_{1}-\textbf{R},-\mathrm{i}\tau)X^{\nu}_{i,j}(\textbf{R},\mathrm{i}\tau)$;

12 Calculate
$O^{\nu}_{i,k}(\textbf{R}_{1},\textbf{R},\mathrm{i}\tau)=M^{\nu}_{i,k}(\textbf{R}_{1},\textbf{R},\mathrm{i}\tau)+M^{\nu*}_{i,k}(\textbf{R}_{1},\textbf{R},-\mathrm{i}\tau)+Z^{\nu}_{i,k}(\textbf{R}_{1},\textbf{R},\mathrm{i}\tau)+Z^{\nu*}_{i,k}(\textbf{R}_{1},\textbf{R},-\mathrm{i}\tau)$;

13 Calculate
$\chi^{0}_{\mu,\nu}(\textbf{R},\mathrm{i}\tau)\mathrel{+}=C_{i(\textbf{0}),k(\textbf{R}_{1})}^{\mu(\textbf{0})}O^{\nu}_{i,k}(\textbf{R}_{1},\textbf{R},\mathrm{i}\tau)$;

14

15

16

Algorithm 1 Loop structure of evaluating
$\chi^{0}_{\mu\nu}({\mathbf{R}},\mathrm{i}\tau)$.
$\langle\mathcal{U}(\textbf{0}),\mathcal{V}(\textbf{R})\rangle$ denotes an
atomic pair with the atom ${\mathcal{U}}$ in the original unit cell ${\bf 0}$
and the atom ${\mathcal{V}}$ in the unit cell ${\mathbf{R}}$. The symbol
$\mathcal{N}[\mathcal{U}]$ represents the set of neighboring atoms around the
atom $\mathcal{U}$, and
$K(\textbf{R}_{1})\in\mathcal{N}[\mathcal{U}(\textbf{0})]$ means that the atom
$K$ (in the unit cell ${\mathbf{R}}_{1}$) belongs to the neighborhood of the
atom $\mathcal{U}$ (in the original unit cell).

The core of our algorithm involves breaking down the evaluation of
$\chi^{0}_{\mu,\nu}(\textbf{R},\mathrm{i}\tau)$ into manageable blocks
associated with individual atomic pairs. The evaluation process starts by
iterating over time grid points ${\tau_{j}}$ and lattice vectors R within the
Born-von Karman (BvK) supercell. For each $(\tau,\textbf{R})$ pair, the
response function matrix is partitioned into blocks corresponding to
individual atomic pairs
$\langle\mathcal{U}(\textbf{0}),\mathcal{V}(\textbf{R})\rangle$. By computing
these blocks separately for each atomic pair and assembling them, we retrieve
the full $\chi^{0}_{\mu,\nu}$ matrix. for a given R, the number of such atomic
pairs is $N_{\rm at}^{2}$, with $N_{\rm at}$ being the number of atoms in a
unit cell.

As illustrated in Algorithm 1, within the loop over atomic pairs, we run
through the atoms $K$ in the unit cell of $\textbf{R}_{1}$ and $L$ in the unit
cell of $\textbf{R}_{2}$ to compute the intermediate quantities like
$N^{\nu}_{j,k}(\textbf{R}_{1},\textbf{R},\mathrm{i}\tau)$ and
$X^{\nu}_{i,j}(\textbf{R},\mathrm{i}\tau)$ (defined in Algorithm 1, and
further details can be found in Ref. [55]). The crucial point is that the
atoms $K$ and $L$ must be neighbors of $\mathcal{U}(\textbf{0})$ and
$\mathcal{V}(\textbf{R})$, respectively. Otherwise, the LRI expansion
coefficients are negligible, so that the contribution from distant atoms
become insignificant. The number of neighboring atoms is determined by the
cutoff radii of NAO basis functions and does not increase with system size,
leading to a computational cost that approaches a constant for each block of
the response function matrix associated with an atomic pair. Consequently, the
total computational cost scales as $O(N_{\rm at}^{2}N_{\mathbf{R}}N_{\tau})$,
which is effectively $O(N^{2})$ for evaluating the response function matrix.
If one further takes into account of the spatial decay behavior of the Green’s
function, the range of the lattice vector ${\mathbf{R}}$ will also be limited,
and then the algorithm becomes asymptotically $O(N)$ [55].

### RPA Correlation Energy ($E_{c}^{\text{RPA}}$)

So far, we have developed the algorithm to compute the real-space imaginary-
time response-function matrix represented in terms of ABFs. However, for the
calculation of the RPA correlation energy, it is more convenient to transform
the response function $\chi^{0}$ to the $\mathbf{k}$ space and the imaginary-
frequency domain. Then, the computation of the RPA correlation energy is
executed through an integral across the imaginary frequency spectrum, which
incorporates the non-interacting response function $\chi^{0}$ and the Coulomb
potential $v$,

$E^{\text{RPA}}_{c}=\frac{1}{2\pi}\int^{\infty}_{0}\mathrm{d}\omega\,\text{Tr}\left[\ln(1-\chi^{0}(\mathrm{i}\omega)v)+\chi^{0}(\mathrm{i}\omega)v\right]\,,$
(10)

where the trace is taken over the matrix form of $\chi^{0}(\mathrm{i}\omega)$
and $v$, which are usually represented in terms of the ABF basis set.

The transformation naturally involves two steps. In the first step, $\chi^{0}$
is converted from the imaginary-time domain to the imaginary-frequency domain,
which can be achieved through a cosine Fourier transform due to the symmetry
properties of $\chi^{0}$ [61, 62, 47],

$\chi^{0}_{\mu,\nu}(\mathbf{R},\mathrm{i}\omega_{k})=-\mathrm{i}\sum_{j=1}^{N_{\tau}}\gamma_{jk}\chi^{0}_{\mu,\nu}(\mathbf{R},\mathrm{i}\tau_{j})\,\cos(\tau_{j}\omega_{k})\,.$
(11)

Here, $\gamma_{jk}$ are coefficients determined through an $L^{2}$
minimization procedure, and the use of nonuniform grids for $\tau_{j}$ and
$\omega_{k}$ from the GreenX library[63, 64, 65] ensures an efficient and
accurate representation of the frequency-dependent response. The second step
consists in the Fourier transform from real space to reciprocal space, which
is straightforward due to the periodicity of the system:

$\chi^{0}_{\mu,\nu}(\mathbf{q},\mathrm{i}\omega)=\sum_{\mathbf{R}}\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{R}}\chi^{0}_{\mu,\nu}(\mathbf{R},\mathrm{i}\omega)\,.$
(12)

Now, with the $\chi^{0}$ matrix ready in the reciprocal space and imaginary-
frequency domain, we can proceed to compute the RPA correlation energy. The
Coulomb matrix in the basis of ABFs, crucial for this computation, is given
by:

$V_{\mu\nu}({\mathbf{q}})=\sum_{{\mathbf{R}}}\mathrm{e}^{\mathrm{i}{\mathbf{q}}\cdot{\mathbf{R}}}V_{\mu\nu}({\mathbf{R}})=\sum_{{\mathbf{q}}}\mathrm{e}^{\mathrm{i}{\mathbf{q}}\cdot{\mathbf{R}}}\iint\frac{P_{\mu}({\mathbf{r}}-{\bm{\tau}}_{\cal
U})P_{\nu}({\mathbf{r}^{\prime}}-{\bm{\tau}}_{\cal
V}-{\mathbf{R}})}{|{\mathbf{r}}-{\mathbf{r}^{\prime}}|}\textrm{d}{\mathbf{r}}\textrm{d}{\mathbf{r}^{\prime}}\,,$
(13)

which represents the Fourier transform of the real-space Coulomb matrix
elements between ABFs. Combining the transformed $\chi^{0}$ and the Coulomb
matrix, we introduce an intermediate quantity,
$\Pi(\mathbf{q},\mathrm{i}\omega)=\chi^{0}(\mathbf{q},\mathrm{i}\omega)V(\mathbf{q})$,
to facilitate the computation of $E_{c}^{RPA}$,

$\displaystyle E^{\text{RPA}}_{\text{c}}$
$\displaystyle=\frac{1}{2\pi}\frac{1}{N_{\mathbf{q}}}\sum_{\mathbf{q}}\int^{\infty}_{0}\mathrm{d}\omega\,\text{Tr}\left[\ln(1-\Pi(\mathbf{q},\mathrm{i}\omega))+\Pi(\mathbf{q},\mathrm{i}\omega)\right]$
$\displaystyle=\frac{1}{2\pi}\frac{1}{N_{\mathbf{q}}}\sum_{\mathbf{q}}\int^{\infty}_{0}\mathrm{d}\omega\,\ln\left[\det(1-\Pi(\mathbf{q},\mathrm{i}\omega))\right]+\text{Tr}[\Pi(\mathbf{q},\mathrm{i}\omega)]\,,$
(14)

where the property $\textrm{Tr}[\textrm{ln}(A)]=\textrm{ln}[\textrm{det}(A)]$
is used. The integration over the imaginary frequency domain, coupled with the
summation over the ${\mathbf{k}}$ points in the Brillouin zone, yields the RPA
correlation energy per unit cell.

##  3 Software Architecture

###  3.1 Architecture

LibRPA is an independent computer program designed for the implementation of
the real-space low-scaling RPA algorithm within the NAO basis set framework.
Fig. 1 presents the overall architecture of LibRPA, illustrating its
interaction with external DFT software and its efficient internal
computational process.

Figure 1: The architecture of LibRPA. The LibRPA package contains: 1) An
interface which allows to integrate the LibRPA with external DFT software via
file-based or API methods; 2) core computational layer that performs the
actual RPA calculations, including response function and correlation energy;
3) the utilities layer that provides matrix manipulation, parallel computing,
data I/O, and format conversion tools; and 4) dependencies layer that links
external libraries like MPI/OpenMP, ScaLapack, LibRI, and GreenX for optimized
performance.

At the interface level, LibRPA can be used in two ways, namely, file-based
integration through standalone program and integration through application
programming interface (API), depending on how it interacts with external DFT
software. In the former case, LibRPA runs as a post-processing program
(driver) that reads data files produced by a DFT software and performs
subsequent RPA calculations. In the latter case, LibRPA functions, implemented
in a library file (such as librpa.so on Linux system), are called by a DFT
software for runtime data transfer and processing.

The core computational layer of LibRPA encapsulates a series of highly
optimized algorithmic functions that are pillars of RPA computations. They are
responsible for performing various complex mathematical operations, including
the evaluation of the response function and the computation of the RPA
correlation energy, ensuring both sufficient accuracy and high efficiency.

The utilities layer provides essential matrix manipulation and parallel
computing services, supporting LibRPA’s high-performance operation. These
tools are developed for efficient processing of large datasets and enhancing
the management efficiency of multithreading computing tasks, which are crucial
for ensuring the high overall computational speed and efficiency of LibRPA.

Finally, the layer of dependencies contains a suite of external libraries
required for LibRPA’s operation, such as MPI/OpenMP and Scalapack, providing
necessary parallel computing support and optimized numerical operation.
Additionally, the GreenX library is employed for generating imaginary
frequency/time integration grids [64, 65], and the LibRI library (optional for
RPA correlation energy) can be used for efficient tensor contractions [66].
The integration of these libraries ensures that LibRPA can run stably on
various hardware and operating system platforms, while efficiently utilizing
the computational power of modern multicore processors.

In summary, LibRPA’s architecture design balances high efficiency and
flexibility, offering a top-tier RPA computational solution, whether running
as an independent software or closely integrated with DFT softwares.

###  3.2 Massive parallelization

LibRPA is massively parallelized, leveraging the power of MPI/OpenMP hybrid
parallelism to accelerate the computational tasks. Throughout the LibRPA
workflow, the calculation of the non-interacting response function matrix
$\chi^{0}$ represents the most time-consuming and critical component. To
streamline this process, we have implemented two task distribution schemes,
which are chosen depending on the system size and the density of
${\mathbf{k}}$ points.

To explain the parallelization schemes, we define two sets of indices
associated with the $\chi^{0}$: ${\cal P}$ represents all possible atomic
pairs $(\mathcal{U,V})$, where the number of such pairs is given by
$N_{\text{atom}}(N_{\text{atom}}+1)/2$ for a system with $N_{\text{atom}}$
atoms in the unit cell. ${\cal Q}$ represents all possible combinations of the
lattice vectors ($N_{{\mathbf{R}}}$) in the BvK supercell and time grid points
($N_{\tau}$). The complete response function can be represented as a union of
all computed $\chi^{0}_{\mathcal{U,V}}({\mathbf{R}},\tau)$ values across the
atomic pairs $(\mathcal{U,V})$ in ${\cal P}$ and the $({\mathbf{R}},\tau)$
grid points in ${\cal Q}$:

$\chi^{0}_{\text{total}}=\bigcup_{(\mathcal{U,V})\in{\cal
P}}\;\bigcup_{({\mathbf{R}},\tau)\in{\cal
Q}}\chi^{0}_{\mathcal{U,V}}({\mathbf{R}},\tau)\;.$ (15)

The first scheme, being atom-pair leading, distributes subtasks based on
atomic pairs. It partitions the global $\chi^{0}$ matrix into blocks
associated with these pairs, making it particularly advantageous for large-
scale systems. In this context, each parallel subtask is responsible for
building the block matrices for all $({\mathbf{R}},\tau)$ points associated
with a given atomic pair, thus optimizing the distribution of computational
loads across multiple processors. Notably, the atom-pair-leading disribution
scheme does not require communication until the calculation of the $\Pi$
matrix is complete.

Alternatively, the second scheme, being $({\mathbf{R}},\tau)$-leading, is
optimized for smaller systems or those with a densely sampled Brillouin zone.
It allocates tasks based on $({\mathbf{R}},\tau)$ points, whereby each
parallel subtask calculates the whole $\chi^{0}_{\mu,\nu}$ matrix on specific
grid points. In this scheme, communication is required during the Fourier
transform. However, since the $\Pi$ matrix for each $({\mathbf{R}},\tau)$
point is complete (not distributed) on a single process, the LU decomposition
can be performed independently without further communication.

These two parallelization schemes, different in their data distributions and
the requirement of interprocess communications, are designed to address the
computational demands in different scientific scenarios. In practice, they are
automatically chosen by LibRPA based on specific system characteristics, and
by so doing, one can maximize the use of parallel computing resources, ensure
efficient load balancing, and minimize interprocess communication overhead.
The flow diagrams in A provide a more detailed explanation of the distinct
subtask distributions and the specific points at which inter-process
communication occurs within the two parallelization schemes.

LibRPA further enhances computational efficiency through OpenMP, which refines
the subtasks assigned to MPI processes by dynamically distributing them among
threads. This finer subtask granularity utilizes the shared memory
capabilities of OpenMP, reducing overall memory consumption while maintaining
computational throughput. By optimizing memory usage across threads, OpenMP
enables LibRPA to handle large-scale computations more efficiently, ensuring
that resources are utilized effectively within the multi-threaded environment.

###  3.3 Interface

LibRPA’s interface design is geared towards consistency and compatibility
across AO-based DFT software packages written in different programming
languages. This applies both to the driver and to the API interfaces.

For the driver interface, we use standardized data formats for input files of
essential physical quantities, including the KS eigenvectors and the expansion
coefficients $C_{i(\mathbf{0}),j(\mathbf{R})}^{\mu(\mathbf{0})}$, as well as
the Coulomb matrix $V_{\mu,\nu}$ in the ABF representation. Export of
essential data in standardized format has been implemented in both FHI-aims
[67] and ABACUS [68]. This standardization facilitates consistent data
exchange and integration across different DFT software packages.

For integration through API, the API functions are exposed with C linkage and
use one-dimensional array pointers for transferring matrices and tensors, as
illustrated in B. This makes the binding to other programming languages easy.
Furthermore, multiple API functions for the transfer of Coulomb matrix are
provided to address varied data distribution and parallel schemes of the
hosting programs. The target quantity, e.g. RPA correlation energy, can be
accessed using dedicated function and further processed within the execution
of the host program.

##  4 Installation and Usage of LibRPA

As mentioned above, LibRPA can work in two modes: either as a standalone
program or as a dynamic library (librpa.so) to be integrated into DFT software
packages that utilize the NAO basis sets. To facilitate its accessibility, we
have streamlined the installation and usage process of LibRPA. In the
following, we outline the necessary steps to install LibRPA and explain how to
incorporate it into a computational workflow. For detailed usage instructions,
including concrete steps to configure and execute RPA calculations, we refer
the readers to the documentation page of the LibRPA GitHub repository [69].

###  4.1 Installation

To embark on the installation, one should first clone the LibRPA repository
from GitHub:

    git clone https://github.com/Srlive1201/LibRPA.git

With the LibRPA package available in a local directory, it is imperative to
configure the cmake.inc file to align with one’s specific computational
environment. This configuration is crucial for setting the compiler flags and
specifying the paths to dependencies, thereby ensuring compatibility and
achieving optimal performance.

The installation proceeds with the CMake-based build system. To compile LibRPA
and generate the librpa.so dynamic library, execute the following commands:

    cmake -B build
    cmake --build build

###  4.2 Usage

As mentioned above, LibRPA offers dual operational modes: as an independent
executable for direct RPA calculations and as a dynamic library for
integration with DFT software such as FHI-aims and ABACUS. This versatility
enhances its utility in different computational scenarios.

####  4.2.1 Standalone Operation

In its standalone mode, LibRPA employs a driver program to invoke
functionalities from librpa.so, allowing it to perform RPA calculations based
on external input data files, generated beforehand by NAO-based DFT program.
This mode separates the execution of the DFT program and LibRPA, which is
particularly convenient for debugging or developing new features within
LibRPA.

####  4.2.2 Integration within DFT Software

When integrating LibRPA as a dynamic library and using LibRPA through API
calls, it is necessary to configure the DFT software to recognize and utilize
librpa.so. This requires adjusting the DFT software’s build process and/or
setting environment variables to include the path to the LibRPA library.

##  5 Performance

Figure 2: Strong scaling of the wall time as a function of the number of CPU
cores for calculating the RPA correlation energy using LibRPA. The test
systems are the C diamond supercells of different sizes, containing 128 (red
circles), 250 (blue diamonds) and 432 (green squares) atoms. The calculations
were performed using Intel Xeon-Silver-4208 CPUs, each with 8 cores, a base
frequency of 2.10 GHz, and a maximum frequency of 3.20 GHz.

Figure 2 presents the strong-scaling behavior of LibRPA with respect to the
computational resources. What is plotted are the wall-clock timings for RPA
correlation energy calculations for diamond supercells containing 128, 250,
and 432 carbon atoms, respectively. The test calculations were performed by
interfacing LibRPA with ABACUS using the DZP basis set and a single
$\Gamma$-point in the Brillouin-zone sampling. In the LibRPA calculations, 12
minimax time/frequency grid points are used. An ideal scaling of
parallelization follows a straight line when plotting the computation time
against the reciprocal of the number of CPU cores. Figure 2 shows the nearly
ideal scaling behavior of LibRPA with respect to computational resources,
indicating efficient expansion with increased numbers of CPU cores. For the
specific test examples in Fig. 2, the calculation took 237 minutes for the
largest system (432 atoms), using 192 CPUs (2 nodes, each with 96 CPU cores).
Doubling the computing resources to 384 CPUs (4 nodes) reduced the computation
time to 121 minutes; and further doubling the resources to 768 CPUs (8 nodes)
decreased the time to 62 minutes. A polynomial fitting of the curve of the
computation time (in minutes) for the 432-atom test case is
$t=44736/{N_{\text{cores}}}+4.383(R^{2}=0.9999)$. This reflects the excellent
balance in task distribution achieved in the design and implementation of the
parallel algorithms.

##  6 Case Study: H2O Adsorption Energy on Graphene

As a simple case study, we use LibRPA to compute the adsorption energy of a
water molecule on graphene. Understanding the adsorption behavior of H2O
molecules on graphene surfaces is an important theme in studying graphene’s
surface properties, including the mechanisms behind its hydrophilicity and
hydrophobicity [70]. Additionally, the adsorption and diffusion behavior of
H2O molecules on graphene surfaces governs the stability and reactivity of
graphene in various environments [71, 72]. From an application perspective,
graphene’s unique physicochemical properties make it highly promising in
fields such as energy storage, catalysis, and separation technologies. For
instance, the adsorption characteristics of H2O molecules on graphene directly
impact its electrochemical performance when used as an electrode material in
supercapacitors and batteries [73, 74].

A H2O molecule can adsorb on graphene surfaces in different configurations,
among which the one-leg and two-leg configurations are of particular interest.
In the one-leg configuration, one O-H bond of the water molecule points
perpendicularly to the graphene plane, while in the two-leg configuration, the
C2 axis of rotation of the water molecule is perpendicular to the plane, with
both hydrogen atoms oriented towards the graphene surface. Lower-rung DFT
studies often suggested that the one-leg configuration of water adsorption on
graphene is more stable [75, 71]. In contrast, previous RPA calculations have
shown that the water molecule binds more strongly with the graphene surface
under the two-leg configuration [75, 76], which has been verified by
experiment [76]. These findings underscore the critical role of RPA in
providing accurate and reliable results for water-graphene interactions.

Here we use LibRPA, interfaced with FHI-aims, to calculate the RPA adsorption
energies of H2O on graphene in the one-leg and two-leg configurations. The
adsorption energy is defined as

$\displaystyle
E_{\text{ads}}=E_{\ce{H2O}/\text{graphene}}-E_{\ce{H2O}}-E_{\text{graphene}}\,,$
(16)

where $E_{\ce{H2O}/\text{graphene}}$, $E_{\text{graphene}}$, and
$E_{\ce{H2O}}$ are the total energies of the H2O-graphene system, clean
graphene, and isolated H2O molecule, respectively. In our calculations, a
$4\times 4$ graphene supercell was used, with a $4\times 4\times 1$
$\mathbf{k}$-point grid. FHI-aims provides the PBE reference states and
corresponding LRI expansion coefficients and Coulomb matrix, while the RPA
correlation energy is calculated by LibRPA. For imaginary frequency and time
grids, 16 minimax grid points are used, and for Green’s function-based
screening, a threshold of $10^{-3}$ is used. These settings provide adequate
accuracy.

(a) Two leg

(b) One leg

Figure 3: Adsorption energies of a H2O molecule adsorbed on graphene as a
function of the height of the O atom above the graphene, obtained using LibRPA
interfaced with FHI-aims. The left and right panels show the results for the
adsorbed water molecule in its two-leg and one-leg configurations,
respectively. The three curves correspond to the results obtained using the
loc-NAO-VCC-3Z and loc-NAO-VCC-4Z basis sets, and the extrapolated results at
the CBS limit using the CBS(3,4) method.

Figure 3 presents the adsorption energy of H2O molecule on graphene as a
function of the distance between the O atom and the graphene surface under the
two-leg (left) and one-leg (right) configurations. The results obtained using
the loc-NAO-VCC-3Z and loc-NAO-VCC-4Z [60] basis sets, as well as the
extrapolated complete basis set (CBS) results from these two basis sets
(denoted CBS(3,4)), are shown. The basis convergence behavior shows that the
adsorption energy increases in magnitude (indicating stronger adsorption) as
the basis size grows. These results highlight the importance of using high-
quality basis sets and the extrapolation method to attain reliable RPA
energies. Furthermore, we note that the basis set superposition errors (BSSEs)
are accounted for in our calculations of the RPA adsorption energies. The way
in which the BSSEs are corrected in the present calculation is explained in C.
Comparing the results in the left and right panels, one finds that the water
molecule binds stronger to graphene in the two-leg configuration than in the
one-leg configuration, consistent with previous studies [75, 76].

Table 1: Adsorption energy $E_{\text{ads}}$ and height of H2O in the two-leg and one-leg configurations on graphene calculated by RPA and other methods. The RPA calculations are based on the PBE reference. The plane-wave results are taken from from Ref. [75], whereby the PBE and RPA results were obtained using the VASP code, and the diffusion Monte Carlo (DMC) results are obtained using the CASINO code.  |  | Two leg | One leg
---|---|---|---
Approach | Basis set |  $E_{\text{ads}}$ (meV) | Height (Å) |  $E_{\text{ads}}$ (meV) | Height (Å)
DMCa | Plane waves | $-70\pm 10$ | 3.4-4.0 | $-70\pm 10$ | 3.4-4.0
RPAa | Plane waves | $-98$ | 3.42 | $-82$ | 3.55
PBEa | Plane waves | $-27$ | 3.65 | $-31$ | 3.65
This work |  |  |
RPA | loc-NAO-VCC-3Z | $-72$ | 3.50 | $-67$ | 3.55
| loc-NAO-VCC-4Z | $-97$ | 3.40 | $-84$ | 3.50
| CBS(3,4) | $-116$ | 3.40 | $-98$ | 3.50
PBE | loc-NAO-VCC-4Z | $-31$ | 3.68 | $-34$ | 3.66
aReference [75] |  |  |  |  |

Table 1 presents the actual adsorption energies and equilibrium heights of H2O
molecules on graphene in the two configurations and compares them with the
results obtained in a previous study using VASP [75]. The table shows that our
RPA results using LibRPA and the loc-NAO-VCC-$n$Z basis sets are in overall
satisfactory agreement with those reported in the literature [75] with respect
to adsorption energies and equilibrium heights. Specifically, the results
obtained using loc-NAO-VCC-4Z basis set are in quantitative agreement with the
previous plane-wave results, while the CBS(3,4) extrapolated results suggest a
stronger binding by an amount of 20 meV in magnitude. This discrepancy might
be due to the pseudopotentials used in the plane-wave calculations, or to the
remaining uncertainties in the basis set extrapolation scheme, or to the
particular procedure for correcting the BSSEs employed in the present work
(cf. C). Nevertheless, the present case study demonstrates the capability of
LibRPA to be applied to surface adsorption problems.

##  7 Summary

In this work, we presented the underlying theoretical foundation and
algorithmic details behind the LibRPA software package, which aims to
calculate the KS response function and the RPA correlation energy within the
AO basis set framework. In addition, we described the architecture of the
LibRPA, the design of its interface with the AO-based DFT softwares, its
MPI/OpenMP hybrid parallelization scheme, as well as the installation and
usage guidelines. Such information is highly useful for the potential users
and developers of this software.

By combining the LRI framework and the spatial locality of the AO basis
functions, LibRPA enables quadratic-scaling computation of the KS response
function matrix, the essential quantity behind the RPA calculations. Further
taking into account of the decaying behavior of the Green’s function, LibRPA
achieves linear scaling for very large systems. LibRPA is massively
parallelized with the hybrid MPI/OpenMP scheme, and test calculations show its
nearly ideal scalability with respect to the computational resources. To
demonstrate the practical usability of LibRPA, we applied it to calculate the
adsorption energy and equilibrium height of a H2O molecule adsorbed on
graphene. The results obtained validate the accuracy and reliability of
LibRPA, showcasing its capability in handling realistic computational
materials science problems.

LibRPA is still under active development, and its capabilities are being
continously expanded. For example, we are currently working on the low-scaling
GW implementation in LibRPA which supports efficient calculation of the quasi-
particle energies. Details about the algorithm and implementation of the GW
module in LibRPA will be described in a future publication. In short, by
utilizing the locality of the AO basis functions and the Green’s functions, as
well as efficient parallelization schemes, LibRPA is being developed into a
useful tool for performing RPA calculations and beyond.

##  Appendix A Flow diagrams of the parallel schemes

The flow diagrams of the parallel distribution schemes of LibRPA are shown in
Fig. 4. The left and right panels illustrate the two distribution schemes
based on atomic pairs ${\\{\mathcal{U},\mathcal{V}\\}}$ and the lattice
vector-imaginary time pairs ${\\{{\mathbf{R}},\tau\\}}$, respectively. In the
first scheme, no communication is needed for the Fourier transform from the
real-space imaginary-time domain to the ${\mathbf{k}}$-space imaginary-
frequency domain, but heavy communication is necessary for the LU
decomposition of a globally distributed response function matrix. On the other
hand, in the second scheme, the Fourier transform requires communications
across different processes, but the LU decomposition can be performed locally
in each process.

(a) atom-pair leading

(b) $({\mathbf{R}},\tau)$-leading

Figure 4: Illustration of the parallel distribution schemes implemented in
LibRPA: the atom-pair leading parallel scheme (left) and the
$({\mathbf{R}},\tau)$-leading parallel scheme (right). The distribution of
subtasks across processes using different strategies
(${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\\{\mathcal{U},\mathcal{V}\\}}$
or
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\\{{\mathbf{R}},\tau\\}}$),
which is automatically chosen by LibRPA depending on the type of the system to
be calculated. The yellow boxes indicate the points at which inter-process
communication occurs.

##  Appendix B API functions of LibRPA

API functions are provided to allow data exchange between the DFT code and
LibRPA at the memory level. This approach eliminates the need for large disk
files for input data and enables post-processing of LibRPA output within a
single execution of the DFT program. The API functions are categorized into
three main groups: environment setup, input parsing and computation. Listing 1
gives a list of API functions as an example.

Listing 1: API functions of LibRPA.

⬇

void initialize_librpa_environment(

MPI_Comm comm_global_in, int is_fortran_comm,

int redirect_stdout, const char *output_filename);

void finalize_librpa_environment();

void set_dimension(int nspins, int nkpts, int nstates, int nbasis, int
natoms);

void set_wg_ekb_efermi(int nspins, int nkpts, int nstates, double* wg, double*
ekb, double efermi);

void set_ao_basis_wfc(int ispin, int ik, double* wfc_real, double* wfc_imag);

void set_latvec_and_G(double lat_mat[9], double G_mat[9]);

void set_kgrids_kvec_tot(int

nk1, int nk2, int nk3, double* kvecs);

void set_ibz2bz_index_and_weight(const int nk_irk, const int* ibz2bz_index,
const double* wk_irk);

void set_ao_basis_aux(int I, int J, int nbasis_i, int nbasis_j, int naux_mu,
int* R, double* Cs_in, int insert_index_only);

void set_aux_bare_coulomb_k_atom_pair(int ik, int I, int J, int naux_mu, int
naux_nu, double* Vq_real_in, double* Vq_imag_in);

void set_aux_bare_coulomb_k_2D_block(int ik, int max_naux, int mu_begin, int
mu_end, int nu_begin, int nu_end, double* Vq_real_in, double* Vq_imag_in);

void get_rpa_correlation_energy(double *rpa_corr, double
*rpa_corr_irk_contrib);

##  Appendix C Correcting the BSSE in adsorption energy calculations

In RPA calculations of adsorption energies using AO basis sets, a significant
challenge is the basis set superposition errors (BSSEs). This error arises
from the artificial stabilization due to overlapping basis sets of interacting
fragments. While the counterpoise (CP) correction scheme is standard for
cluster systems, such a correction scheme is not implemented in FHI-aims for
bulk systems. As a workaround, we estimate the magnitude of BSSE based on a
cluster model, where the $\ce{H2O}$ molecule is adsorbed on the a graphene
cluster containing 27 C atoms. We check that further increasing the cluster
size does not lead to a noticeable change of the value of the BSSE. This
scheme provides a pragmatic estimate of the BSSE for periodic bulk systems.
The results presented in Fig. 3 and Table 1 are corrected for the BSSE using
the above-described approach.

## References

  * [1] D. Bohm, D. Pines, A collective description of electron interactins: Iii. coulomb interaction in a degenerate electron gas, Phys. Rev. 92 (1953) 609.
  * [2] M. Gell-Mann, K. A. Brueckner, Correlation energy of an electron gas at high density, Phys. Rev. 106 (1957) 364.
  * [3] J. Hubbard, The description of collective motions in terms of many-body perturbation theory. ii. the correlation energy of a free-electron gas, Proc. Roy. Soc. (London) 243 (1234) (1958) 336–352.
  * [4] F. Furche, Molecular tests of the random phase approximation to the exchange-correlation energy functional, Phys. Rev. B 64 (2001) 195120.
  * [5] M. Fuchs, X. Gonze, Accurate density functionals: Approaches using the adiabatic-connection fluctuation-dissipation theorem, Phys. Rev. B 65 (2002) 235109.
  * [6] H. Eshuis, J. E. Bates, F. Furche, Electron correlation methods based on the random phase approximation, Theor. Chem. Acc. 131 (2012) 1084.
  * [7] X. Ren, P. Rinke, C. Joas, M. Scheffler, Random-phase approximation and its applications in computational chemistry and materials science, J. Mater. Sci. 47 (2012) 7447.
  * [8] D. C. Langreth, J. P. Perdew, Exchange-correlation energy of a metal surface: Wave-vector analysis, Phys. Rev. B 15 (1977) 2884.
  * [9] O. Gunnarsson, B. I. Lundqvist, Exchange and correlation in atoms, molecules, and solids by the spin-density-functional formalism, Phys. Rev. B 13 (1976) 4274.
  * [10] J. P. Perdew, K. Schmidt, Jacob’s ladder of density functional approximations for the exchange-correlation energy, in: V. Van Doren, C. Van Alsenoy, P. Geerlings (Eds.), Density Functional Theory and its Application to Materials , AIP, Melville, NY, 2001.
  * [11] P. Hohenberg, W. Kohn, Phys. Rev. 136 (1964) B864.
  * [12] W. Kohn, L. J. Sham, Phys. Rev. 140 (1965) A1133.
  * [13] X. Ren, P. Rinke, M. Scheffler, Exploring the random phase approximation: Application to co adsorbed on cu(111), Phys. Rev. B 80 (2009) 045402.
  * [14] L. Schimka, J. Harl, A. Stroppa, A. Grüneis, M. Marsman, F. Mittendorfer, G. Kresse, Accurate surface and adsorption energies from many-body perturbation theory, Nat. Mater. 9 (2010) 741.
  * [15] J. Paier, X. Ren, P. Rinke, G. E. Scuseria, A. Grüneis, G. Kresse, M. Scheffler, New J. Phys. 14 (2012) 043002.
  * [16] X. Ren, P. Rinke, G. E. Scuseria, M. Scheffler, Renormalized second-order perturbation theory for the electron correlation energy: Concept, implementation, and benchmarks, Phys. Rev. B 88 (2013) 035120.
  * [17] S. Lebègue, J. Harl, T. Gould, J. G. Ángyán, G. Kresse, J. F. Dobson, Cohesive properties and asymptotics of the dispersion interaction in graphite by the random phase approximation, Phys. Rev. Lett. 105 (2010) 196401.
  * [18] Y. Zhang, J. W. Furness, B. Xiao, J. Sun, Subtlety of TiO2 phase stability: Reliability of the density functional theory predictions and persistence of the self-interaction error, J. Chem. Phys. 150 (2019) 014105.
  * [19] N. Sengupta, J. E. Bates, A. Ruzsinszky, From semilocal density functionals to random phase approximation renormalized perturbation theory: A methodological assessment of structural phase transitions, Phys. Rev. B 97 (2018) 235136.
  * [20] C. Cazorla, T. Gould, Polymorphism of bulk boron nitride, Sci. Adv. 5 (1) (2019) eaau5832.
  * [21] S. Yang, X. Ren, Phase stability of the argon crystal: first-principles study based on random phase approximation plus renormalized single excitation corrections, New J. Phys. 24 (2022) 033049.
  * [22] J. Chedid, N. Jocelyn, H. Eshuis, Energies, structures, and harmonic frequencies of small water clusters from the direct random phase approximation, J. Chem. Phys. 155 (8) (2021) 084303.
  * [23] M. N. Tahir, T. Zhu, H. Shang, J. Li, V. Blum, X. Ren, Localized resolution of identity approach to the analytical gradients of random-phase approximation ground-state energy: algorithm and benchmarks, J. Chem. Theory Comput. 18 (2022) 5297.
  * [24] I. Y. Zhang, X. Xu, Simultaneous attenuation of both self-interaction error and nondynamic correlation error in density functional theory: A spin- pair distinctive adiabatic-connection approximation, J. Phys. Chem. Lett. 10 (2019) 2617.
  * [25] H. Eshuis, J. E. Bates, F. Furche, Electron correlation methods based on the random phase approximation, Theor. Chem. Acc. 131 (2012) 1–18.
  * [26] X. Ren, P. Rinke, V. Blum, J. Wieferink, A. Tkatchenko, A. Sanfilippo, K. Reuter, M. Scheffler, Resolution-of-identity approach to Hartree–Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions, New J. Phys. 14 (5) (2012) 053020.
  * [27] D. Rocca, Z. Bai, R.-C. Li, G. Galli, A block variational procedure for the iterative diagonalization of non-hermitian random-phase approximation matrices, J. Chem. Phys. 136 (3) (2012).
  * [28] H.-V. Nguyen, S. de Gironcoli, Efficient calculation of exact exchange and rpa correlation energies in the adiabatic-connection fluctuation-dissipation theory, Phys. Rev. B 79 (2009) 205114.
  * [29] J. Harl, G. Kresse, Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipartion theory, Phys. Rev. B 77 (2008) 045136.
  * [30] D. Lu, H.-V. Nguyen, G. Galli, Power series expansion of the random phase approximation correlation energy: The role of the third- and higher-order contributions, J. Chem. Phys. 133 (2010) 154110.
  * [31] J. Harl, G. Kresse, Accurate bulk properties from approximate many-body techniques, Phys. Rev. Lett. 103 (5) (2009) 056401.
  * [32] T. Olsen, K. S. Thygesen, Random phase approximation applied to solids, molecules, and graphene-metal interfaces: From van der Waals to covalent bonding, Phys. Rev. B 87 (7) (2013) 075111.
  * [33] D. Rocca, Random-phase approximation correlation energies from lanczos chains and an optimal basis set: Theory and applications to the benzene dimer, J. Chem. Phys. 140 (18) (2014).
  * [34] H. Eshuis, J. Yarkony, F. Furche, Fast computation of molecular random phase approximation correlation energies using resolution of the identity and imaginary frequency integration, J. Chem. Phys. 132 (2010) 234114.
  * [35] M. Del Ben, O. Schütt, T. Wentz, P. Messmer, J. Hutter, J. VandeVondele, Enabling simulation at the fifth rung of DFT: Large scale RPA calculations with excellent time to solution, Comput. Phys. Commun. 187 (2015) 120–129.
  * [36] D. Neuhauser, E. Rabani, R. Baer, Expeditious stochastic calculation of random-phase approximation energies for thousands of electrons in three dimensions, J. Phys. Chem. Lett. 4 (2013) 1172.
  * [37] J. E. Moussa, Cubic-scaling algorithm and self-consistent field for the random-phase approximation with second-order screened exchange, J. Chem. Phys. 140 (1) (2014) 014107.
  * [38] M. Kaltak, J. Klimeš, G. Kresse, Cubic scaling algorithm for the random phase approximation: Self-interstitials and vacancies in si, J. Chem. Theory Comput. 10 (2014) 2498.
  * [39] M. Kállay, Linear-scaling implementation of the direct random-phase approximation, J. Chem. Phys. 142 (20) (2015) 204105.
  * [40] J. Wilhelm, P. Seewald, M. Del Ben, J. Hutter, Large-scale cubic-scaling random phase approximation correlation energy calculations using a gaussian basis, J. Chem. Theory Comput. 12 (12) (2016) 5851–5859.
  * [41] D. Graf, M. Beuerle, H. F. Schurkus, A. Luenser, G. Savasci, C. Ochsenfeld, Accurate and efficient parallel implementation of an effective linear-scaling direct random phase approximation method, J. Chem. Theory Comput. 14 (2018) 2505.
  * [42] A. Luenser, H. F. Schurkus, C. Ochsenfeld, Vanishing-overhead linear-scaling random phase approximation by cholesky decomposition and an attenuated coulomb-metric, J. Chem. Theory Comput. 13 (4) (2017) 1647–1655.
  * [43] J. Lu, K. Thicke, Cubic scaling algorithms for rpa correlation using interpolative separable density fitting, J. Comput. Phys. 351 (15) (2017) 187–202.
  * [44] I. Duchemin, X. Blase, Separable resolution-of-the-identity with all-electron gaussian bases: Application to cubic-scaling rpa, J. Chem. Phys. 150 (17) (2019) 174120.
  * [45] H. N. Rojas, R. W. Godby, R. J. Needs, Space-time method for ab initio calculations of self-energies and dielectric response functions of solids, Phys. Rev. Lett. 74 (1995) 1827.
  * [46] I. D. White, R. W. Godby, M. M. Rieger, R. J. Needs, Dynamic image potential at an al(111) surface, Phys. Rev. Lett. 80 (1998) 4265.
  * [47] M. M. Rieger, L. Steinbeck, I. White, H. Rojas, R. Godby, The GW space-time method for the self-energy of large systems, Comput. Phys. Commun. 117 (3) (1999) 211–228.
  * [48] P. Merlot, T. Kjrgaard, T. Helgaker, R. Lindh, F. Aquilante, S. Reine, T. B. Pedersen, Attractive electron-electron interactions within robust local fitting approximations, J. Comput. Chem. 34 (2013) 1486.
  * [49] L. N. Wirz, S. S. Reine, T. B. Pedersen, On resolution-of-the-identity electron repulsion integral approximations and variational stability, J. Chem. Theory Comput 13 (2017) 4897.
  * [50] D. S. Hollman, H. F. Schaefer, E. F. Valeev, Semi-exact concentric atomic density fitting: Reduced cost and increased accuracy compared to standard density fitting, J. Chem. Phys. 140 (6) (2014) 064109.
  * [51] X. Wang, C. A. Lewis, E. F. Valeev, Efficient evaluation of exact exchange for periodic systems via concentric atomic density fitting, J. Chem. Phys. 153 (12) (2020) 124116.
  * [52] S. V. Levchenko, X. Ren, J. Wieferink, R. Johanni, P. Rinke, V. Blum, M. Scheffler, Hybrid functionals for large periodic systems in an all-electron, numeric atom-centered basis framework, Comp. Phys. Comm. 192 (2015) 60.
  * [53] P. Lin, X. Ren, L. He, Accuracy of localized resolution of the identity in periodic hybrid functional calculations with numerical atomic orbitals, J. Phys. Chem. Lett. 11 (8) (2020) 3082–3088.
  * [54] P. Lin, X. Ren, L. He, Efficient Hybrid Density Functional Calculations for Large Periodic Systems Using Numerical Atomic Orbitals, J. Chem. Theory Comput. 17 (1) (2021) 222–239.
  * [55] R. Shi, P. Lin, M.-Y. Zhang, L. He, X. Ren, Subquadratic-scaling real-space random phase approximation correlation energy calculations for periodic systems with numerical atomic orbitals, Phys. Rev. B 109 (3) (2024) 035103.
  * [56] S. Kokott, F. Merz, Y. Yao, C. Carbogno, M. Rossi, V. Havu, M. Rampp, M. Scheffler, V. Blum, Efficient all-electron hybrid density functionals for atomistic simulations beyond 10 000 atoms, J. Chem. Phys. 161 (2) 024112.
  * [57] X. Ren, F. Merz, H. Jiang, Y. Yao, M. Rampp, H. Lederer, V. Blum, M. Scheffler, All-electron periodic ${G}_{0}{W}_{0}$ implementation with numerical atomic orbital basis functions: Algorithm and benchmarks, Phys. Rev. Materials 5 (2021) 013807.
  * [58] E. Spadetto, P. H. T. Philipsen, A. Förster, L. Visscher, Toward pair atomic density fitting for correlation energies with benchmark accuracy, J. Chem. Theory Comput. 19 (5) (2023) 1499–1516.
  * [59] A. C. Ihrig, J. Wieferink, I. Y. Zhang, M. Ropo, X. Ren, P. Rinke, M. Scheffler, V. Blum, Accurate localized resolution of identity approach for linear-scaling hybrid density functionals and for many-body perturbation theory, New J. Phys. 17 (2015) 093020.
  * [60] I. Y. Zhang, A. J. Logsdail, X. Ren, S. V. Levchenko, L. Ghiringhelli, M. Scheffler, Test set for materials science and engineering with user-friendly graphic tools for error analysis: Systematic benchmark of the numerical and intrinsic errors in state-of-the-art electronic-structure approximations, New J. Phys. 21 (2019) 013025.
  * [61] M. Kaltak, J. Klimeš, G. Kresse, Low scaling algorithms for the random phase approximation: Imaginary time and laplace transformations, J. Chem. Theory Comput. 10 (6) (2014) 2498–2507.
  * [62] M. Kaltak, J. Klimeš, G. Kresse, Cubic scaling algorithm for the random phase approximation: Self-interstitials and vacancies in si, Phys. Rev. B 90 (2014) 054115.
  * [63] CP2K github repository: https://github.com/cp2k/cp2k.
  * [64] GreenX github repository: https://github.com/nomad-coe/greenX.
  * [65] M. Azizi, J. Wilhelm, D. Golze, M. Giantomassi, R. L. Panadés-Barrueta, F. A. Delesma, A. Buccheri, A. Gulans, P. Rinke, C. Draxl, et al., Time-frequency component of the GreenX library: minimax grids for efficient RPA and GW calculations, J. Open Source Softw. 8 (90) (2023) 5570.
  * [66] LibRI github repository: https://github.com/abacusmodeling/LibRI.
  * [67] V. Blum, F. Hanke, R. Gehrke, P. Havu, V. Havu, X. Ren, K. Reuter, M. Scheffler, Ab-initio molecular simulations with numeric atom-centered orbitals, Comp. Phys. Comm. 180 (2009) 2175.
  * [68] P. Li, X. Liu, M. Chen, P. Lin, X. Ren, L. Lin, C. Yang, L. He, Large-scale ab initio simulations based on systematically improvable atomic basis, Comput. Mater. Sci. 112 (2016) 503.
  * [69] LibRPA github repository: https://github.com/Srlive1201/LibRPA.
  * [70] T. O. Wehling, A. I. Lichtenstein, M. I. Katsnelson, First-principles studies of water adsorption on graphene: The role of the substrate, Appl. Phys. Lett. 93 (20) (2008).
  * [71] I. Hamada, Adsorption of water on graphene: A van der Waals density functional study, Phys. Rev. B 86 (19) (2012) 195436.
  * [72] M. Sacchi, A. Tamtögl, Water adsorption and dynamics on graphene and other 2D materials: computational and experimental advances, ADV PHYS-X 8 (1) (2023) 2134051.
  * [73] L. Yang, L. Zhang, X. Jiao, Y. Qiu, W. Xu, The electrochemical performance of reduced graphene oxide prepared from different types of natural graphites, RSC Adv. 11 (7) (2021) 4042–4052.
  * [74] P. K. Sahoo, N. Kumar, A. Jena, S. Mishra, C.-P. Lee, S.-Y. Lee, S.-J. Park, Recent progress in graphene and its derived hybrid materials for high-performance supercapacitor electrode applications, RSC Adv. 14 (2) (2024) 1284–1303.
  * [75] J. Ma, A. Michaelides, D. Alfè, L. Schimka, G. Kresse, E. Wang, Adsorption and diffusion of water on graphene from first principles, Phys. Rev. B 84 (2011) 033402.
  * [76] J. G. Brandenburg, A. Zen, M. Fitzner, B. Ramberger, G. Kresse, T. Tsatsoulis, A. Grüneis, A. Michaelides, D. Alfè, Physisorption of Water on Graphene: Subchemical Accuracy from Many-Body Electronic Structure Methods, J. Phys. Chem. Lett 10 (3) (2019) 358–368.