# A machine learning model to classify dynamic processes in liquid water

Jie Huang  Jie Huang, Prof. Dr. Shiben Li
Department of Physics, Wenzhou University, Wenzhou, Zhejiang 325035, China;
E-mail<EMAIL_ADDRESS>Gang Huang  Dr. Gang Huang
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190,
China; E-mail<EMAIL_ADDRESS>Shiben Li  Jie Huang, Prof. Dr. Shiben Li
Department of Physics, Wenzhou University, Wenzhou, Zhejiang 325035, China;
E-mail<EMAIL_ADDRESS>
###### Abstract

The dynamics of water molecules plays a vital role in understanding water. We
combined computer simulation and deep learning to study the dynamics of
H-bonds between water molecules. Based on _ab initio_ molecular dynamics
simulations and a newly defined directed Hydrogen (H-) bond population
operator, we studied a typical dynamic process in bulk water: interchange, in
which the H-bond donor reverses roles with the acceptor. By designing a
recurrent neural network-based model, we have successfully classified the
interchange and breakage processes in water. We have found that the ratio
between them is approximately 1:4, and it hardly depends on temperatures from
280 to 360 K. This work implies that deep learning has the great potential to
help distinguish complex dynamic processes containing H-bonds in other
systems.

††preprint: AIP/123-QED

##  I Introduction

As one of the big questions in the 21st century (1), the structure of water is
essential for understanding cells, biological processes, and ecosystems (2, 3,
4, 5). Water’s surprising properties (6, 7, 8), such as increased density on
melting, high surface tension, maximum density at 4 ∘C, are closely related to
the H-bonds(9, 10, 11). Despite the fact that it is tough to capture the
ultrafast motion of atoms during dynamic processes (12), watching water
molecules as they dance is the key to understand the dynamic properties of
water (13) from the molecular level. Many methods over the past three decades
were used to study water molecules’ motion, such as scanning tunneling
microscopy (STM) (14, 15), femtosecond pump-probe (16), infrared (IR)
spectroscopy (17, 18, 12, 19), X-rays (20, 13, 21), neutron scattering (22),
and computer simulations (23, 24, 9, 25).

In this work, we focus on one specific dynamic process in bulk water:
interchange (26), in which the H-bond donor reverses roles with the acceptor
in the same H-bond. This process was observed in the gas-phase water dimer by
Saykally and coworkers. The interchange process, which involves the quantum
tunneling effect (27, 26, 23, 15, 25), is essential for understanding water
molecules’ dynamics. Also, since interchange processes are closely related to
the H-bond network dynamics, it is likely to play a critical role in
biological processes, like proton transfer (28, 29, 30). So far, interchange
processes have been found in water dimer adsorbed on metal surfaces (15).
Using _ab initio_ molecular dynamics (AIMD) simulations (31), Ranea et al.
(23) found that the interchange process can be used to explain the rapid
diffusion behavior of water dimer on the Pd(111) surface. Fang et al. (25)
found that interchange process is a mechanism of the rapid movement of water
dimers on metal surfaces. As for bulk water, Lagge and Hynes found that the
redirection of water molecules involves large-angle jumps (24), which involves
the redirection of multiple water molecules referred to as H-bond exchange,
and it is supported by the subsequent experiments (32, 33).

The interchange process involves the concerted rotation of both water
molecules engaged in a H-bonded pair. This mechanism is important in small
clusters where the future hydrogen-bond donor OH group is typically initially
dangling. There are some simulation studies on the interchange process in
water clusters (27, 26, 34, 35, 36), as far as we know, the question of the
ratio of interchange to other dynamic processes related to H-bonds in bulk
water has not been discussed. To determine the proportion of interchange
processes, we simulated bulk water in a canonical (NVT) ensemble using a
specific AIMD simulation method: the density functional molecular dynamics
(DFTMD) simulation. We observed interchange processes in bulk water by
analyzing the dynamic trajectory.

As it’s tough to quantify interchange processes in a large number of ultrafast
dynamic processes in liquid water, we have designed a recurrent neural network
(RNN)-based model to classify the H-bond dynamic processes. Unlike general
classification methods, this model has the capability of classifying the
_dynamic processes_ related to H-bonds in bulk water. Using this model, we
have obtained the relative ratio of interchange and breakage processes in bulk
water and explored the effect of temperature on this ratio. Our work presents
the great capacity to use the RNN-based deep learning method to study the
dynamic properties of liquid water.

The aims of this work is to provide a machine learning-based model to classify
dynamic processes and to determine the proportion of interchange processes in
water as one useage of the model. The organization of the paper is the
following. We present the results and discussion in Sec. II. At first, the
dynamic graph representation of H-bond networks is introduced in II.1 and the
main characteristics of interchange processes are obtained in II.2. Then we
implement the RNN-based classifier for different types of dynamic processes in
liquid water in II.3 and explore the temperature dependence of the relative
ratios of interchange and breakage processes in II.4; The discussion of two
factors, the mean number of H-bonds and the rate of breakage and formation of
H-bonds in liquid water, related closely to the temperature dependence are
discussed in II.5. Finally, we present the methods details and conclusions of
our study in Sec. III and IV, respectively.

##  II Results and discussion

###  II.1 Dynamic graph representation of H-bond networks

As shown in Fig. 1, a directed _dynamic graph_ is used to describe the bulk
water system of $N$ water molecules. Each water molecule may form an H-bond
with any of the remaining $N$-$1$ molecules. For convenience, we call any pair
of water molecules $(i,j)$ a _quasi_ -hydrogen bond (Q-bond), denoted as
$b_{ij}$ and represented as a dashed line in Fig. 1.

Figure 1: Dynamic graph representation of the H-bond network in simulated bulk
water. Nodes represent water molecules; solid red or green arrows represent
H-bonds; and dashed grey lines represent Q-bonds. The colors red, grey, and
green indicate $\tilde{h}_{ij}$=$1$, $\tilde{h}_{ij}$=$0$, and
$\tilde{h}_{ij}$=$-1$, respectively. From the time sequence of
$\tilde{h}_{ij}$, we know how the H-bond configuration of $b_{ij}$ changes
over time. Four typical H-bond configuration change processes are illustrated
for $b_{13}$, $b_{15}$, $b_{45}$, and $b_{24}$, corresponding to interchange,
breakage, formation, and no change, respectively.

Inspired by Luzar and Chandler’s H-bond population operator (37), we define a
_directed_ H-bond population operator $\tilde{h}_{ij}$ for $b_{ij}\ (i<j)$ at
time $t$ as Eq. 4.

$\displaystyle\tilde{h}_{ij}(t)=\left\\{\begin{array}[]{rcl}1&&{\text{H-bonded,
$i$ is the donor}}\\\ 0&&{\text{Not H-bonded}}\\\ -1&&{\text{H-bonded, $j$ is
the donor}}\end{array}\right.$ (4)

We know from $\tilde{h}_{ij}$ whether an H-bond exists in $b_{ij}$ and the
donor-acceptor pair of the formed H-bond. At the bottom of Fig. 1, we
demonstrate four typical H-bond configuration change processes by using the
sequences of $\tilde{h}$: interchange, breakage, formation, and no change.
Besides, the Q-bonds likely to form H-bonds are the most relevant water
molecule pairs to the breakage and reforming of H-bond networks. The following
geometric criteria (38, 39, 40) of an H-bond is used: O-O distance
$R_{\text{OO}^{\prime}}<R_{\text{cutoff}}=3.5\text{ \AA}$ and angle
$\text{O-H}\cdots\text{O}>\theta_{\text{cutoff}}=120^{\circ}$.

Figure 2:  Scheme of the geometric coordinates. $R_{\text{OO}^{\prime}}$ is
the O-O distance. Four angles $\widehat{\mathrm{O}\mathrm{H1O^{\prime}}}$,
$\widehat{\mathrm{O}\mathrm{H2O^{\prime}}}$ ,
$\widehat{\mathrm{O^{\prime}}\mathrm{H^{\prime}1O}}$, and
$\widehat{\mathrm{O^{\prime}}\mathrm{H^{\prime}2O}}$ are represented as
$\theta_{\text{a}}$, $\theta_{\text{b}}$, $\theta_{\text{c}}$, and
$\theta_{\text{d}}$, respectively. If $R_{\text{OO}^{\prime}}<3.5\text{ \AA}$
, and any angle $\theta>120^{\circ}$
($\theta\in\\{\theta_{\text{a}}$,$\theta_{\text{b}}$, $\theta_{\text{c}}$,
$\theta_{\text{d}}\\}$), then an H-bond exists in this Q-bond. Here, the
oxygen atom O as a donor donates the hydrogen atom H2 to the acceptor
$\text{O}^{\prime}$. Since $R_{\text{OO}^{\prime}}<3.5\text{ \AA}$ and
$\theta_{\text{b}}>120^{\circ}$, we describe this state of
$b_{\text{OO}^{\prime}}$ at this time $t$ by
$\tilde{h}_{\text{OO}^{\prime}}(t)=1$.

As shown in Fig. 2, $R_{\text{OO}^{\prime}}$, $\theta_{\text{a}}$,
$\theta_{\text{b}}$, $\theta_{\text{c}}$, and $\theta_{\text{d}}$ are
monitored for Q-bonds to study the reorientation and breakage mechanism of
H-bonds.

###  II.2 Interchange process

Figure 3:  Interchange ($\text{I}_{1}$), breakage ($\text{I}_{2}$), and
bifurcation rearrangement ($\text{I}_{3}$) process for one typical Q-bond in
bulk water. When an H-bond exists in a Q-bond if
$\theta_{\text{a}}>\theta_{\text{cutoff}}$ or
$\theta_{\text{b}}>\theta_{\text{cutoff}}$, then the oxygen atom O is the
donor; else, if $\theta_{\text{c}}>\theta_{\text{cutoff}}$ or
$\theta_{\text{d}}>\theta_{\text{cutoff}}$, then the oxygen atom O′ is the
donor. Three typical processes are interchange, where the water molecule pairs
exchange their roles as H-bond donor and acceptor; breakage, where the H-bond
is breaking as the distance increase of this water molecule pair; and
bifurcation rearrangement, where the donated hydrogen atom of the H-bond donor
exchanged. Through $\tilde{h}$, we can see whether an H-bond exists between a
Q-bond, also know the donor and acceptor if an H-bond exists. In panel (C),
the grey, red, and green lines indicate the $\tilde{h}_{{\text{OO}}^{\prime}}$
states.

The AIMD simulation trajectory allows us to observe the details of the H-bond
dynamics. Figure 3 demonstrates the dynamics of the distance, angles, and
directed H-bond population for $b_{\text{OO}^{\prime}}$. Intervals
$\text{I}_{1}$, $\text{I}_{2}$, and $\text{I}_{3}$ correspond to three typical
H-bond dynamic processes. interchange ($\text{I}_{1}$): We notice
$\theta_{\text{b}}>\theta_{\text{cutoff}}$ in the first half and
$\theta_{\text{d}}>\theta_{\text{cutoff}}$ in the second half. Besides,
$\tilde{h}_{\text{OO}^{\prime}}$ changes from $1$ to $-1$, indicating that the
donor and acceptor have exchanged. _breakage_ ($\text{I}_{2}$):
$\tilde{h}_{\text{O}\text{O}^{\prime}}=-1$ in the first half of
$\text{I}_{2}$, and $\tilde{h}_{\text{O}\text{O}^{\prime}}=0$ for most of the
second half. There is no H-bond in the second half because
$R_{\text{OO}^{\prime}}>R_{\text{cutoff}}$, i.e., the increase of distance
$R_{\text{OO}^{\prime}}$ causes the H-bond to break. _bifurcation
rearrangement_ motion ($\text{I}_{3}$)(26, 41): At first, the hydrogen atom
$\text{H}^{\prime}1$ is donated to form an H-bond as
$\theta_{\text{c}}>\theta_{\text{cutoff}}$. Then $\theta_{c}$ decreases and
$\theta_{d}$ increases until $\theta_{\text{d}}>\theta_{\text{cutoff}}$, i.e.,
the other hydrogen atom $\text{H}^{\prime}2$ of the donor is donated to form
the H-bond. Therefore, the hydrogen atom contributed by the donor is changed.
Because of the identity of hydrogen atoms, it is impossible to distinguish the
configuration of water molecules before and after the process. However, during
the interchange process, the direction of the water molecules’ dipole moment
will change, indicating that the water molecules’ microscopic configuration
will change. So in the rest of the article, we focus on the interchange
process.

Figure 4: A typical interchange process, where two water molecules exchange
their roles as H-bond donor and acceptor via water molecules’ reorientation in
an concerted manner. The donor oxygen atom has changed from the original O to
O′ (color of dashed line changed from red to green). Besides, we have also
noticed that the H-bond briefly breaks during the interchange process, causing
the fluctuation of the $\tilde{h}$ sequence.

Figure 4 shows a typical interchange process in water (see Fig. S4 in SI Sec.
3 and movies in supplementary material for more H-bond configuration change
processes). A dashed line represents an H-bond, and its color (red or green)
indicates its direction. Using $\tilde{h}$, we can describe the H-bond
configuration change progress without paying attention to the distance and
angles. Therefore, $\tilde{h}$ dramatically simplifies the description for the
H-bond configuration change process. Nevertheless, during dynamic processes,
the fluctuations of $\tilde{h}$ that can result from the vibration of water
molecules will bring a huge challenge for the classification of H-bond
configuration change processes. In addition, due to a large number of Q-bonds
in the simulated bulk water, finding a specific H-bond configuration change
process in 60 ps is like finding a needle in a haystack. Therefore, we design
an RNN-based model that recognizes the dynamic processes related to H-bonds
and uses it to determine various processes in water, thereby determining the
ratio of interchanges.

###  II.3 RNN-based classifier for H-bond configuration change process

We can see the interchange and breakage processes intuitively from
$\tilde{h}$. Specifically, in the interchange process, $\tilde{h}$ changes
from $\pm 1$ to $\mp 1$; in the breakage process, $\tilde{h}$ changes from
$\pm 1$ to 0. Therefore, in principle, by observing the sequence of
$\tilde{h}$ within a time window, we can classify the H-bond configuration
changes during this period. Although we can see some change patterns in
interchange and breakage processes, it is still challenging to distinguish
different $\tilde{h}$ sequences due to the fluctuation. Therefore, we have
designed a processing flow to classify the H-bond configuration change process
based on RNN, as shown in Fig. 5.

Figure 5: The processing flow of the H-bond configuration change classifier
based on RNN. (i). Different types of $\tilde{h}$ sequences: $\text{T}_{1}$:
Formation or no H-bond; $\text{T}_{2}$, $\text{T}_{3}$: No change;
$\text{T}_{4}$: Negative sequence; $\text{T}_{5}$, $\text{T}_{6}$: Diffusion;
$\text{T}_{7}$, $\text{T}_{8}$: interchange. We refer to the sequences of
breakage and interchange as positive sequences. (ii). The preprocessor filters
out the high-frequency components of $\tilde{h}$ and excludes $\text{T}_{1}$,
$\text{T}_{2}$, and $\text{T}_{3}$. (iii). The classifier consists of a BLSTM
AE to separate the positive and negative sequences and a final classifier to
distinguish breakage and interchange sequences.

In the preprocessor, we use a low-pass filter to filter out the high-frequency
fluctuations of $\tilde{h}$ sequences. As we focus on the configuration change
processes of H-bonds, we exclude sequences without H-bonds at the beginning
($\text{T}_{1}$) and the sequences whose H-bond configuration are unchanged
($\text{ T}_{2}$, $\text{T}_{3}$) according to the initial value and the
variance of $\tilde{h}$ sequences (see Methods section). After preprocessing,
the task we need to deal with is a time series classification problem: In
addition to interchange and breakage processes, there are also many irregular
and complicated processes. We call the sequences of interchange and breakage
positive and all sequences other than these two types negative for
convenience. Negative sequences ($\text{T}_{4}$) do not have any particular
pattern. We do not expect that general supervised learning can be used to
distinguish them. Nevertheless, we can teach a machine to learn to recognize
positive sequences. Due to the need to classify time series, we use a typical
method for modeling ordered data (42, 43, 44), recurrent neural network (RNN)
(45, 46). Specifically, we have designed a bidirectional long short-term
memory (BLSTM) autoencoder (AE), whose goal is to reconstruct the input
sequences as much as possible. We have trained this AE using positive
sequences only and evaluated how well the AE reconstructs for an input
sequence $\mathbf{x}$ using reconstruction error $\mathcal{L}(\mathbf{x})$
(see Methods Section III.3, Eq. 8). After training, the autoencoder can
reconstruct positive processes very well. However, when we input negative
sequences into the AE, likely, it would not be able to reconstruct them well,
leading to the reconstruction errors of these negative sequences greater than
that of the positive sequences. Through the reconstruction error, we can
determine whether a $\tilde{h}$ sequence is positive or negative.

Figure 6: (A) Densities of reconstruction error $\mathcal{L}$ for interchange,
breakage, and negative sequences. (i) BLSTM AE can reconstruct positive
sequences well. Hence, the reconstruction errors for interchange and breakage
sequences are relatively small, mainly less than $\mathcal{L}_{\text{T}}$.
(ii) Since negative sequences are not used to train BLSTM AE, it is much more
difficult for the autoencoder to reconstruct them. Therefore, the
reconstruction errors are relatively large, mainly greater than
$\mathcal{L}_{\text{T}}$. (iii) Once $\mathcal{L}_{\text{T}}$ is determined,
we use it as the threshold to distinguish positive and negative sequences. (B)
Densities of the range $\delta$ for interchange and breakage sequences. The
two densities are significantly different from each other.

Finally, we use a final classifier to distinguish sequences between
interchange and breakage processes from positive sequences. We use the range
of a positive sequence $\mathbf{x}$ to determine whether it is interchange or
breakage, which is defined as
$\delta(\mathbf{x})=\max{\mathbf{x}}-\min{\mathbf{x}}$.

Figure 6 (A) shows the densities of the reconstruction errors for interchange,
breakage, and negative sequences. Since BLSTM AE can reconstruct positive
sequences well, the reconstruction errors of interchange and breakage
sequences are small, most of which are smaller than the reconstruction error
threshold $\mathcal{L}_{\text{T}}$ ($\mathcal{L}_{\text{T}}$ determination and
corresponding accuracy analysis are described in SI Sec. 3, Fig. S3). Negative
sequences are not used to train the autoencoder, so it is much more difficult
to reconstruct them. Hence, the reconstruction errors are relatively large,
most of which are greater than $\mathcal{L}_{\text{T}}$. As long as we find a
suitable reconstruction error threshold, we can get a classifier for positive
and negative sequences. Figure 6 (B) shows the densities for the range of
normalized interchange and breakage sequences. The two distributions are
significantly different from each other. Therefore, the final classifier can
distinguish interchange and breakage sequences very well via
$\delta_{\text{T}}=0.66$, as shown in the dashed line (see the classification
process in SI Sec. 3, Fig. S4-S5 ). Therefore, we have obtained an H-bond
configuration change classifier based on an RNN autoencoder.

###  II.4 Proportions of interchange at different temperatures

Figure 7: The number (A) and proportion (B) of interchange and breakage
processes determined by the RNN-based classifier at different temperatures.
(i) With the temperature increasing, the number of interchange and breakage
processes increases first and then decreases on the whole. (ii) The relative
ratio of interchange to breakage basically does not depend on temperature.

To explore the effect of temperature on the H-bond configuration change
process, we have simulated nine bulk water systems containing $N=64$ water
molecules. The temperature ranges from 280 to 360 K every 10 K. Using the RNN-
based model, we classify $\tilde{h}$ sequence, count the number of interchange
and breakage sequences at each temperature. As shown in Fig. 7, the number of
interchange and breakage processes shows a "rising first, then decreasing"
trend as the temperature increases. In other words, there is an overall upward
trend from 280 to 330 K. However, as the temperature continues to rise, the
number of detected interchange and breakage processes tends to decrease. As we
use the method of width-fixed sliding window, the absolute number of
interchanges and breakages would change along the step size of sliding window.
These numbers would increase as we decrease the step size. Therefore, we
focused on the trend of the detected processes over the temperatures. On the
other hand, although the number of interchange and breakage processes vary at
different temperatures, the relative ratios between the two are almost
unchanged, which is still about 1:4 (see SI Sec.4 for the step size effect of
the sliding window). This result indicates that the relative ratio is almost
not dependent on temperature, and the interchange process is another important
mechanism in bulk water besides the breakage process. Next, we will explain
this trend of the number of interchange and breakage processes from the
following two aspects: the number of H-bond per molecule and the change rate
of the coarse-grained H-bond network configuration.

###  II.5 The trend of interchange and breakage process number

To understand the trend in Fig. 7 (A), we first calculate the number of
H-bonds per molecule ($n_{\text{HB}}$) in the simulated system. At time $t$,
$n_{\text{HB}}$ can be expressed as Eq. 5,

$\displaystyle
n_{\text{HB}}(t)=\frac{2}{N}\sum_{i=1}^{N}\sum_{j>i}^{N}|\tilde{h}_{ij}(t)|$
(5)

where $N=64$ is the number of water molecules in bulk water systems, and
$|\tilde{h}_{ij}(t)|$ is the absolute value of $\tilde{h}_{ij}(t)$, i.e., the
H-bond direction is ignored. The factor 2 is derived from the fact that one
H-bond in water is shared by two water molecules. For a certain trajectory at
one temperature, by counting $n_{\text{HB}}$ at each time $t$, we get the
distribution of $n_{\text{HB}}$ (one density plot in Fig. 8 A). Then we use an
$L$-dimensional vector $\bm{\tilde{h}}(t)$ to represent the coarse-grained
H-bond network configuration for the simulated bulk water system at time $t$
in Eq. 6,

$\displaystyle\bm{\tilde{h}}(t)=(\tilde{h}_{12}(t),\tilde{h}_{13}(t),\cdots,\tilde{h}_{ij}(t),\cdots,\tilde{h}_{N-1,N}(t))$
(6)

where $L=N(N-1)/2$ is the number of Q-bonds in the system. So in a unit time,
we get a set $H$ of $\bm{\tilde{h}}(t)$ in Eq. 7,

$\displaystyle H=\\{\bm{\tilde{h}}(t)\mid t=t_{0}+k\Delta t,k=0,1,\cdots,M\\}$
(7)

where $t_{0}$ represents the start time of the unit time window, $\Delta t$ is
the time interval between two adjacent frames, and $M$ is the length of the
unit time window. In a unit time $t_{\text{w}}=M\Delta t$, the number of graph
configuration can be expressed as $\Omega=|H|$, where $|H|$ is the size of the
set $H$, i.e., the number of different $\bm{\tilde{h}}$ vectors in this unit
time. The number $\Omega$ of graph configuration per unit time characterizes
the rate of breakage and reforming of the H-bonds in bulk water. The
theoretical upper bound of $\Omega$ in $t_{\text{w}}$ is $M+1$; in this case,
all $\bm{\tilde{h}}$ vectors are different. By changing the start point
$t_{0}$, we get the distribution of $\Omega$.

Figure 8:  The temperature dependence of (A) The distributions of the number
$n_{\text{HB}}$ of H-bonds per molecule. (B) The distributions of $\ln\Omega$
characterizing the rate of H-bond breakage and reforming. The dashed line
denotes the upper bound of $\ln\Omega$ in the unit time of 1 ps.

Figure 8 shows the temperature dependence of the distributions of
$n_{\text{HB}}$ and $\ln\Omega$. The width of a density plot indicates the
probability of $n_{\text{HB}}$ or $\ln\Omega$ at the corresponding
temperature. From the medians (white dots) of violin plots in Fig. 8 (B), we
see $\Omega$ is relatively smaller at lower temperatures, indicating fewer
changes of H-bond configuration in the unit time. This result explains why the
number of interchange and breakage processes at lower temperatures in Fig. 7
(A) are smaller. Besides, the direct reason for the decrease in the number of
interchange and breakage processes at higher temperatures is that thermal
motions tend to break H-bonds (thus reducing $n_{\text{HB}}$). Therefore, the
number of interchange and breakage processes in Fig. 7 (A) is determined by
$n_{\text{HB}}$ and $\Omega$ together.

##  III Methods

###  III.1 AIMD simulations

AIMD simulations were carried out for bulk water of 64 water molecules within
the canonical NVT ensemble using CP2K/QUICKSTEP (v7.1) (47). The number $N$ of
water molecules was 64 for all bulk water systems at different temperatures
from 280 to 360 K. The length of the periodic cubic box was 12.4295 Å. The
discretized integration time step $\Delta t$ was set to 0.5 fs. The simulation
time was 60 ps. The BLYP functional, which consists of Becke non-local
exchange (48) and Lee-Yang-Parr correlation (49), was used; Interactions
between the valence electrons and the ionic cores were described by GTH
pseudopotentials (50, 51); Valence electrons were expanded in a basis set
consisting of double-zeta Gaussian functions (52) and plane waves with a
cutoff energy of 280 Ry (47). The Nos$\acute{\text{e}}$-Hoover chain
thermostat (53) was used to conserve temperature. DFT-D3 correction (54) for
the dispersion interaction was used to obtain a more accurate description of
the vibrational properties. It is worth mentioning that the analysis method we
proposed can be used on various simulation data, and the AIMD simulation used
here is one of the options. The graph-based analysis method is independent of
simulation data, so this method can be used to analyze more accurate
simulation data in the future.

###  III.2 Sequence collection and preprocessing

The sequence length of $\tilde{h}$ was 200 corresponding to 8 ps simulation
time. Positive sequences in which only one interchange or breakage process
occurred were collected. Negative $\tilde{h}$ sequences used to evaluate the
BLSTM AE classifier were also collected. BLSTM AE was trained by 6786 positive
sequences, of which the interchange and breakage processes each accounted for
half (754 positive sequences at each temperature). There were 18,931 negative
sequences for evaluating the BLSTM AE classifier. The filtered sequence
$\tilde{h}_{f}[n]$ was obtained by second-order Butterworth filter implemented
by Scipy (55). In addition, if $\tilde{h}_{f}[0]-0.5<0.15$, indicating no
H-bond at the beginning ($\text{ T}_{1}$). If the standard deviation $\sigma$
of $\tilde{h}_{f}[n]$ satisfy $\sigma<0.1$, then we consider the H-bond
configuration in the Q-bond has not changed ($\text{ T}_{2}$, $\text{
T}_{3}$).

###  III.3 Bidirectional LSTM autoencoder classifier

The encoder and the decoder of BLSTM AE can be expressed as two
transformations, $\phi:\mathcal{X}\rightarrow\mathcal{F}$ and
$\psi:\mathcal{F}\rightarrow\mathcal{X}$, where $\mathcal{X}$ and
$\mathcal{F}$ are the input space and the feature space, respectively. The
dimension of $\mathcal{F}$ is smaller than that of $\mathcal{X}$, and the
feature vector $\phi(\mathbf{x})$ is the compressed representation of input
$\mathbf{x}$. The input $\mathbf{x}$ of BLSTM AE is the normalized and
filtered directed H-bond population operator sequence $\tilde{h}_{f}[n]$. The
reconstruction error of BLSTM AE for a sequence $\mathbf{x}=\tilde{h}_{f}[n]$
is defined as

$\displaystyle\mathcal{L}_{\omega,\omega^{\prime}}\left(\mathbf{x}\right)=\left\|\mathbf{x}-\psi_{\omega^{\prime}}\left(\phi_{\omega}\left(\mathbf{x}\right)\right)\right\|^{2}$
(8)

where $\omega$, $\omega^{\prime}$ represent the parameters of the encoder and
decoder respectively. The purpose of training is to obtain the optimal
$\omega$, $\omega^{\prime}$,

$\displaystyle\omega^{*},\omega^{\prime*}=\arg\min_{\omega,\omega^{\prime}}\frac{1}{m}\sum_{i=1}^{m}\mathcal{L}_{\omega,\omega^{\prime}}(\mathbf{x}^{i})$
(9)

where $\mathbf{x}^{i}$ represents the $i$-th sequence (SI, Fig. S1-S2).

##  IV Summary

In summary, we have designed and trained a deep learning-based model to
recognize different types of processes related to H-bonds. The priority of
this model are its remarkable ability to classify different dynamic processes
of water molecules and its wide range of applications to different kinds of
simulation methods. The model can be transfered to other dynamic systems
containning H-bonds with the form of $\text{O-H}\cdots\text{O}$. As a feasible
example, combined with AIMD simulations, we have found that the relative ratio
of interchange and breakage processes in bulk water is approximately 1:4, and
this ratio hardly depends on temperature.

Moreover, the key concepts used in this work are the dynamic graph and the
newly defined directed H-bond population. This reasonable coarse-grained
description of the H-bond network simplifies the analysis of H-bond dynamics
dramatically. This work demonstrates that the semi-supervised RNN-based model
has an outstanding capability of classifying the dynamic processes related to
H-bonds in bulk water, which implies the great potential to extend our present
scheme to distinguish more complex dynamic processes in other systems like the
water-vapor interface and electrolyte solutions.

###### Acknowledgements.

This research was supported by the National Natural Science Foundation of
China (NSFC) (Grant No. 21973070) and the Graduate Scientific Research
Foundation of Wenzhou University. The simulations were performed on the
cluster in the College of Mathematics and Physics at Wenzhou University.

## References

  * Kennedy (2005) D. Kennedy, “What don’t we know?” Science 309, 75–75 (2005).
  * Franks (2000) F. Franks,  _Water: a matrix of life_ , 2nd ed. (Royal Society of Chemistry, 2000).
  * Pal and Zewail (2004) S. K. Pal and A. H. Zewail, “Dynamics of water in biological recognition,” Chemical Reviews 104, 2099–2124 (2004).
  * Chaplin (2006) M. Chaplin, “Do we underestimate the importance of water in cell biology?” Nature Reviews Molecular Cell Biology 7, 861–866 (2006).
  * Ball (2017) P. Ball, “Water is an active matrix of life for cell and molecular biology,” Proceedings of the National Academy of Sciences 114, 13327–13335 (2017).
  * Stillinger (1980) F. H. Stillinger, “Water revisited,” Science 209, 451–457 (1980).
  * Errington and Debenedetti (2001) J. R. Errington and P. G. Debenedetti, “Relationship between structural order and the anomalies of liquid water,” Nature 409, 318–321 (2001).
  * Dumé (2020) I. Dumé, “Second critical point appears in two models of water,” Physics World 33, 7i–7i (2020).
  * Kumar, Schmidt, and Skinner (2007) R. Kumar, J. R. Schmidt, and J. L. Skinner, “Hydrogen bonding definitions and dynamics in liquid water,” The Journal of Chemical Physics 126, 204107 (2007).
  * Nilsson and Pettersson (2015) A. Nilsson and L. G. M. Pettersson, “The structural origin of anomalous properties of liquid water,” Nature Communications 6 (2015), 10.1038/ncomms9998.
  * Wilhelmsen and Deutsches Elektronen-Synchrotron (2020) U. Wilhelmsen and D. Deutsches Elektronen-Synchrotron,  _The strangest liquid in the world: water amazes scientists time and again_, edited by T. Mundzeck, Vol. 20 (Deutsches Elektronen Synchrotron, DESY, Hamburg, 2020) p. 44 pages.
  * Karamatskos  _et al._ (2019) E. T. Karamatskos, S. Raabe, T. Mullins, A. Trabattoni, P. Stammer, G. Goldsztejn, R. R. Johansen, K. Długołecki, H. Stapelfeldt, M. J. J. Vrakking, S. Trippel, A. Rouzée, and J. Küpper, “Molecular movie of ultrafast coherent rotational dynamics of OCS,” Nature Communications 10 (2019), 10.1038/s41467-019-11122-y.
  * Perakis  _et al._ (2018) F. Perakis, G. Camisasca, T. J. Lane, A. Späh, K. T. Wikfeldt, J. A. Sellberg, F. Lehmkühler, H. Pathak, K. H. Kim, K. Amann-Winkel, S. Schreck, S. Song, T. Sato, M. Sikorski, A. Eilert, T. McQueen, H. Ogasawara, D. Nordlund, W. Roseker, J. Koralek, S. Nelson, P. Hart, R. Alonso-Mori, Y. Feng, D. Zhu, A. Robert, G. Grübel, L. G. M. Pettersson, and A. Nilsson, “Coherent x-rays reveal the influence of cage effects on ultrafast water dynamics,” Nature Communications 9 (2018), 10.1038/s41467-018-04330-5.
  * Mitsui (2002) T. Mitsui, “Water diffusion and clustering on pd(111),” Science 297, 1850–1852 (2002).
  * Kumagai  _et al._ (2008) T. Kumagai, M. Kaizu, S. Hatta, H. Okuyama, T. Aruga, I. Hamada, and Y. Morikawa, “Direct observation of hydrogen-bond exchange within a single water dimer,” Physical Review Letters 100 (2008), 10.1103/physrevlett.100.166101.
  * Woutersen, Emmerichs, and Bakker (1997) S. Woutersen, U. Emmerichs, and H. J. Bakker, “Femtosecond mid-ir pump-probe spectroscopy of liquid water: Evidence for a two-component structure,” Science 278, 658–660 (1997), https://science.sciencemag.org/content/278/5338/658.full.pdf .
  * Bakker and Nienhuys (2002) H. J. Bakker and H.-K. Nienhuys, “Delocalization of protons in liquid water,” Science 297, 587–590 (2002).
  * Fecko (2003) C. J. Fecko, “Ultrafast hydrogen-bond dynamics in the infrared spectroscopy of water,” Science 301, 1698–1702 (2003).
  * ichi Inoue  _et al._ (2020) K. ichi Inoue, M. Ahmed, S. Nihonyanagi, and T. Tahara, “Reorientation-induced relaxation of free OH at the air/water interface revealed by ultrafast heterodyne-detected nonlinear spectroscopy,” Nature Communications 11 (2020), 10.1038/s41467-020-19143-8.
  * Iwashita  _et al._ (2017) T. Iwashita, B. Wu, W.-R. Chen, S. Tsutsui, A. Q. R. Baron, and T. Egami, “Seeing real-space dynamics of liquid water through inelastic x-ray scattering,” Science Advances 3, e1603079 (2017).
  * Loh  _et al._ (2020) Z.-H. Loh, G. Doumy, C. Arnold, L. Kjellsson, S. H. Southworth, A. A. Haddad, Y. Kumagai, M.-F. Tu, P. J. Ho, A. M. March, R. D. Schaller, M. S. B. M. Yusof, T. Debnath, M. Simon, R. Welsch, L. Inhester, K. Khalili, K. Nanda, A. I. Krylov, S. Moeller, G. Coslovich, J. Koralek, M. P. Minitti, W. F. Schlotter, J.-E. Rubensson, R. Santra, and L. Young, “Observation of the fastest chemical processes in the radiolysis of water,” Science 367, 179–182 (2020).
  * Head-Gordon and Hura (2002) T. Head-Gordon and G. Hura, “Water structure from scattering experiments and simulation,” Chemical Reviews 102, 2651–2670 (2002).
  * Ranea  _et al._ (2004) V. A. Ranea, A. Michaelides, R. Ramírez, P. L. de Andres, J. A. Vergés, and D. A. King, “Water dimer diffusion on pd$\\{$111$\\}$ assisted by an h-bond donor-acceptor tunneling exchange,” Physical Review Letters 92 (2004), 10.1103/physrevlett.92.136104.
  * Laage and Hynes (2006) D. Laage and J. T. Hynes, “A Molecular Jump Mechanism of Water Reorientation,” Science 311, 832–835 (2006).
  * Fang  _et al._ (2020) W. Fang, J. Chen, P. Pedevilla, X.-Z. Li, J. O. Richardson, and A. Michaelides, “Origins of fast diffusion of water dimers on surfaces,” Nature Communications 11 (2020), 10.1038/s41467-020-15377-8.
  * Keutsch and Saykally (2001) F. N. Keutsch and R. J. Saykally, “Water clusters: Untangling the mysteries of the liquid, one molecule at a time,” Proceedings of the National Academy of Sciences 98, 10533–10540 (2001).
  * Fellers  _et al._ (1999) R. S. Fellers, C. Leforestier, L. B. Braly, M. G. Brown, and R. J. Saykally, “Spectroscopic determination of the water pair potential,” Science 284, 945–948 (1999).
  * Agmon (1995) N. Agmon, “The grotthuss mechanism,” Chemical Physics Letters 244, 456–462 (1995).
  * Thomaston  _et al._ (2017) J. L. Thomaston, R. A. Woldeyes, T. Nakane, A. Yamashita, T. Tanaka, K. Koiwai, A. S. Brewster, B. A. Barad, Y. Chen, T. Lemmin, M. Uervirojnangkoorn, T. Arima, J. Kobayashi, T. Masuda, M. Suzuki, M. Sugahara, N. K. Sauter, R. Tanaka, O. Nureki, K. Tono, Y. Joti, E. Nango, S. Iwata, F. Yumoto, J. S. Fraser, and W. F. DeGrado, “XFEL structures of the influenza m2 proton channel: Room temperature water networks and insights into proton conduction,” Proceedings of the National Academy of Sciences 114, 13357–13362 (2017).
  * Gelenter  _et al._ (2021) M. D. Gelenter, V. S. Mandala, M. J. M. Niesen, D. A. Sharon, A. J. Dregni, A. P. Willard, and M. Hong, “Water orientation and dynamics in the closed and open influenza b virus m2 proton channels,” Communications Biology 4, 338 (2021).
  * Kühne  _et al._ (2020) T. D. Kühne, M. Iannuzzi, M. D. Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack, and J. Hutter, “CP2k: An electronic structure and molecular dynamics software package - quickstep: Efficient and accurate electronic structure calculations,” The Journal of Chemical Physics 152, 194103 (2020).
  * Moilanen  _et al._ (2009) D. E. Moilanen, D. Wong, D. E. Rosenfeld, E. E. Fenn, and M. D. Fayer, “Ion–water hydrogen-bond switching observed with 2d ir vibrational echo chemical exchange spectroscopy,” Proceedings of the National Academy of Sciences 106, 375–380 (2009).
  * Ji, Odelius, and Gaffney (2010) M. Ji, M. Odelius, and K. J. Gaffney, “Large angular jump mechanism observed for hydrogen bond exchange in aqueous perchlorate solution,” Science 328, 1003–1005 (2010).
  * Schulz  _et al._ (2018) R. Schulz, Y. von Hansen, J. O. Daldrop, J. Kappler, F. Noé, and R. R. Netz, “Collective hydrogen-bond rearrangement dynamics in liquid water,” The Journal of Chemical Physics 149, 244504 (2018).
  * Samala and Agmon (2019) N. R. Samala and N. Agmon, “Thermally induced hydrogen-bond rearrangements in small water clusters and the persistent water tetramer,” ACS Omega 4, 22581–22590 (2019).
  * Méndez and Laria (2020) E. Méndez and D. Laria, “Nuclear quantum effects on the hydrogen bond donor–acceptor exchange in water–water and water–methanol dimers,” The Journal of Chemical Physics 153, 054302 (2020).
  * Luzar and Chandler (1996) A. Luzar and D. Chandler, “Hydrogen-bond kinetics in liquid water,” Nature 379, 55–57 (1996).
  * Sciortino and Fornili (1989) F. Sciortino and S. L. Fornili, “Hydrogen bond cooperativity in simulated water: Time dependence analysis of pair interactions,” The Journal of Chemical Physics 90, 2786–2792 (1989).
  * Balasubramanian, Pal, and Bagchi (2002) S. Balasubramanian, S. Pal, and B. Bagchi, “Hydrogen-bond dynamics near a micellar surface: Origin of the universal slow relaxation at complex aqueous interfaces,” Physical Review Letters 89 (2002), 10.1103/physrevlett.89.115505.
  * Michaud-Agrawal  _et al._ (2011) N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein, “MDAnalysis: A toolkit for the analysis of molecular dynamics simulations,” Journal of Computational Chemistry 32, 2319–2327 (2011).
  * Brown, Keutsch, and Saykally (1998) M. G. Brown, F. N. Keutsch, and R. J. Saykally, “The bifurcation rearrangement in cyclic water clusters: Breaking and making hydrogen bonds,” The Journal of Chemical Physics 109, 9645–9647 (1998).
  * Hughes  _et al._ (2019) T. W. Hughes, I. A. D. Williamson, M. Minkov, and S. Fan, “Wave physics as an analog recurrent neural network,” Science Advances 5 (2019), 10.1126/sciadv.aay6946.
  * Rank  _et al._ (2020) N. Rank, B. Pfahringer, J. Kempfert, C. Stamm, T. Kühne, F. Schoenrath, V. Falk, C. Eickhoff, and A. Meyer, “Deep-learning-based real-time prediction of acute kidney injury outperforms human predictive performance,” npj Digital Medicine 3 (2020), 10.1038/s41746-020-00346-8.
  * Tsai, Kuo, and Tiwary (2020) S.-T. Tsai, E.-J. Kuo, and P. Tiwary, “Learning molecular dynamics with simple language model built upon long short-term memory neural network,” Nature Communications 11 (2020), 10.1038/s41467-020-18959-8.
  * Hopfield (1982) J. J. Hopfield, “Neural networks and physical systems with emergent collective computational abilities,” Proceedings of the National Academy of Sciences 79, 2554–2558 (1982).
  * Hochreiter and Schmidhuber (1997) S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural Computation 9, 1735–1780 (1997).
  * VandeVondele  _et al._ (2005) J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing, and J. Hutter, “Quickstep: Fast and accurate density functional calculations using a mixed gaussian and plane waves approach,” Comput. Phys. Commun. 167, 103–128 (2005).
  * Becke (1988) A. D. Becke, “Density-functional exchange-energy approximation with correct asymptotic behavior,” Phys. Rev. A 38, 3098 (1988).
  * Lee, Yang, and Parr (1988) C. Lee, W. Yang, and R. G. Parr, “Development of the colic-salvetti correlation-energy formula into a functional of the electron density,” Phys. Rev. B 37, 785 (1988).
  * Hartwigsen, Goedecker, and Hutter (1998) C. Hartwigsen, S. Goedecker, and J. Hutter, “Relativistic separable dual-space gaussian pseudopotentials from H to Rn,” Phys. Rev. B 58, 3641–3662 (1998).
  * G. Lippert and Parrinello (1999) J. H. G. Lippert and M. Parrinello, “The gaussian and augmented-plane-wave density functional method for ab initio molecular dynamics simulations,” Theor. Chem. Acc. 103, 124 (1999).
  * VandeVondele and Hutter (2007) J. VandeVondele and J. Hutter, “Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases,” The Journal of Chemical Physics 127, 114105 (2007).
  * Martyna, Klein, and Tuckerman (1992) G. J. Martyna, M. L. Klein, and M. Tuckerman, “Nosé-hoover chains: The canonical ensemble via continuous dynamics,” J. Chem. Phys. 97, 2635 (1992).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, “A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu,” J. Chem. Phys. 132, 154104 (2010).
  * Virtanen  _et al._ (2020) P. Virtanen, , R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, and P. van Mulbregt, “SciPy 1.0: fundamental algorithms for scientific computing in python,” Nature Methods 17, 261–272 (2020).


# Efficient Hartree-Fock Exchange Algorithm with Coulomb Range Separation and
Long-Range Density Fitting

Qiming Sun<EMAIL_ADDRESS>Quantum Engine LLC, Washington 98516, US

###### Abstract

Separating the Coulomb potential into short-range and long-range components
enables the use of different electron repulsion integral algorithms for each
component. The short-range part can be efficiently computed using the
analytical algorithm due to the locality in both Gaussian-type orbital basis
and the short-range Coulomb potentials. The integrals for the long-range
Coulomb potential can be approximated with the density fitting method. A very
small auxiliary basis is sufficient for the density fitting method to
accurately approximate the long-range integrals. This feature significantly
reduces the computational efforts associated with the $N^{4}$ scaling in
density fitting algorithms. For large molecules, the range separation and
long-range density fitting method outperforms the conventional analytical
integral evaluation scheme employed in Hartree-Fock calculations and provides
more than twice the overall performance. Additionally, this method yields
higher accuracy compared to regular density fitting methods. The error in the
Hartree-Fock energy can be easily reduced to 0.1 $\mu E_{h}$ per atom, which
is significantly more accurate than the typical error of 10 $\mu E_{h}$ per
atom observed in regular density fitting methods.

##  I Introduction

A well-known bottleneck in Hartree-Fock and Kohn-Sham density functional
theory is the computation of the two-electron repulsion integrals (ERIs)

$(\mu\nu|\kappa\lambda)=\iint\mu(\mathbf{r}_{1})\nu(\mathbf{r}_{1})\frac{1}{r_{12}}\kappa(\mathbf{r}_{2})\lambda(\mathbf{r}_{2})d^{3}\mathbf{r}_{1}d^{3}\mathbf{r}_{2}$
(1)

for four atomic orbitals basis functions $\mu(\mathbf{r})$, $\nu(\mathbf{r})$,
$\kappa(\mathbf{r})$, $\lambda(\mathbf{r})$. To improve the computation speed
of ERIs, the density fitting (DF) or resolution of identity methods were
developed to approximate ERIs Whitten (1973); Dunlap, Connolly, and Sabin
(1979); Vahtras, Almlöf, and Feyereisen (1993); Eichkorn  _et al._ (1995);
Dunlap (2000). This technique has been widely and successfully applied in
quantum chemistry simulations for both molecules and crystalste Velde and
Baerends (1991); Jaffe and Hess (1996); Čársky, Čurík, and Varga (2012);
Burow, Sierka, and Mohamed (2009); Feyereisen, Fitzgerald, and Komornicki
(1993); Komornicki and Fitzgerald (1993); Manby, Knowles, and Lloyd (2001);
Jung  _et al._ (2005); Hollman, Schaefer, and Valeev (2014); Sun  _et al._
(2017); Ye and Berkelbach (2021); Lee  _et al._ (2022). Despite the great
success that DF has achieved, it still faces several challenges in terms of
computation costs and accuracy.

The DF method exhibits notable performance advantages when computing Hartree-
Fock exchange (HFX) in small and medium-size systems. However, these
advantages become less pronounced for larger molecules with more than 10,000
basis functions. When applying the DF method to calculate HFX, the integral
problem is transformed into a tensor contraction problem with a scaling of
$N^{4}$. The performance advantage of the DF method can be attributed to the
computational benefits provided by high-performance tensor contraction
software, such as MKL, OpenBLAS, BLIS, Cutlass, as well as hardware such as
GPUs. For sufficiently large systems, the conventional analytical ERI
computation for HFX may actually be faster than DF HFX. Due to the locality of
Gaussian-type orbital (GTO) basis and the sparsity of the GTO integral tensor,
the conventional method has an asymptotic scaling between $N^{1}$ and $N^{2}$.
Besides the scaling wall, another challenge that restricts the application of
the DF method to large systems is the issue of I/O and storage. The DF
integral tensor is a three index tensor, with two indices representing the
orbital basis and one index representing the auxiliary basis. Typically, the
number of auxiliary basis functions is more than three times the number of
orbital basis functions. In systems with 10,000 atomic orbitals, the size of a
DF tensor can quickly reach tens or even hundreds of terabytes. Even
considering the sparsity of the DF tensor, the compressed tensor can still
occupy a significant amount of storage.

The effectiveness of the DF approximation heavily relies on the quality of the
auxiliary basis. To find a balance between computational cost and accuracy,
the auxiliary basis must be optimized in advance for the specific orbital
basis and computation tasks. Different computation tasks, such as the
computation of Coulomb energy (denoted as JFIT on basis exchangePritchard  _et
al._ (2019)), the Coulomb and exchange terms (JKFIT), the post-HF methods, and
relativistic methods, have different requirements on the auxiliary basis. If
an inappropriate auxiliary basis is used, the DF approximation may result in
reduced accuracy, or increased computational costs, or both. To address this
issue, various approaches have been proposed to automate the construction of
auxiliary basis setsYang, Rendell, and Frisch (2007); Boström  _et al._
(2009); Aquilante  _et al._ (2009); Kállay (2014); Stoychev, Auer, and Neese
(2017); Lehtola (2021, 2023); De Chavez and Lindh (2023). The automatic basis
generation procedure can simplify the selection of auxiliary basis sets in
calculations, although it may slightly increases computational costs.

In our previous workSun (2023a), we developed the algorithm of ERI computation
in crystal systems with the separation of short-range (SR) and long-range (LR)
in Coulomb potential

$\frac{1}{r_{12}}=\frac{\mathrm{erfc}(\omega
r_{12})}{r_{12}}+\frac{\mathrm{erf}(\omega r_{12})}{r_{12}}.$ (2)

SR ERIs are computed analytically in real space, while LR ERIs are computed
numerically on reciprocal grids. In this decomposition treatment, we observed
that the rank of the LR ERI tensor is significantly lower than that of the
full-range ERIs. Depending on the value of the parameter $\omega$ in the
range-separation operators, the rank of LR ERIs can be reduced to 1/100 of the
rank of full-range ERIs. This feature can be leveraged to construct low-rank
DF integral tensors, which effectively solves the problems of auxiliary basis
dependence, I/O pressure and performance penalty in regular DF methods. In
this work, we will test the combination of range-separation treatment and the
long-range density fitting (LRDF) ERIs for the computation of Coulomb and HFX
terms in molecules.

##  II Theory

###  II.1 Auxiliary Gaussian basis for long-range density fitting

The range separation of the Coulomb potential enables us to decompose the ERI
tensor into SR and LR components,

$(\mu\nu|\kappa\lambda)=(\mu\nu|\kappa\lambda)_{\text{SR}}+(\mu\nu|\kappa\lambda)_{\text{LR}}.$
(3)

The SR ERI tensor contains details of the ERI structure. Its rank is similar
to that of the full-range ERI tensor. The LR ERI tensor is a low-rank tensor.
It is efficient to use DF technique to construct this tensor with auxiliary
basis $\phi_{P}(\mathbf{r})$

$(\mu\nu|\kappa\lambda)_{\text{LR}}=\sum_{PQ}(\mu\nu|\phi_{P})_{\text{LR}}[M^{-1}]_{PQ}(\kappa\lambda|\phi_{Q})_{\text{LR}}.$
(4)

where

$\displaystyle(\mu\nu|\phi_{P})_{LR}=\iint\mu(\mathbf{r}_{1})\nu(\mathbf{r}_{1})\frac{\mathrm{erf}(\omega
r_{12})}{r_{12}}\phi_{P}(\mathbf{r}_{2})d^{3}\mathbf{r}_{1}d^{3}\mathbf{r}_{2}.$
(5) $\displaystyle
M_{PQ}=\iint\phi_{P}(\mathbf{r}_{1})\frac{\mathrm{erf}(\omega
r_{12})}{r_{12}}\phi_{Q}(\mathbf{r}_{2})d^{3}\mathbf{r}_{1}d^{3}\mathbf{r}_{2}.$
(6)

By performing pivoted Cholesky decomposition (CD)Lehtola (2019) or eigenvalue
decomposition (ED) on the metric matrix $\mathbf{M}$,

$\mathbf{M}=\mathbf{U}\lambda\mathbf{U}^{\dagger},$ (7)

we can obtain the long-range density fitting (LRDF) tensor $\mathbf{L}$ with
reduced auxiliary dimension

$L_{\mu\nu,P}=\sum_{Q}\lambda_{Q}^{-1/2}U_{PQ}(\mu\nu|\phi_{Q})_{\text{LR}}.$
(8)

The transformed coefficients $\mathbf{U}$ can be truncated to form a short and
wide matrix, with fewer rows than columns. It only needs to retain the rows
that correspond to non-zero values of $\lambda$.

A straightforward choice for the auxiliary basis in LRDF would be the regular
DF auxiliary basis that was previously optimized for the full-range Coulomb
potential. However, these types of auxiliary basis sets are typically over-
complete and too large, compared to the rank of the LR ERI tensor. Since the
LR Coulomb potential screens the interactions between nearby density
distributions, the auxiliary basis for LRDF only needs to ensure an accurate
description of the effective potential experienced by the remote electron
density distributions, which are the products of two orbital basis functions
$\mu(\mathbf{r})\nu(\mathbf{r})$. This requirement is much lower than the
requirements for the auxiliary basis in regular DF methods. Therefore, a
smaller auxiliary basis with enough flexibility to represent potentials in the
multipole expansion is sufficient for the decomposition of LR ERIs. The
numerical tests in Section III suggest that primitive Gaussian basis functions
up to angular momentum $l=2$ with exponent of 1.0 on each atom are enough to
represent the LR potential. In practice, using only $s$-type Gaussian
functions for Hydrogen atoms and Gaussian functions with $l=0..2$ for the
remaining atoms can achieve an accuracy of micro-Hartree per atom.

###  II.2 Computing LR ERIs with PW basis

The LR ERI can be Fourier transformed to an integral in reciprocal space

$(\mu\nu|\kappa\lambda)_{\text{LR}}=\frac{1}{(2\pi)^{3}}\int\frac{4\pi}{G^{2}}e^{-\frac{G^{2}}{4\omega^{2}}}\rho_{\mu\nu}(\mathbf{G})\rho_{\lambda\kappa}(\mathbf{G})^{*}d\mathbf{G}^{3}.$
(9)

Eq. (9) is essentially a one-particle Gaussian integrals. The integrand
$\rho_{\mu\nu}(\mathbf{G})$ can be obtained through a Fourier transformation
of the orbital product $\rho_{\mu\nu}(\mathbf{r})$

$\rho_{\mu\nu}(\mathbf{G})=\int
e^{-i\mathbf{G}\cdot\mathbf{r}}\rho_{\mu\nu}(\mathbf{r})d\mathbf{r}^{3}=e^{-\frac{G^{2}}{4\alpha_{\mu\nu}}}e^{-i\mathbf{G}\cdot\mathbf{R}_{\mu\nu}}P_{\mu\nu}(\mathbf{G}).$
(10)

Here, $P_{\mu\nu}(\mathbf{G})$ is a polynomial associated with the Fourier
transform of the polynomial part in Gaussian product theoremBesalú and Carbó-
Dorca (2011). By changing the coordinates system to spherical coordinates,

$d^{3}\mathbf{G}=G^{2}\sin\theta d\theta d\varphi,$ (11)

we can transform the integral (9) to

$\frac{1}{4\pi^{2}}\int
e^{-\frac{G^{2}}{4\omega^{2}}}\rho_{\mu\nu}(G,\theta,\varphi)\rho_{\lambda\kappa}(G,\theta,\varphi)^{*}\sin\theta
dGd\theta d\varphi.$ (12)

Assuming that we have already evaluated the angular part of the integrals, the
remaining integral for the radial part $G$ can be expressed formally as

$\displaystyle\int_{0}^{\infty}e^{-\frac{G^{2}}{4\theta}}S(G)P(G)dG,$ (13)
$\displaystyle\theta=(\alpha_{\mu\nu}^{-1}+\alpha_{\kappa\lambda}^{-1}+\omega^{-2})^{-1},$
(14) $\displaystyle\alpha_{\mu\nu}=\alpha_{\mu}+\alpha_{\nu}.$ (15)

$S(G)$ is a factor associated to the coordinates of the four Gaussian orbital
basis functions.

$\displaystyle S(G)\sim\int
e^{-i\mathbf{G}|\mathbf{R}_{\mu\nu}-\mathbf{R}_{kl}|}d\theta d\varphi,$ (16)
$\displaystyle\mathbf{R}_{\mu\nu}=\frac{\alpha_{\mu}\mathbf{R}_{\mu}+\alpha_{\nu}\mathbf{R}_{\nu}}{\alpha_{\mu\nu}}.$
(17)

The other terms in the radial part are collected in the polynomial $P(G)$.

When the value of $\omega$ is selected such that $\omega^{2}\ll\alpha_{ij}$
and $\omega^{2}\ll\alpha_{kl}$, the exponential term in the integral (13) can
be approximated as $e^{-u^{2}}$, where $u=\frac{G}{2\omega}$. The integral can
then be numerically evaluated using Gauss-Hermite quadrature with respect to
the weight function $e^{-u^{2}}$. To generate the grids $G_{i}$ and weights
$w(G)_{i}$ for the radial part integral, we use the roots and weights of
Gauss-Hermite quadrature $u_{i}$ and $w(u)_{i}$ in the following manner

$\displaystyle G_{i}=2\omega u_{i},$ (18) $\displaystyle w(G)_{i}=2\omega
w(u)_{i}.$ (19)

The angular part of the integral involves a two-dimensional integration for
$\theta$ and $\varphi$,

$\int_{0}^{\pi}\sin\theta d\theta\int_{0}^{2\pi}d\varphi.$ (20)

It is convenient to change the variable to $t=\cos\theta$, and transform the
limits of angular integration to

$\int_{-1}^{1}dt\int_{0}^{2\pi}d\varphi.$ (21)

The numerical integration for the angular part can then be performed using
grids and weights on a spherical surface. One possible choice for the angular
grids is the Lebedev scheme, which is widely used in numerical integration for
density functional theory. It is also possible to generate the angular grids
by applying Gauss-Legendre quadratures for $t$ and $\varphi$ separately.
Gauss-Legendre quadratures are preferred in this case. They are more flexible
in the number of grids that can be used and they are often more accurate than
the Lebedev grids when using the same number of angular grids.

Combining the Gauss-Hermite quadratures for the radial part and Gauss-Legendre
quadratures for the angular part, we obtain the integration grids (planewave
basis) and weights in reciprocal space

$\displaystyle G_{x}$ $\displaystyle=G_{i}\sin\arccos t_{j}\cos\varphi_{k},$
$\displaystyle G_{y}$ $\displaystyle=G_{i}\sin\arccos t_{j}\sin\varphi_{k},$
$\displaystyle G_{z}$ $\displaystyle=G_{i}\cos\arccos t_{j},$ $\displaystyle
w$ $\displaystyle=\frac{1}{2\pi^{2}}G_{i}^{2}w(G)_{i}w(t)_{j}w(\varphi)_{k},$

where $w(t)_{j}$ and $w(\varphi)_{k}$ are the weights of the Gauss-Legendre
quadratures for $t$ and $\varphi$, respectively. It should be noted that the
structure factor $S(G)$ in (16) can result in a polynomial of infinite order.
It is necessary to have high resolution in both the radial and angular parts
to fully account for this effect. The number of integration grids required for
accurate integration is system dependent.

By using these integration grids and weights, the LR integral (9) can be
computed

$(\mu\nu|\kappa\lambda)_{\text{LR}}=\frac{1}{2\pi^{2}}\sum_{n}\frac{w_{n}}{|\mathbf{G}_{n}|^{2}}\rho_{\mu\nu}(\mathbf{G}_{n})\rho_{\lambda\kappa}(\mathbf{G}_{n})^{*}$
(22)

The numerical integration in reciprocal space have several advantages.

  * •

PW basis does not respond to the displacement of atoms. It makes the
computation of nuclear gradients and hessian with PW basis much simpler than
the procedure with Gaussian auxiliary fitting basisBoström  _et al._ (2013).

  * •

This method does not require I/O access. The Fourier transform (10) can be
performed on the fly. The runtime generation is more efficient than reading
data from disk.

  * •

Parallelization is straightforward. The LR potential in reciprocal space has a
diagonal structure, which allows us to isolate the FT of orbital products for
each grid point. The computation of FT on different grids can be easily
distributed to different computing nodes.

###  II.3 SCF algorithm

In the SCF program, the SR and LR parts of the Coulomb and exchange matrices
are computed separately. The SR part can be incrementally constructed during
SCF iterations using the AO-direct algorithmAlmlöf, Faegri Jr., and Korsell
(1982). Integral screening is a crucial factor for the performance of SR
integralsIzmaylov, Scuseria, and Frisch (2006); Thompson and Ochsenfeld
(2017). We employed the three hierarchical screening prescription developed in
our earlier worksSun (2023b, a). Given integral screening threshold
$\varepsilon$ and density matrix $\gamma$, the SR integral
$(\mu\nu|\kappa\lambda)_{\text{SR}}$ can be discarded if it satisfies

$\displaystyle\max(\gamma_{\mu\nu},\gamma_{\kappa\lambda},\gamma_{\mu\kappa},\gamma_{\mu\lambda},\gamma_{\nu\kappa},\gamma_{\nu\lambda})$
(23)
$\displaystyle\cdot\min(\sqrt{(\mu\nu|\mu\nu)_{\text{SR}}(\kappa\lambda|\kappa\lambda)_{\text{SR}}},$
$\displaystyle\sqrt{(\mu\mu|\kappa\kappa)_{\text{SR}}(\nu\nu|\lambda\lambda)_{\text{SR}}},$
$\displaystyle\sqrt{(\mu\mu|\lambda\lambda)_{\text{SR}}(\nu\nu|\kappa\kappa)_{\text{SR}}},$
$\displaystyle\frac{Q_{\mu\nu}Q_{\kappa\lambda}e^{-\theta R^{2}}}{R^{2}}$
$\displaystyle)<\varepsilon.$

In the above inequality, $Q_{\mu\nu}$ is a factor that behaves like the
overlap between two primitive Gaussian functions, as defined in reference Sun,
2023b. $R$ is the distance between two orbital pairs

$R=|\mathbf{R}_{\mu\nu}-\mathbf{R}_{\kappa\lambda}|.$ (24)

For the LR part HFX, we can leverage the orbital coefficients $\mathbf{C}$ and
occupation numbers $\Lambda$ to reduce computational costs

$\displaystyle L_{\mu i,P}=\sum_{\nu}L_{\mu i,P}C_{\nu i},$ (25)
$\displaystyle K_{\mu\kappa}=\sum_{iP}\Lambda_{i}L_{\mu i,P}L_{\kappa i,P}.$
(26)

The computational cost of LR HFX can be further reduced if one employs the
occ-RI-K algorithm Manzer, Epifanovsky, and Head-Gordon (2015).

##  III Numerical tests and discussions

We have implemented the LRDF approach based on the PySCF software package
v2.3.0Sun  _et al._ (2020). The SR integrals in current work are computed
using Libcint v6.0Sun (2015). Previous versions of the Libcint library can
also be used, but the performance is significantly slower for $s$-type and
$p$-type GTOs.

In the first set of tests, we focus on the impact of various factors on the
LRDF ERIs, including the auxiliary basis, the range-separation parameter
$\omega$, the threshold of ED, and the numerical integration grids in
reciprocal space. To assess the accuracy of LRDF ERIs, we compute the HF
energy using the conventional analytical 4-center ERIs, and compare this
energy with the one-shot HF energies of the LRDF methods in terms of various
parameter settings. The one-shot HF energy is computed using the density
matrix obtained from the converged HF with the analytical ERIs.

To simulate the possible scenarios in real applications, we employ two systems
in these tests: a water cluster consisting of 12 water molecules, and a linear
molecule made up of 6 glycine monomers (Gly6). They can represent two extreme
cases of basis distributions in real space. In the water cluster, orbital
basis and auxiliary basis have significant overlap, making it relatively easy
to expand the orbital product $\rho_{\mu\nu}$ with the auxiliary basis. In the
linear molecule Gly6, it is more challenging to use auxiliary basis to
describe the orbital products as the auxiliary basis is sparsely distributed
in space. With the same parameter settings, linear molecules should generally
produce larger errors. Results are collected in Table 1 \- 3.

In the second set of tests, we utilize the parameters determined from the
previous tests to evaluate the SCF performance of the LRDF method on large
molecules. The test systems consist of a water cluster with 128 water
molecules and a linear molecule with 30 glycine monomers (Gly30). Orbital
basis set cc-pVDZ is employed throughout all tests. The results of these tests
are presented in Table 4.

Table 1: The overall errors in HF energy of (H2O)12 with auxiliary Gaussian basis. The reference HF energy is -912.40220404. Auxiliary Basis | $N_{aux}$ | $\omega=0.1$ | $\omega=0.15$ | $\omega=0.2$
---|---|---|---|---
O | H |  | $N_{ED}$ | $\Delta E/\mu E_{h}$ | $N_{ED}$ | $\Delta E/\mu E_{h}$ | $N_{ED}$ | $\Delta E/\mu E_{h}$
ED threshold $10^{-4}$ |  |  |  |
$s$ | $s$ | 36 | 30 | 91.66 | 36 | 2848. | 36 | 18380.
$sp$ | $s$ | 72 | 36 | -120. | 56 | -9.97 | 70 | 899.5
$sp$ | $sp$ | 144 | 46 | -4.59 | 80 | -40.8 | 107 | -88.8
$spd$ | $s$ | 132 | 39 | -124. | 72 | -71.5 | 104 | -96.0
$spd$ | $spd$ | 324 | 49 | 0.31 | 96 | -38.0 | 152 | -57.4
weigend | 852 | 64 | -37.7 | 113 | -72.4 | 169 | -13.1
cc-pvdz-jkfit | 1392 | 68 | -36.3 | 124 | -25.6 | 193 | -7.63
ED threshold $10^{-5}$ |  |  |  |
$s$ | $s$ | 36 | 35 | 110.8 | 36 | 2848. | 36 | 18380.
$sp$ | $s$ | 72 | 49 | 2.151 | 68 | 78.31 | 72 | 930.0
$sp$ | $sp$ | 144 | 62 | -0.22 | 99 | -22.9 | 128 | 83.44
$spd$ | $s$ | 132 | 56 | -4.71 | 96 | -8.69 | 123 | 3.119
$spd$ | $spd$ | 324 | 71 | -0.66 | 131 | -5.98 | 196 | -20.9
weigend | 852 | 86 | -4.39 | 146 | -2.09 | 224 | -2.48
cc-pvdz-jkfit | 1392 | 91 | -2.67 | 164 | -1.79 | 253 | -4.43
ED threshold $10^{-6}$ |  |  |  |
$s$ | $s$ | 36 | 36 | 117.1 | 36 | 2848. | 36 | 18380.
$sp$ | $s$ | 72 | 59 | 1.313 | 72 | 79.44 | 72 | 930.0
$sp$ | $sp$ | 144 | 80 | -0.20 | 117 | 4.703 | 141 | 76.77
$spd$ | $s$ | 132 | 73 | -0.21 | 115 | -2.67 | 131 | 2.951
$spd$ | $spd$ | 324 | 93 | -0.25 | 167 | -1.43 | 235 | -4.75
weigend | 852 | 108 | -0.49 | 187 | -0.26 | 291 | -0.16
cc-pvdz-jkfit | 1392 | 117 | -0.19 | 209 | -0.22 | 322 | -0.65
Table 2: The overall errors in HF energy of (H2O)12 with PW basis. The reference HF energy is -912.40220404. PW Basis | $\omega=0.1$ | $\omega=0.15$ | $\omega=0.2$
---|---|---|---
$N_{r}$ | $N_{\Theta}$ | $N_{G}$ | $\Delta E/mE_{h}$ | $N_{G}$ | $\Delta E/mE_{h}$ | $N_{G}$ | $\Delta E/mE_{h}$
15 | 100 | 944 | -10.190 | 960 | -30.050 | 964 | 3.496
15 | 200 | 1852 | 3.323 | 1888 | 58.706 | 1896 | 278.077
15 | 300 | 2738 | -0.546 | 2786 | -16.148 | 2806 | -103.990
15 | 400 | 3772 | 0.057 | 3836 | 2.994 | 3884 | 26.264
20 | 100 | 1096 | -10.190 | 1100 | -30.050 | 1116 | 3.496
20 | 200 | 2144 | 3.323 | 2152 | 58.706 | 2184 | 278.077
20 | 300 | 3159 | -0.546 | 3176 | -16.148 | 3236 | -103.990
20 | 400 | 4372 | 0.057 | 4388 | 2.994 | 4464 | 26.264
30 | 400 | 5488 | 0.058 | 5508 | 2.994 | 5544 | 26.264
Table 3: The overall errors in HF energy of Gly6 with auxiliary Gaussian basis. The reference HF energy is -1316.96782986. Auxiliary Basis | $N_{aux}$ | $\omega=0.1$ | $\omega=0.15$ | $\omega=0.2$
---|---|---|---|---
C,N,O | H |  | $N_{ED}$ | $\Delta E/\mu E_{h}$ | $N_{ED}$ | $\Delta E/\mu E_{h}$ | $N_{ED}$ | $\Delta E/\mu E_{h}$
ED threshold $10^{-4}$ |  |  |  |
$s$ | $s$ | 45 | 33 | 867.9 | 44 | 7779. | 45 | 35595.
$sp$ | $s$ | 120 | 46 | 20.34 | 71 | 175.5 | 97 | 1224.
$sp$ | $sp$ | 180 | 55 | -3.81 | 96 | 74.34 | 131 | 286.1
$spd$ | $s$ | 245 | 50 | 14.04 | 87 | -23.5 | 131 | -13.9
$spd$ | $spd$ | 405 | 59 | -8.64 | 115 | -27.4 | 184 | -56.5
ED threshold $10^{-5}$ |  |  |  |
$s$ | $s$ | 45 | 42 | 751.6 | 45 | 7842. | 45 | 35595.
$sp$ | $s$ | 120 | 60 | 9.802 | 91 | 152.5 | 106 | 1120.
$sp$ | $sp$ | 180 | 73 | -4.04 | 122 | 11.82 | 152 | 149.9
$spd$ | $s$ | 245 | 68 | 0.747 | 119 | -6.92 | 165 | -51.4
$spd$ | $spd$ | 405 | 85 | 0.052 | 156 | -1.42 | 233 | -12.4
ED threshold $10^{-6}$ |  |  |  |
$s$ | $s$ | 45 | 44 | 744.0 | 45 | 7842. | 45 | 35595.
$sp$ | $s$ | 120 | 70 | 7.205 | 99 | 147.3 | 117 | 1103.
$sp$ | $sp$ | 180 | 94 | 0.290 | 139 | 10.05 | 170 | 219.8
$spd$ | $s$ | 245 | 90 | -0.33 | 148 | -6.35 | 195 | -49.2
$spd$ | $spd$ | 405 | 111 | -0.23 | 201 | -0.79 | 280 | -5.20
Table 4: Convergence tests for HF with LRDF. Convergence threshold for HF energy is set to $10^{-7}$ | $N_{AO}$ | Auxiliary Basis | $\omega$ | ED thresh | $N_{ED}$ | HF energy/Eh | HF iters. | CPU time/s
---|---|---|---|---|---|---|---|---
|  | C,N,O | H |  |  |  |  |  |
(H2O)128 | 3072 | conventional analytical ERIs | -9732.741296 | 11 | 64357
(H2O)128 | 3072 | $sp$ | $sp$ | 0.1 | $10^{-4}$ | 265 | -9732.741301 | 11 | 34156
(H2O)128 | 3072 | $spd$ | $s$ | 0.15 | $10^{-4}$ | 517 | -9732.742271 | 11 | 29299
gly30 | 2154 | conventional analytical ERIs | -6280.889857 | 13 | 19398
gly30 | 2154 | $sp$ | $sp$ | 0.1 | $10^{-4}$ | 235 | -6280.889816 | 13 | 7479
gly30 | 2154 | $spd$ | $s$ | 0.15 | $10^{-4}$ | 382 | -6280.890081 | 13 | 6864

###  III.1 Auxiliary basis

To evaluate the impact of auxiliary basis, two groups of auxiliary basis are
considered in these tests:

  * •

The regular density fitting auxiliary basis (cc-pvdz-jkfitWeigend (2002,
2006), Weigend Coulomb-fit basisWeigend  _et al._ (1998));

  * •

A set of light-weight auxiliary basis. Each light-weight auxiliary basis only
consists of primitive Gaussian orbitals with distinct angular momentum. The
exponents of the Gaussian auxiliary basis are uniformly set to 1.0.

In the DF approximation, the main purpose of the auxiliary basis is to
describe the multipoles of the orbital products. To capture the effects of
charge (monopole) on each atom, at least one $s$-type Gaussian function is
required on each atom, which constitutes the smallest basis set used in the
tests. This smallest basis set represents the worst-case scenario for the LRDF
method. Nevertheless, the error is less than 10 $\mu E_{h}$ per atom in the
low $\omega$ region (such as $\omega=0.1$). To achieve a more accurate
description of the high order multipoles, polarized functions such as $p$-type
and $d$-type Gaussian functions can be added to the auxiliary basis. The
inclusion of these functions leads to a significant improvement in accuracy
Even in the worst-case scenario where $\omega=0.2$ and the ED threshold is
$10^{-4}$. The $spd$ basis, which consists of one $s$, one $p$, and one
$d$-type Gaussian on each atom, reduces the error to $0.1\mu E_{h}$ per atom
for the water cluster and $1\mu E_{h}$ per atom for the linear Gly6 molecule.
When $\omega=0.1$ and ED threshold $10^{-5}$ are utilized, the errors in both
molecules are approximately $0.01\mu E_{h}$ per atom, which is negligible.

To compare the present LRDF method with the traditional DF method, we also
performed HF energy calculations using the traditional DF method. The errors
of the traditional DF method are significantly larger than the errors of the
LRDF method. For the water cluster, the traditional DF method yields an error
of 10003 $\mu E_{h}$ with the Weigend Coulomb-fit basis and 348 $\mu E_{h}$
with the cc-pvdz-jkfit basis. These errors are near the errors of the LRDF
method with the smallest auxiliary basis.

Let us now consider the size of the auxiliary basis. The size of the auxiliary
basis is denoted as $N_{aux}$, while the number of linearly independent ED
vectors is denoted as $N_{ED}$. It is worth to note that the regular DF basis
is actually over-complete when applied to the LR Coulomb potential metric.
When we employ the ED method to eliminate linear dependency, we find that 70%
to 95% of the auxiliary basis can be discarded. Even in the case of
lightweight auxiliary bases, the ED process still reduces the problem size by
70%. This suggests that the lightweight auxiliary basis may still be over-
complete, and there are room for further improvement.

###  III.2 Range separation parameter $\omega$

The range separation parameter, $\omega$ controls the exposure of the local
structure of the electron density in the LR DF integrals. By choosing a small
value of $\omega$, fewer short-range details are exposed to the LR DF part,
leading to more accurate results and reduced requirements on the auxiliary
basis. Moreover, a small $\omega$ can effectively reduce the size of linearly
independent ED vectors. As shown in equation (26), the computational cost of
DF HFX is directly proportional to the size of linearly independent ED
vectors. Therefore, a relatively small $\omega$ is favorable due to the
$N^{4}$ scaling in the computation of DF HFX. However, a smaller $\omega$
results in a slower decay of the SR potential and increased computational
efforts for the analytical SR ERIs. It is not clear in current study what is
the optimal value of $\omega$ that balances accuracy and computational costs.

###  III.3 Threshold of eigenvalue decomposition

As shown by the $N_{ED}$ columns in the results tables, the ED threshold has
an impact on the size of linearly independent ED vectors, which in turn
affects the computational costs of DF HFX. In the case of dense systems such
as the water cluster, tightening the ED threshold can significantly improve
the accuracy. In the linear molecule, ED threshold is not the dominant impact
compared to the influence of the auxiliary basis. With a small value of
$\omega$ such as $\omega=0.1$, an ED threshold of $10^{-4}$ is generally
sufficient to achieve sub-micro Hartree accuracy.

###  III.4 Integration with PW basis

In Table 2, the number of radial grids and the number of angular grids used
for the reciprocal space integral (9) are represented by the symbols $N_{r}$
and $N_{\Theta}$ respectively. The two parameters are combined to generate all
potential plane wave (PW) basis functions, which are then trimmed based on the
resulting integration weights to reduce computational expenses. The number of
remaining PW functions is denoted as $N_{G}$ in the table.

The accuracy of LR ERIs in this approach relies heavily on the density of
integration grids. As a rough estimation, when $\omega=0.1$, selecting 2
radial grids per Bohr relative to the largest atomic distance in the molecule,
and approximately $N_{\Theta}\approx N_{r}^{2}$ angular grids, can yield an
accuracy of approximately $1\mu E_{h}$ per atom. For instance, the water
cluster requires $18\times 300$ grids, while the linear Gly6 molecule requires
$70\times 5000$ grids. After eliminating the PW basis functions that have
negligible integration weights, 3000 PWs are still required for the water
molecule, while the linear molecule requires 100k PWs. This is 100 - 1000
times more computationally demanding than the case with Gaussian auxiliary
basis.

###  III.5 SCF performance

The SCF convergence achieved with the LRDF method closely matches that of the
SCF with conventional analytical ERIs. There is no accumulation of errors in
larger systems during the SCF iterations. When $\omega=0.1$, the errors in the
converged HF energy are approximately 5 $\mu E_{h}$ for the 128 water cluster
and 40 $\mu E_{h}$ for the linear molecule Gly30. The average error is only
around $0.1\mu E_{h}$ per atom.

The LRDF method significantly reduce computational time. The water cluster
poses a challenge for screening SR integrals since many atoms are close to
each other. Despite, we observe a reduction in computational time of 50-60% in
the water cluster. The atoms in the Gly30 molecule are more widely
distributed. The integral screening is more effective in this system. We
obtain a 60-70% reduction in computational time. It can be expected that the
computational cost should be further reduced in larger systems.

The performance of the SR ERI plays a crucial role in this algorithm. In the
Libcint library, ERIs are evaluated with the Rys-quadrature algorithm. For low
angular momentum GTO functions, the quadrature roots and weights of Rys
polynomials are fitted using Chebyshev polynomials when computing the full-
range Coulomb integrals. This polynomial fitting technique is particularly
efficient for the full-range ERIs. However, there are some challenges to
develop efficient polynomial fitting schemes for Rys-quadrature in SR ERIs. In
previous versions of Libcint (v5.*), the quadratures were computed from
scratch at runtime, which significantly slows down the performance of SR ERIs.
In Libcint v6.0, SR ERIs are computed as the difference between full-range
ERIs and long-range ERIs for low angular momentum GTO functions. Although this
treatment requires the evaluation of twice the number of integrals, it is more
efficient than calculating the Rys-quadrature of SR ERIs from scratch. This
implies that there is room for improving the algorithm of SR ERIs, which could
further accelerate the SCF performance.

##  IV Conclusions

In this work, we explore the idea of range-separation with long-range density
fitting to compute Coulomb and exchange terms of Hartree-Fock method for
molecules. By combining the analytical algorithm for short-range ERIs and
long-range density fitting ERIs, we find an efficient approach to approximate
ERIs for large molecule systems. Compared to Hartree-Fock with conventional
analytical ERIs, the new method achieves more than two-fold acceleration in
performance without compromising accuracy. The error is around 0.1 $\mu E_{h}$
per atom. Our preliminary study shows that the long-range density fitting
component has a lower requirement on the auxiliary basis compared to the
regular density fitting method. Both Gaussian basis and plane-wave basis can
be used as the auxiliary basis. We can achieve errors of 0.1 $\mu E_{h}$ per
atom with a very small auxiliary Gaussian basis, which is only 1/10 the size
of orbital basis. However, to obtain a similar level of accuracy, it is
necessary to use a larger number of plane-wave basis functions, approximately
100 times more than the auxiliary Gaussian basis.

## References

  * Whitten (1973) J. L. Whitten, “Coulombic potential energy integrals and approximations,” J. Chem. Phys. 58, 4496–4501 (1973).
  * Dunlap, Connolly, and Sabin (1979) B. I. Dunlap, J. W. D. Connolly, and J. R. Sabin, “On some approximations in applications of xa theory,” J. Chem. Phys. 71, 3396–3402 (1979).
  * Vahtras, Almlöf, and Feyereisen (1993) O. Vahtras, J. Almlöf, and M. Feyereisen, “Integral approximations for lcao-scf calculations,” Chem. Phys. Lett. 213, 514 – 518 (1993).
  * Eichkorn  _et al._ (1995) K. Eichkorn, O. Treutler, H. Öhm, M. Häser, and R. Ahlrichs, “Auxiliary basis sets to approximate coulomb potentials,” Chem. Phys. Lett. 240, 283 – 290 (1995).
  * Dunlap (2000) B. I. Dunlap, “Robust and variational fitting,” Phys. Chem. Chem. Phys. 2, 2113–2116 (2000).
  * te Velde and Baerends (1991) G. te Velde and E. J. Baerends, “Precise density-functional method for periodic structures,” Phys. Rev. B 44, 7888–7903 (1991).
  * Jaffe and Hess (1996) J. E. Jaffe and A. C. Hess, “Gaussian basis density functional theory for systems periodic in two or three dimensions: Energy and forces,” J. Chem. Phys. 105, 10983–10998 (1996).
  * Čársky, Čurík, and Varga (2012) P. Čársky, R. Čurík, and v. Varga, “Efficient evaluation of coulomb integrals in a mixed gaussian and plane-wave basis using the density fitting and cholesky decomposition,” J. Chem. Phys. 136, 114105 (2012).
  * Burow, Sierka, and Mohamed (2009) A. M. Burow, M. Sierka, and F. Mohamed, “Resolution of identity approximation for the coulomb term in molecular and periodic systems,” J. Chem. Phys. 131, 214101 (2009).
  * Feyereisen, Fitzgerald, and Komornicki (1993) M. Feyereisen, G. Fitzgerald, and A. Komornicki, “Use of approximate integrals in ab initio theory. an application in {MP2} energy calculations,” Chem. Phys. Lett. 208, 359 – 363 (1993).
  * Komornicki and Fitzgerald (1993) A. Komornicki and G. Fitzgerald, “Molecular gradients and hessians implemented in density functional theory,” J. Chem. Phys. 98, 1398–1421 (1993).
  * Manby, Knowles, and Lloyd (2001) F. R. Manby, P. J. Knowles, and A. W. Lloyd, “The poisson equation in density fitting for the kohn-sham coulomb problem,” J. Chem. Phys. 115, 9144–9148 (2001).
  * Jung  _et al._ (2005) Y. Jung, A. Sodt, P. M. W. Gill, and M. Head-Gordon, “Auxiliary basis expansions for large-scale electronic structure calculations,” Proc. Natl. Acad. Sci. U. S. A. 102, 6692–6697 (2005), http://www.pnas.org/content/102/19/6692.full.pdf .
  * Hollman, Schaefer, and Valeev (2014) D. S. Hollman, H. F. Schaefer, and E. F. Valeev, “Semi-exact concentric atomic density fitting: Reduced cost and increased accuracy compared to standard density fitting,” J. Chem. Phys. 140, 064109 (2014).
  * Sun  _et al._ (2017) Q. Sun, T. C. Berkelbach, J. D. McClain, and G. K.-L. Chan, “Gaussian and plane-wave mixed density fitting for periodic systems,” J. Chem. Phys. 147, 164119 (2017).
  * Ye and Berkelbach (2021) H.-Z. Ye and T. C. Berkelbach, “Fast periodic gaussian density fitting by range separation,” J. Chem. Phys. 154, 131104 (2021).
  * Lee  _et al._ (2022) J. Lee, A. Rettig, X. Feng, E. Epifanovsky, and M. Head-Gordon, “Faster exact exchange for solids via occ-ri-k: Application to combinatorially optimized range-separated hybrid functionals for simple solids with pseudopotentials near the basis set limit,” J. Chem. Theory Comput. 18, 7336–7349 (2022).
  * Pritchard  _et al._ (2019) B. P. Pritchard, D. Altarawy, B. Didier, T. D. Gibson, and T. L. Windus, “New basis set exchange: An open, up-to-date resource for the molecular sciences community,” J. Chem. Inf. Model. 59, 4814–4820 (2019).
  * Yang, Rendell, and Frisch (2007) R. Yang, A. P. Rendell, and M. J. Frisch, “Automatically generated coulomb fitting basis sets: Design and accuracy for systems containing h to kr,” J. Chem. Phys. 127, 074102 (2007).
  * Boström  _et al._ (2009) J. Boström, F. Aquilante, T. B. Pedersen, and R. Lindh, “Ab initio density fitting: Accuracy assessment of auxiliary basis sets from cholesky decompositions,” J. Chem. Theory Comput. 5, 1545–1553 (2009).
  * Aquilante  _et al._ (2009) F. Aquilante, L. Gagliardi, T. B. Pedersen, and R. Lindh, “Atomic cholesky decompositions: A route to unbiased auxiliary basis sets for density fitting approximation with tunable accuracy and efficiency,” J. Chem. Phys. 130, 154107 (2009).
  * Kállay (2014) M. Kállay, “A systematic way for the cost reduction of density fitting methods,” J. Chem. Phys. 141, 244113 (2014).
  * Stoychev, Auer, and Neese (2017) G. L. Stoychev, A. A. Auer, and F. Neese, “Automatic generation of auxiliary basis sets,” J. Chem. Theory Comput. 13, 554–562 (2017).
  * Lehtola (2021) S. Lehtola, “Straightforward and accurate automatic auxiliary basis set generation for molecular calculations with atomic orbital basis sets,” J. Chem. Theory Comput. 17, 6886–6900 (2021).
  * Lehtola (2023) S. Lehtola, “Automatic generation of accurate and cost-efficient auxiliary basis sets,” J. Chem. Theory Comput. (2023), 10.1021/acs.jctc.3c00670.
  * De Chavez and Lindh (2023) D. De Chavez and R. Lindh, “One-centre corrected two-electron integrals in inner projection-based integral evaluations*,” Mol. Phys. , e2228431 (2023).
  * Sun (2023a) Q. Sun, “Exact exchange with range-separated algorithm for thermodynamic limit of periodic hartree-fock theory,” J. Chem. Phys. 159, 024108 (2023a).
  * Lehtola (2019) S. Lehtola, “Curing basis set overcompleteness with pivoted cholesky decompositions,” J. Chem. Phys. 151, 241102 (2019).
  * Besalú and Carbó-Dorca (2011) E. Besalú and R. Carbó-Dorca, “The general gaussian product theorem,” J. Math. Chem. 49, 1769–1784 (2011).
  * Boström  _et al._ (2013) J. Boström, F. Aquilante, T. B. Pedersen, and R. Lindh, “Analytical gradients of hartree-fock exchange with density fitting approximations,” J. Chem. Theory Comput. 9, 204–212 (2013).
  * Almlöf, Faegri Jr., and Korsell (1982) J. Almlöf, K. Faegri Jr., and K. Korsell, “Principles for a direct scf approach to licao–moab-initio calculations,” J. Comput. Chem. 3, 385–399 (1982), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.540030314 .
  * Izmaylov, Scuseria, and Frisch (2006) A. F. Izmaylov, G. E. Scuseria, and M. J. Frisch, “Efficient evaluation of short-range hartree-fock exchange in large molecules and periodic systems,” J. Chem. Phys. 125, 104103 (2006).
  * Thompson and Ochsenfeld (2017) T. H. Thompson and C. Ochsenfeld, “Distance-including rigorous upper bounds and tight estimates for two-electron integrals over long- and short-range operators,” J. Chem. Phys. 147, 144101 (2017), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4994190/14837196/144101_1_online.pdf .
  * Sun (2023b) Q. Sun, “Various integral estimations and screening schemes for extended systems in pyscf,” (2023b), arXiv:2302.11307 [physics.chem-ph] .
  * Manzer, Epifanovsky, and Head-Gordon (2015) S. F. Manzer, E. Epifanovsky, and M. Head-Gordon, “Efficient implementation of the pair atomic resolution of the identity approximation for exact exchange for hybrid and range-separated density functionals,” J. Chem. Theory Comput. 11, 518–527 (2015), http://dx.doi.org/10.1021/ct5008586 .
  * Sun  _et al._ (2020) Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, “Recent developments in the pyscf program package,” J. Chem. Phys. 153, 024109 (2020), https://doi.org/10.1063/5.0006074 .
  * Sun (2015) Q. Sun, “Libcint: An efficient general integral library for gaussian basis functions,” J. Comput. Chem. 36, 1664–1671 (2015).
  * Weigend (2002) F. Weigend, “A fully direct ri-hf algorithm: Implementation, optimised auxiliary basis sets, demonstration of accuracy and efficiency,” Phys. Chem. Chem. Phys. 4, 4285–4291 (2002).
  * Weigend (2006) F. Weigend, “Accurate coulomb-fitting basis sets for h to rn,” Phys. Chem. Chem. Phys. 8, 1057–1065 (2006).
  * Weigend  _et al._ (1998) F. Weigend, M. Häser, H. Patzelt, and R. Ahlrichs, “Ri-mp2: optimized auxiliary basis sets and demonstration of efficiency,” Chem. Phys. Lett. 294, 143 – 152 (1998).