#  Effective Many-body Interactions in Reduced-Dimensionality Spaces Through
Neural Network Models

Senwei Liang<EMAIL_ADDRESS>Applied Mathematics and Computational
Research Division
Lawrence Berkeley National Laboratory  Karol Kowalski<EMAIL_ADDRESS>Physical Sciences Division, Pacific Northwest National Laboratory, Richland,
Washington 99354, USA  Chao Yang<EMAIL_ADDRESS>Applied Mathematics and
Computational Research Division
Lawrence Berkeley National Laboratory  Nicholas P. Bauman
<EMAIL_ADDRESS>Physical Sciences Division, Pacific Northwest
National Laboratory, Richland, Washington 99354, USA

###### Abstract

Accurately describing properties of challenging problems in physical sciences
often requires complex mathematical models that are unmanageable to tackle
head-on. Therefore, developing reduced-dimensionality representations that
encapsulate complex correlation effects in many-body systems is crucial to
advance the understanding of these complicated problems. However, a numerical
evaluation of these predictive models can still be associated with a
significant computational overhead. To address this challenge, in this paper,
we discuss a combined framework that integrates recent advances in the
development of active-space representations of coupled cluster (CC) downfolded
Hamiltonians with neural network approaches. The primary objective of this
effort is to train neural networks to eliminate the computationally expensive
steps required for evaluating hundreds or thousands of Hugenholtz diagrams,
which correspond to multidimensional tensor contractions necessary for
evaluating a many-body form of downfolded/effective Hamiltonians. Using small
molecular systems (the H2O and HF molecules) as examples, we demonstrate that
training neural networks employing effective Hamiltonians for a few nuclear
geometries of molecules can accurately interpolate/extrapolate their forms to
other geometrical configurations characterized by different intensities of
correlation effects. We also discuss differences between effective
interactions that define CC downfolded Hamiltonians with those of bare
Hamiltonians defined by Coulomb interactions in the active spaces.

††preprint: APS/123-QED

##  I Introduction

The development of methodologies and approximations that effectively capture
complex correlation effects is crucial in computational chemistry and physics
for comprehending the diverse properties of quantum systems. Various
representations of quantum mechanics, such as wave function,[1, 2, 3] density
matrix, Green’s function, and electron density-based formulations have been
employed over the last few decades to enable these approaches. For example,
the Density Functional Theory (DFT) has evolved into the driving engine for
numerous electronic structure simulations in chemistry and materials sciences,
aiming to approximate exchange-correlation potential. On the other end of the
spectrum, a hierarchical family of approximations provided by a broad class of
coupled cluster (CC) formalisms offers a highly accurate characterization of
ground- and excited-state processes, converging in the asymptotic limit to the
exact theory. In many instances, formulations used in the electronic structure
theory provide a universal computational framework that works across length,
energy, and time scales, as demonstrated by the universality of the CC
formalism used in chemistry, materials [4, 5, 6, 7, 8], and nuclear structure
theory simulations [9, 10, 11, 12]. However, the main challenge that precludes
the widespread utilization of these methods is the polynomial scaling with the
system size, which significantly narrows the corresponding area of
applications, particularly in the context of high-accuracy methodologies.

To address this challenge, several approaches to handle numerical overhead
have been pursued, including reducing numerical costs of methodologies and
embracing high-performance conventional computing, as well as attempting to
utilize quickly evolving quantum computing technologies. While high-
performance computing (HPC) is currently a necessary component of nearly all
computational packages used in simulations, the best effect is achieved when
HPC technologies are integrated with the theoretical frameworks to either
reduce the scaling or the dimensionality of the problem. Reduced-scaling
formulations, commonly referred to as local formulations, enforce sparsity
conditions on the dense representations of equations [13, 14, 15, 16]. On the
other hand, reduced-dimensionality approaches involve representing the problem
in reduced-dimensionality subspace(s), also known as active space(s). In
general, reduced-dimensionality representations provide an alternative means
of incorporating sparsity compared to typical reduced-scaling approaches.

In recent years, significant inroads have been made in developing downfolded
Hamiltonians for quantum chemistry and quantum computing applications. These
Hamiltonians are designed to accurately reproduce the system’s ground-state
energy in active spaces by integrating external correlation effects in the
many-body form of the Hamiltonians. The CC downfolding approach results from
the exponential parametrization (ansatz) for the ground-state wave function in
the single-reference parametrization. It can be considered a many-body
embedding or renormalization procedure when applying multi-step downfolding
procedures. These formulations provide a rigorous theoretical and
computational framework for calculating effective Hamiltonians, promising
prospects for further advancements in quantum chemistry, quantum embedding
theory, and quantum computing.[17, 18, 19, 20]

When it comes to utilizing CC downfolding methods in practical applications,
the most significant challenge is the numerical effort required to evaluate
many tensor contractions associated with the Hugenholtz diagrams that define
the many-body structure of an effective Hamiltonian. Even HPC technologies may
only provide a viable solution within a particular system size limit, and
therefore a computational/evaluation algorithm change is needed for larger and
more complex applications. One potential solution to these issues is using
neural networks (NN). This choice of computational methodology is motivated by
recent applications of deep NNs (DNN) to solve many-body Schrödinger equation
and identify corresponding wave function [21, 22, 23, 24] (see also Ref.[25]
for the related effort to solve integro-differential equations). Our approach
is driven by a different paradigm closely related to the ability of downfolded
Hamiltonians to provide highly accurate estimates of the total ground-state
energy. Instead of attempting to determine wave functions, which generally
require evaluating a massive number of parameters and properly accounting for
the anti-symmetry for Fermionic systems, we use the DNN to evaluate a limited
set of interactions contained within small-dimensionality active spaces. Even
though downfolded Hamiltonians generally contain interactions beyond pairwise
interactions, the number of these interactions (using tensors and second-
quantized representations to represent them) is significantly smaller than the
number of wave function parameters for general dense systems typically
encountered in quantum chemistry.

The proposed framework can be regarded as an instance of physics-informed
neural networks (PINNs) [26], which leverage the fundamental principles of the
lower-dimensional representations of quantum systems to design machine
learning (ML) models that take into account special structures of interactions
in diverse physical scenarios.

In this paper, for the first time, we put forward a DNN framework (referred to
as the VNet) to evaluate effective interactions in correlated systems. To
assess the feasibility of the VNet, we focus entirely on molecular benchmark
systems (the H2O and HF molecules) and two-body effective interactions. The
VNet DNN takes four indices of the tensor representing two-body interactions
of the effective Hamiltonians in the active space, as well as the geometry of
the molecule, as the input, and predicts the value of the corresponding tensor
element as the output. We assume that $P$, $Q$, $R$, and $S$ spin-orbital
indices designate active spin-orbitals. The output is used in a loss function
that measures the discrepancy between the predicted and the true value of the
tensor elements $k^{PQ}_{RS}$ contained in a training dataset that consists of
easily obtainable elements of the active-space effective interactions
associated with the downfolded Hamiltonians for a few selected geometries.
Also included in the training dataset are two-body bare interactions
$v^{PQ}_{RS}$ associated with the original Hamiltonian in the same active
space. The DNN parameters in VNet are optimized in a batched training
procedure that minimizes the loss function in an iterative scheme. The
training is carried out in two steps. In the first step, several bare
interaction tensors contained in the training dataset are used to initialize
the VNet parameters. These parameters are finetuned in a subsequent step in
which selected effective interaction tensors associated with a few selected
geometries are used to finetune the VNet parameters. An exciting feature of
VNet is that the DNN optimization cost (excluding the formation cost of
training data - downfolded Hamiltonians) depends only on the number of active
orbitals.

The paper is organized as follows. In section II, we briefly describe the
coupled cluster downfolding theory that can be used to derive an effective
downfolded Hamiltonian within an active space. The second quantized
representation of the downfolded effective Hamiltonian is given in section
III. The main ideas of using machine learning ML techniques to predict
effective interactions for different molecular geometries using a carefully
designed and physics-informed neural network trained with bare and effective
interactions at a few selected molecular geometries are presented in section
IV. We demonstrate the effectiveness of the ML approach in section V using the
H2O and HF molecules in different geometric configurations as test cases. We
also show that the ML approach not only allows us to predict effective
interactions but also provides information that can possibly be used to
characterize the nature of screening effects in the effective interactions.
Some main conclusions and future directions are discussed in section VI.

##  II Coupled cluster downfolding

In contrast to the numerous effective Hamiltonian formulations discussed in
physics and chemistry, the CC downfolding discussed in Refs.[17, 27] is a
consequence of the ground-state wave function’s single-reference exponential
representation (ansatz). The existing two classes of CC downfolding, Hermitian
and non-Hermitian CC downfolding procedures, utilize the standard single-
reference CC (SR-CC) and double unitary CC (DUCC) representations of the
ground-state wave functions where the ansatz can be naturally factorized into
products of exponential ansatzen defined by internal and external cluster
operators generating correlation effects inside/outside the active spaces,
respectively. Below, we provide a brief review of these formulations.

###  II.1 Non-Hermitian CC Downfolding

The non-Hermitian CC downfolding is a consequence of the SR-CC exponential
ansatz, where the ground-state wave function $|\Psi\rangle$ in exact or
approximate representations is expressed as

$|\Psi\rangle=e^{T}|\Phi\rangle\;,$ (1)

where $T$ is the cluster operator and $|\Phi\rangle$ designates the so-called
reference function, usually chosen as a Hartree–Fock (HF) determinant. We also
introduce the notion of the complete active space (CAS), which is generated by
the sets of active occupied ($R$) and virtual ($S$) orbitals, where all
excited Slater determinants spanning CAS can be generated as the action of
elements of the commutative sub-algebra
$\mathfrak{h}\equiv\mathfrak{g}^{(N)}(R,S)$ (expressed in the particle-hole
representation with respect to the reference Slater determinant
$|\Phi\rangle$) on the reference function $|\Phi\rangle$. The sub-algebra
$\mathfrak{h}$ in a natural induces the partitioning of the cluster operator
$T$ into internal, $T_{\rm int}(\mathfrak{h})$, producing excited
configurations within CAS when acting on $|\Phi\rangle$, and external, $T_{\rm
ext}(\mathfrak{h})$, producing excited configurations in the orthogonal
complement of CAS when acting on $|\Phi\rangle$, parts, i.e.,

$T=T_{\rm int}(\mathfrak{h})+T_{\rm ext}(\mathfrak{h})$ (2)

and

$|\Psi\rangle=e^{T_{\rm ext}(\mathfrak{h})}e^{T_{\rm
int}(\mathfrak{h})}|\Phi\rangle\;.$ (3)

If the $T_{\rm int}(\mathfrak{h})$ operator generates all possible excited
determinants when acting on $|\Phi\rangle$ we call the $\mathfrak{h}$ sub-
algebra a sub-system embedding sub-algebra (SES) for a given representation of
the cluster operator $T$.

The partitioning of the cluster operator into internal and external parts has
been initially introduced in the context of the state-selective CC formalism
in Refs.[28, 29, 30, 31] and was invoked as a selection mechanism for high-
rank cluster amplitudes. An analogous representation of the cluster operator
was used to develop various tailored CC methods [32, 33].

Employing the idea of SESs one can demonstrate that the CC can be obtained in
an alternative way to the standard text-book formula

$E=\langle\Phi|e^{-T}He^{T}|\Phi\rangle,$ (4)

namely as an eigenvalue of effective or downfolded Hamiltonian $H^{\rm
D}(\mathfrak{h})$ in the $\mathfrak{h}$-generated CAS

$H^{\rm D}(\mathfrak{h})e^{T_{\rm int}(\mathfrak{h})}|\Phi\rangle=Ee^{T_{\rm
int}(\mathfrak{h})}|\Phi\rangle\;,$ (5)

where

$H^{\rm D}(\mathfrak{h})=(P+Q_{\rm int}(\mathfrak{h}))e^{-T_{\rm
ext}(\mathfrak{h})}He^{T_{\rm ext}(\mathfrak{h})}(P+Q_{\rm
int}(\mathfrak{h}))\;,$ (6)

and $P$ and $Q_{\rm int}(\mathfrak{h})$ refer to the projection operator onto
the reference function and all excited Slater determinants belonging to
$\mathfrak{h}$-generated CAS. The standard energy expression (4) is a special
case of (6) for $R$ and $S$ being empty sets (in such a case the CAS is
generated by a reference function $|\Phi\rangle$).

Since in the effective Hamiltonian (6) the external degrees of freedom are
integrated out, the CC formalism can be viewed as a special example of an
embedding procedure (see Refs.[34]) or renormalization algorithm. It is also
instructive to scrutinize the many-body structure of the $H^{\rm
D}(\mathfrak{h})$ operator and ensuing properties:

  * •

Due to the commutative character of all components defining $T_{\rm
ext}(\mathfrak{h})$ the exact form of the effective Hamiltonian can be
constructed. As in the standard SR-CC formulations, the effective Hamiltonian
will involve terms that are, at most, fourth powers of cluster operators.

  * •

The effective Hamiltonians can be constructed for each SES and the
corresponding active space. In Ref.[34] the number of SESs defined at the
orbital level for the CCSD formalism (CC with singles and doubles of Ref.[35])
amounts to

$n_{o}(2^{n_{v}}-1)+n_{v}(2^{n_{o}}-1)-n_{o}n_{v}\;,$ (7)

where $n_{o}$ and $n_{v}$ refer to the number of occupied and unoccupied
orbitals, respectively.

  * •

The SES active-space problems can be integrated into the flows, where total
cluster operator $T_{\rm F}$ is defined through the unique excitation defining
internal excitations for active spaces involved in the flow. The so-called
Equivalence Theorem, see Ref.[36], states that the solution of the flow
equations defined in such a way is equivalent to the solution of standard SR-
CC equations defined by the operator $T_{\rm F}$ and therefor can be viewed as
a ”correlated” variant of Aufbau principle (in the sense of adding a sub-
system described by its active space). The SR-CC flow equations have been
tested for ground- and excited-state simulations [17].

  * •

The non-Hermitian CC downfolding and flow equations can be extended to the
time-domain [37, 36]. In Ref.[37] it has been shown that when the fast-varying
part of the wave function defined by the $T_{\rm ext}(\mathfrak{h},t)$ is
known then the slow-varying dynamics of the entire system can be described by
downfolded Hamiltonian $H^{\rm D}(\mathfrak{h},t)$.

  * •

The effective Hamiltonian and its many-body structure provide a theoretical
framework to describe a sub-system defined by active space as an open system
interacting with the surrounding environment. In Fig. 1, we provided a
specific example demonstrating that although the effective Hamiltonian
preserves the number of active electrons, the intermediate state (the $S_{2}$
state) can involve the exchange of electrons between the sub-system and
surrounding environment.

Figure 1:  Schematic representation of many-body processes encapsulated in the
many-body structure of $H^{\rm D}(\mathfrak{h})$. The leftmost panel
represents a one-body Hugenholtz diagram contributing to the effective
Hamiltonians $H^{\rm D}(\mathfrak{h})$. The black circle represents the two-
body interactions of the bare Hamiltonians, the white oval represents the
doubly excited external cluster amplitude, whereas the black rectangle
corresponds to the singly excited Slater determinant
$|\Phi_{I}^{A}\rangle=a^{\dagger}_{A}a_{I}|\Phi\rangle$
($a_{p}^{\dagger}$/$a_{p}$ represent creation/annihilation operators for
electron in $p$-th spin-orbital). The intermediate states $S_{1}$, $S_{2}$,
and $S_{3}$ correspond to the following Slater determinants:
$|\Phi_{I}^{A}\rangle$, $|\Phi_{IJ\bar{m}}^{A\bar{e}\bar{f}}\rangle$ and
$|\Phi_{J}^{A}\rangle$. The indices $I,J,...$ ($A,B,...$) correspond to the
active occupied (unoccupied) spin-orbital indices (occupied/unoccupied in
$|\Phi\rangle$), whereas indices $\bar{m},\bar{n}$ ($\bar{e},\bar{f},...$)
correspond to the inactive occupied (unoccupied) spin-orbitals. The three
panels to the right illustrate the electron excitation processes corresponding
to states $S_{1}$, $S_{2}$, and $S_{3}$. The green horizontal lines in these
panels separate occupied orbitals from unoccupied ones.

In ample numerical tests performed for various levels of CC theory and
published in Ref.[34], it was demonstrated that using Eq.(6), one can
reproduce exact values of CC calculated using the standard formula given by
Eq.(4) even for active spaces that would not provide a sufficient level of
correlation effect in diagonalizing bare electronic Hamiltonians.

The properties of the non-Hermitian CC downfolding offer not only a new way of
interpreting SR-CC theory but also provide guiding principles for the design
of its Hermitian counterpart.

###  II.2 Hermitian CC Downfolding

The Hermitian CC downfolding is based on the utilization of the DUCC ansatz,
[27] which assumes a similar form as Eq.(3), i.e.,

$|\Psi\rangle=e^{\sigma_{\rm ext}(\mathfrak{h})}e^{\sigma_{\rm
int}(\mathfrak{h})}|\Phi\rangle\;,$ (8)

where external and internal general-type anti-Hermitian cluster operators
$\sigma_{\rm ext}(\mathfrak{h})$ and $\sigma_{\rm int}(\mathfrak{h})$ are
expressed in terms of parameters carrying at least one inactive and no
inactive spin-orbital indices, respectively. In general, due to the non-
commutativity of many-body components defining $\sigma_{\rm
ext}(\mathfrak{h})$ and $\sigma_{\rm int}(\mathfrak{h})$ the analysis of the
DUCC-based Hermitian downfolding is more challenging. It requires several
approximate steps in practical applications. As in the non-Hermitian case, we
assume that the active space is defined by sub-algebra $\mathfrak{h}$.

In Ref.[37], we performed the analysis of the expansion (8) using disentangled
unitary coupled cluster formalism [38] in a way that it factorizes into two
parts defined products defined by cluster amplitudes (or parameters mentioned
earlier) carrying only active spin-orbital indices followed by part defined by
cluster amplitudes carrying at least one inactive spin-orbital index. Applying
the Baker–Campbell–Hausdorff formula multiple times to the first and the
second part leads Eq.(8) where $\sigma_{\rm ext}(\mathfrak{h})$ and
$\sigma_{\rm int}(\mathfrak{h})$ are defined by complex and non-terminating
commutator expansions. If these expansions exist, then the DUCC ansatz (8)
represents the exact ground-state wave function. It can be shown that if the
values of general-type anti-Hermitian cluster operator (in the sense of the
discussion of Ref.[39]) is known, then the exact energy can be reproduced by
diagonalizing Hermitian Hamiltonian $H^{\rm D}(\mathfrak{h})$ defined as:

$H^{\rm D}(\mathfrak{h})=(P+Q_{\rm int}(\mathfrak{h}))e^{-\sigma_{\rm
ext}(\mathfrak{h})}He^{\sigma_{\rm ext}(\mathfrak{h})}(P+Q_{\rm
int}(\mathfrak{h}))$ (9)

using any type of solver. It should be noted that, in analogy to non-Hermitian
CC downfolding, only external cluster operators are needed.

The Hermitian downfolding can also be extended to the flow formalisms (e.g.,
Quantum Flow (QFlow) formulation of Ref.[36, 20] introduced in the context of
quantum computing). As in the Hermitian case, the QFlow formalism can be
exploited to capture the sparsity of quantum systems.

##  III Numerical methods

When the $T_{\rm ext}(\mathfrak{h})$ or $\sigma_{\rm ext}(\mathfrak{h})$
operators are known or can be accurately approximated, particularly when
static correlation effects can be effectively incorporated by carefully
selecting the active space, both Hermitian and non-Hermitian effective
Hamiltonians (6) and (9) serve as frameworks for devising approximate methods
to calculate ground-state energies (see for example Refs.[17, 40, 18, 41]). In
this section, our primary emphasis will be on the Hermitian downfolding and
exploration of the usual approximations and computational resources required
to determine the approximate representations of the effective Hamiltonians.
However, many of these considerations can be easily extended to the Hermitian
case.

###  III.1 Approximations for Hermitian downfolded Hamiltonians

In many applications, the second-quantized representation of the $H^{\rm
D}(\mathfrak{h})$ operator

$H^{\rm D}(\mathfrak{h})=\Gamma^{\rm
D}_{0}(\mathfrak{h})+\sum_{PQ}g^{P}_{Q}(\mathfrak{h})a_{P}^{\dagger}a_{Q}+\frac{1}{4}\sum_{P,Q,R,S}k^{PQ}_{RS}(\mathfrak{h})a_{P}^{\dagger}a_{Q}^{\dagger}a_{S}a_{R}+\frac{1}{36}\sum_{P,Q,R,S,T,U}l^{PQR}_{STU}(\mathfrak{h})a_{P}^{\dagger}a_{Q}^{\dagger}a_{R}^{\dagger}a_{U}a_{T}a_{S}+\ldots\;,$
(10)

is required. In the above equations $P$, $Q$, $\ldots$ indices represent
general-type active spin-orbitals and $\Gamma^{\rm D}_{0}(\mathfrak{h})$ is a
scalar. Compared to the second-quantized form of the bare Hamiltonian $H^{\rm
B}(\mathfrak{h})$ in the $\mathfrak{h}$-generated active space

$H^{\rm B}(\mathfrak{h})=\Gamma^{\rm
B}_{0}(\mathfrak{h})+\sum_{PQ}h^{P}_{Q}(\mathfrak{h})a_{P}^{\dagger}a_{Q}+\frac{1}{4}\sum_{P,Q,R,S}v^{PQ}_{RS}(\mathfrak{h})a_{P}^{\dagger}a_{Q}^{\dagger}a_{S}a_{R}$
(11)

the $H^{\rm D}(\mathfrak{h})$ operator contains not only higher-than-pairwise
interactions but its matrix representation in the active space is denser than
the analogous matrix representation of the $H^{\rm B}(\mathfrak{h})$ operator.

Due to the non-commuting character of components defining the $\sigma_{\rm
ext}(\mathfrak{h})$ operator, one has to rely on the finite-rank commutator
expansions, i.e.,

$H^{\rm D}(\mathfrak{h})\simeq(P+Q_{\rm
int}(\mathfrak{h}))(H+\sum_{i=1}^{l}\frac{1}{i!}[\ldots[H,\sigma_{\rm
ext}(\mathfrak{h})],\ldots],\sigma_{\rm ext}(\mathfrak{h})]_{i}(P+Q_{\rm
int}(\mathfrak{h}))\;.$ (12)

where $[\ldots[...]\ldots]_{i}$ designates rank-$i$ commutator, to obtain
approximate many-body form of the $H^{\rm D}(\mathfrak{h})$ operator. The
length $l$ is not only determined by the accuracy requirements but also by the
associated cost. In typical applications in quantum chemistry including
contributions from single, double, and part of the triple commutators provide
a good compromise between the cost and accuracies (see Ref.[40] for more
detailed discussion). In Ref.[27], we provided an overview of algebraic
techniques to derive the second quantized form of Hamiltonian (10). Since the
discussed formalisms are of single-reference type, the most compressed form of
the approximate effective Hamiltonians are provided by utilizing particle-hole
formalism.

The second and third classes of approximations are associated with the rank of
many-body effects included in Eq.(10) and the approximation protocol for the
$\sigma_{\rm ext}(\mathfrak{h})$ operator. In all our applications (which will
be also used throughout this paper) we include one- and two-body effects
defined by $g^{P}_{Q}$ and $k^{PQ}_{RS}$ tensors whereas

$\sigma_{\rm ext}(\mathfrak{h})\simeq T_{\rm ext}^{\rm
CCSD}(\mathfrak{h})-T_{\rm ext}^{\rm CCSD}(\mathfrak{h})^{\dagger}\;,$ (13)

where $T_{\rm ext}^{\rm CCSD}(\mathfrak{h})$ are external part of the cluster
operator $T$ obtained in the CCSD calculations.

Even the most basic non-trivial approximation entails managing a combinatorial
explosion in the quantity of contributing diagrams, equivalent to thousands of
Hugenholtz diagrams associated with the incorporation of higher-rank
commutators. The manual derivation of all diagrams is not only susceptible to
errors but is also excessively time-consuming. To address this issue, we have
created a symbolic algebra system named SymGen that can automatically derive
tensor expressions corresponding to intricate contractions in Wick’s Theorem.

Due to the significant numerical overhead involved in multi-dimensional tensor
contractions, the output of the software tool SymGen has been integrated with
the Tensor Algebra for Many-body Methods (TAMM) parallel tensor library [42]
to address the high numerical scaling associated with the above approximation
scheme. This scaling increases rapidly with the inclusion of higher many-body
effects in the effective Hamiltonian and $\sigma_{\rm ext}(\mathfrak{h})$
operator, not to mention the cost of higher-rank commutators. Therefore, there
is a clear need for new computational paradigms to advance effective
interaction-based computational models.

##  IV Machine learning for evaluating downfolded Hamiltonians

The high cost associated with the evaluation of the one and two-body
interactions in the second quantization representation of the effective
Hamiltonian calls for an alternative way to calculate these interactions
efficiently. Because the bare interactions are relatively easy to compute, and
the one and two-body interactions appearing in the effective Hamiltonian can
be viewed as screened (or dressed) interactions that are related to the bare
interactions by a nonlinear map, it may be possible to construct a neural
network (NN) representation of such a map denoted by (14) and train such a NN
with computed bare interactions and a subset of effective interactions we can
afford to compute.

$H^{\rm B}(\mathfrak{h})\xrightarrow{\text{NN}}H^{\rm D}(\mathfrak{h})\;.$
(14)

In particular, we may use bare and effective interactions computed for certain
molecular geometries as the training data and use the trained NN to predict
effective interactions from bare interactions at other molecular geometries.
We can also view the NN as a computational tool to complete the three- or
five-dimensional (two- or four-dimensional effective interaction tensors with
additional dimension corresponding to set of molecular geometries), from a
subset of the interaction tensor elements. We will illustrate the performance
of NN techniques by using the example of two-body effective interactions.

###  IV.1 Physics-Informed Neural Network

Although it is possible to use a generic deep NN that takes orbital indices
such as $P$, $Q$, $R$, $S$ as well as the numerical value of $v^{PQ}_{RS}$ as
the input and predicts the value of $k^{PQ}_{RS}$ as the output, such a NN is
likely to be less effective than a carefully designed and physics-informed NN
[26] (for the sake of notational simplicity, we will sometimes skip the sub-
algebra index $\mathfrak{h}$ in the interaction tensors). PINNs have recently
proven effective at solving partial differential equations that arise in the
physical sciences and engineering applications [26, 43]. Specifically, PINNs
incorporate physical laws, such as conservation principles or equations of
motion, directly into the training process, enhancing the model generalization
[44]. This approach has been successfully applied across various domains,
including fluid dynamics, quantum mechanics, and materials science [45, 46].
Because one and two-body interactions are naturally represented as integrals
of molecular orbitals connected by an interaction kernel, the PINN we
construct will preserve such a structure. In particular, we will use a NN to
represent screened molecular orbitals and a separate matrix to represent an
effective interaction kernel. These components are connected as shown in
Figure 3 to predict the value of $k^{PQ}_{RS}(\mathfrak{h},\mathbf{R})$ for a
given geometry $\mathbf{R}$ (see also the orbital variant
$(pq|rs)_{D}(\mathbf{R})$ of these tensors defined in next subsection). This
approach allows us to reduce the number of parameters in the PINN compared to
a generic NN. We will describe how such a NN, which we call VNet, is
constructed and trained.

It is important to note that a PINN may not be universally applicable to all
types of chemical transformations. Instead, we intend to train various NNs
based on the nature of the problem to be solved, which is typically determined
by the size of the active spaces and the level of approximation used to
evaluate effective Hamiltonians. For instance, for a wide range of chemical
transformations that involve breaking a single bond, we will utilize NNs
trained for relatively small active spaces. In contrast, for problems
involving the breaking of multiple bonds, larger active spaces will
necessitate the use of NNs trained for such purposes.

The cost and fidelity of NNs are contingent upon the active space size and
rank of the many-body interactions. As a proof of concept, in this paper, we
focus on two-body interactions (described by the $k^{PQ}_{RS}(\mathfrak{h})$
tensor) and active spaces, where all occupied orbitals are defined as active.
It also epitomizes the challenging situation for NN algorithms because entries
$k^{PQ}_{RS}(\mathfrak{h})$ elements where $P$, $Q$, $R$, and $S$ indices
correspond to low-lying orbitals can assume large values.

We will use the Mulliken-type notation for the orbital representation
$(p|q)_{B}$/$(pq|rs)_{B}$ and $(p|q)_{D}$/$(pq|rs)_{D}$ for the one- and two-
body interaction tensors of bare and effective Hamiltonians, respectively,
where indices $p$, $q$, $r$, $s$ designate the active orbitals.

The $(p,q,r,s)$th element of the bare two-body interaction tensor is defined
by

$(pq|rs)_{B}=\int d\mathbf{r}^{\prime}\int
d\mathbf{r}\phi_{p}(\mathbf{r})\phi_{q}(\mathbf{r})\frac{1}{|\mathbf{r}-\mathbf{r}^{\prime}|}\phi_{r}(\mathbf{r}^{\prime})\phi_{s}(\mathbf{r}^{\prime}),$
(15)

where $\phi_{p}(\mathbf{r})$ is the $p$th active molecular orbital. When a
Gaussian basis set is used to represent $\phi_{p}$’s, the integral in (15) can
be performed analytically. A similar integral can be written for the one-body
interaction. For brevity, we will focus on the two-body interaction below.

Our work makes the assumption that the effective two-body interaction for the
downfolded Hamiltonian $H^{\rm D}(\mathfrak{h})$ can also be written as a
double integral except that the bare Coulomb interaction kernel
$1/|\mathbf{r}-\mathbf{r}^{\prime}|$ in (15) should be replaced by a dressed
(effective) interaction kernel $W^{D}(\mathbf{r},\mathbf{r}^{\prime})$, and
the active molecular orbitals $\\{\phi_{p}\\}$ may need to be modified also to
account for change of symmetry in the effective two-body interaction tensor
and the nature of the effective interaction.

In general, we do not know the analytical form of
$W^{D}(\mathbf{r},\mathbf{r}^{\prime})$. Although it may be possible to
construct a parameterized analytic or even a neural network model for such a
kernel and determine the model parameters through a fitting procedure that
minimizes the difference between a subset of the effective interaction tensor
elements predicted by the model and those computed from the downfolding
theory, it is likely to be difficult to evaluate the double integral that
defines the effective interaction accurately using such an ansatz.

###  IV.2 VNet Neural Network

In this work, we propose an alternative approach schematically depicted in
Fig. 2. This approach consists of two steps

  1. 1.

In the first step, we take a bare two-body interaction tensor $(pq|rs)_{B}$
and try to decompose it as

$(pq|rs)_{B}=\left[\phi_{p}\odot\phi_{q}\right]^{T}W^{B}\left[\phi_{r}\odot\phi_{s}\right],$
(16)

where $\phi_{p}$, $\phi_{q}$, $\phi_{r}$, $\phi_{s}$ are vectors of length
$\ell$, chosen from a set of $N_{\rm orb}^{\rm(act)}$ vectors
$\\{\phi_{1},...,\phi_{N_{\rm orb}^{\rm(act)}}\\}$, $W^{B}$ is an
$\ell\times\ell$ symmetric matrix, and $\phi_{p}\odot\phi_{q}$ denotes the
elementwise prodcut of $\phi_{p}$ and $\phi_{q}$. We perform the decomposition
by minimizing the difference between the left and right-hand sides of (16)
with respect to $\\{\phi_{1},\phi_{2},...,\phi_{N_{\rm orb}^{\rm(act)}}\\}$,
and the unique matrix elements in $W^{B}$. The total number of parameters
involved in the optimization problem is $N_{\rm
orb}^{\rm(act)}\times\ell+(\ell+1)\ell/2$ in this case. Although we can
possibly achieve an accurate fit by making $\ell$ sufficiently large, more
flexibility is provided if $\phi$ is represented by a neural network (NN)
parameterized by a set of weights and biases contained in a vector
$\mathbf{\theta}$. The neural network takes the orbital index as one of the
inputs. Because molecular orbitals depend on the geometry of the molecule, the
neural network should take the geometry represented by a parameter
$\mathbf{R}$ as the input also. Moreover, we can introduce another parameter
$\gamma\in\\{0,1\\}$ to account for special symmetry of the $(pq|rs)_{B}$
tensor. We will denote this NN by
$\phi_{\mathbf{\theta}}(i,\mathbf{R},\gamma)$. The output of the NN contains
$\ell$ real numbers. The architecture of this NN is shown in Figure 3. For
bare interaction tensors, we always choose $\gamma=0$ so that the 8-fold
symmetry of the tensor is preserved.

We use the bare interaction tensors for several different geometries
$\\{\mathbf{R}\\}$ to train the NN to obtain an initial set of parameters
$\mathbf{\theta}$ and the $W^{B}$ matrix. For each geometry, the training data
consists of elements of $(pq|rs)_{B}$ that belong to one nonsymmetric unit.
Note that the algebraic decomposition of the bare two-body interaction tensor
in (16) eliminates physical spatial coordinate $\mathbf{r}$ of a molecular
orbital and replaces it with a latent spatial coordinate.

  2. 2.

In the second step, we use the effective two-body interaction tensor elements
$(pq|rs)_{D}(\mathbf{R})$ associated with a few selected geometries
$\mathbf{R}_{i}$, $i=1,2,...,n_{\rm ref}$, to refine the NN parameters
contained in $\mathbf{\theta}$ and obtain an approximation to the effective
interaction kernel $W^{D}$ by solving the following optimization problem

$\min_{(W^{D})^{T}=W^{D},\mathbf{\theta}}\sum_{i=1}^{n_{\rm
ref}}\sum_{\\{p,q,r,s\\}\in\mathcal{S}}\biggl{|}(pq|rs)_{D}(\mathbf{R}_{i})-\left[\phi_{\mathbf{\theta}}(p,\mathbf{R}_{i},0)\odot\phi_{\mathbf{\theta}}(q,\mathbf{R}_{i},1)\right]^{T}W^{D}\left[\phi_{\mathbf{\theta}}(r,\mathbf{R}_{i},0)\odot\phi_{\mathbf{\theta}}(s,\mathbf{R}_{i},1)\right]\biggr{|}^{2},$
(17)

where $\mathcal{S}$ is a subset of the full index set for the $N_{\rm
orb}^{\rm(act)}\times N_{\rm orb}^{\rm(act)}\times N_{\rm
orb}^{\rm(act)}\times N_{\rm orb}^{\rm(act)}$ tensor $(pq|rs)_{D}(\mathbf{R})$
that belongs to a nonsymmetric subunit, and $n_{\rm ref}$ is the number of
selected reference geometries. Note that $(pq|rs)_{D}$ only has a 4-fold
symmetry. Calculating $(pq|rs)_{D}$, in contrast to $(pq|rs)_{B}$, is
associated with a significant numerical overhead.

It is important to point out that both $(pq|rs)_{B}(\mathbf{R})$ and
$(pq|rs)_{D}(\mathbf{R})$ are geometry dependent. This $\mathbf{R}$ dependency
translates into the dependence of $\phi_{\mathbf{\theta}}$ on $\mathbf{R}$.
However, both the bare and effective interaction kernels $W^{B}$ and $W^{D}$
are independent of $\mathbf{R}$. As a result, we can choose $\mathcal{S}$ to
be the entire index subset for a nonsymmetric subunit of
$(pq|rs)_{D}(\mathbf{R})$ and determine $W^{D}$ from the effective two-body
interaction tensor for a few geometric configurations associated with a few
different $\mathbf{R}_{i}$’s. This $W^{D}$ can be combined with the refined NN
$\phi_{\mathbf{\theta}}(i,\mathbf{R}^{\prime},\gamma)$ to predict effective
two-body interaction tensors at other geometries defined by a different set of
$\mathbf{R}^{\prime}$.

Figure 2:  Schematic representation of PINNs for evaluating effective
interactions. Training data is aware of symmetries of bare and interaction
tensors. Moreover, the PINNs are aware of the problem’s complexity, reflected
by the size of the active space used. $N_{\rm orb}$ and $N^{\rm(act)}_{\rm
orb}$ designate the total number of orbitals and the number of active
orbitals, respectively. It is assumed that $N^{\rm(act)}_{\rm orb}\ll N_{\rm
orb}$.

In the forthcoming analysis, we aim to present a proof of concept for the
applicability of the proposed algorithms in small-sized benchmark systems,
namely H2O and HF molecules. Our approach will showcase the potential benefits
of training neural networks (NNs) using datasets acquired from the immediate
vicinity of equilibrium geometries, which is motivated by the fact that SR-CC
formalisms yield the most accurate results in such regions. Subsequently, we
will employ the NN-predicted effective Hamiltonians to explore the other parts
of the configurational space of the molecular systems, for example, those
related to the stretching or breaking of a single bond.

Figure 3:  VNet. (a) The neural network representation of the molecular
orbital $p$. (b) The proposed ansatz to approximate the effective two-body
interaction tensor. (c) The proposed paradigm for pre-training and finetuning.

##  V Results and discussions

We now demonstrate the effectiveness using the bare and effective two-body
interactions of two molecules H2O and HF at a few selected reference
geometries to train VNet and predict (complete) the two-body effective
interactions of these molecules at other geometries. For both H2O and HF, we
first computed the complete active space self-consistent field (CASSCF)
reference with the active space consisting of the eight lowest-energy orbitals
and the cc-pVTZ basis set. The final effective interactions were then computed
using the A7 DUCC approximation that is described in Ref. [40].

For the H2O molecule, we used the bare two-body interaction tensors associated
with different relative O-H bond lengths corresponding to a single bond
breaking, i.e., the ratio between the O-H bond length ($R_{\rm OH}$) and the
O-H bond length at equilibrium ($R_{\rm eq}$; $R_{\rm eq}={\rm
1.84345\;a.u.}$), ranging from 1.1 to 2.5 (a total of 66 different geometries)
to train the VNet to obtain both $W^{B}$ and an initial set of parameters
$\theta$ for $\phi_{\mathbf{\theta}}$. We used the effective two-body
interaction tensors associated with relative bond lengths 1.15, 1.45, 1.95,
and 2.45 to finetune the parameter $\theta$ in the NN representations of the
molecular orbitals $\phi_{\mathbf{\theta}}$ and to obtain an effective
interaction kernel matrix $W^{D}$. The NN used to represent
$\phi_{\mathbf{\theta}}$ consists of $N_{\rm orb}^{\rm(act)}+2$ input neurons.
This includes $N_{\rm orb}^{\rm(act)}$ neurons in a one-hot encoding to
account for the orbital index, one neuron for the geometry $\mathbf{R}$, and
one neuron for the symmetry index $\gamma$. The NN contains 3 hidden layers
with 200 neurons in each layer. The output of the NN consists of 300 neurons
and we represent $W^{B}$ and $W^{D}$ as $300\times 300$ symmetric matrices.
The implementation details, such as optimization specifics, are provided in
Appendix A.

In Figure 6, we plot the maximum absolute error (MAE) and mean square error
(MSE) of the two-body effective interaction tensor elements for several
geometries associated with different relative O-H bond lengths (orange curve
marked with dots). We observe that the MAEs for all geometries are around
$10^{-3}$, which are sufficiently small. Not surprisingly, the MAEs of the
predicted effective tensors associated with geometries close to the geometries
used for training (marked by green crosses) are much smaller. As a comparison,
we also plot the maximum difference between an element of the bare two-body
interaction tensor and the corresponding element in the effective two-body
tensor. We can clearly see this difference is almost two orders of magnitude
larger.

To further evaluate the accuracy of the neural network-derived Hamiltonians,
we computed the reference, correlation, and total energies for the
Hamiltonians and compared them against the expected DUCC Hamiltonians. The
one-body contribution to the Hamiltonian was fixed to the corresponding
Hamiltonian (either bare or downfolded). As seen in Fig. 4, the DUCC
Hamiltonians well capture the correlation along the bond-breaking coordinate
and are a significant improvement over a simple truncation of the Hamiltonian
to the active space, which is illustrated by the bare results. For all
geometries, the neural network generated Hamiltonians provide energies close
to the expected DUCC results. In fact, the neural network Hamiltonians capture
over 97% of the correlation energy for all geometries. The largest source of
error in the total energies is attributed to the two-body interactions of the
largest magnitude, which are represented by the reference energy ($(ii|jj)$
and $(ij|ji)$, where $i$ and $j$ are occupied spin-orbital indices). For these
large-magnitude interactions, even small relative errors can lead to
measurable errors in energies. We plan to investigate this further in later
studies to find approaches to better address this class of interactions.

Figure 4: Top Panel: Comparison of total energies of H2O obtained from the
RHF, DUCC, and neural network Hamiltonians. Bottom Panel: Analysis of the
energy components between the neural network Hamiltonians and corresponding
DUCC Hamiltonians. The x-axis is the $\rm O-H$ bond distance for the stretched
bond along the single $\rm O-H$ bond breaking pathway and is measured relative
to the equilibrium bond length (${\rm R}_{\rm eq}={\rm 1.84345a.u.},{\rm
H-O-H}=110.6^{\circ}$).

For the HF molecule, we used the bare two-body interaction tensors associated
with different relative H-F bond lengths ($R_{\rm HF}/R_{\rm eq}$; $R_{\rm
eq}={\rm 1.7325a.u.}$) ranging from 0.85 to 2.0 (a total of 215 different
geometries) to train the VNet to obtain both $W^{B}$ and an initial set of
parameters $\theta$ for $\phi_{\mathbf{\theta}}$. We used the effective two-
body interaction tensors associated with relative bond lengths 0.95, 1.35,
1.65, and 1.95 to finetune the parameter $\theta$ in the NN representations of
the molecular orbitals $\phi_{\mathbf{\theta}}$ and to obtain an effective
interaction kernel matrix $W^{D}$. The NN architecture remains the same for
HF.

In Figure 7, we plot the MAE and MSE of the two-body effective interaction
tensor elements for several geometries associated with different relative H-F
bond lengths (orange curve marked with dots). We observe that the MAEs for all
geometries are around $10^{-3}$, which are sufficiently small. Not
surprisingly, the MAEs of the predicted effective tensors associated with
geometries close to the geometries used for training (marked by green crosses)
are smaller. As a comparison, we also plot the maximum difference between an
element of the bare two-body interaction tensor and the corresponding element
in the effective two-body tensor. We can clearly see this difference is almost
two orders of magnitude larger.

As with H2O, we computed the reference, correlation, and total energies for
the Hamiltonians of HF and compared them against the expected DUCC
Hamiltonians to evaluate the accuracy of the neural network-derived
Hamiltonians. Fig. 5 shows that the DUCC Hamiltonians well capture the
correlation and are a significant improvement over the truncated bare
Hamiltonian in the active space once again. It is important to emphasize that
any error in DUCC can be systematically improved through the underlying
approximations as shown in Ref. 47. For all geometries, the neural network
generated Hamiltonians provide energies close to the expected DUCC results
capturing over 99% of the correlation energy for all geometries. Again, the
largest source of error in the total energies is attributed to the two-body
interactions of the largest magnitude which are represented by the reference
energy. Compared with H2O, both the reference and correlation energies better
agree with the expected DUCC results in general.

Figure 5: Top Panel: Comparison of total energies of HF obtained from the RHF,
DUCC, and neural network Hamiltonians. Bottom Panel: Analysis of the energy
components between the neural network Hamiltonians and corresponding DUCC
Hamiltonians. The x-axis is the H–F bond distance relative to the equilibrium
bond length (${\rm R}_{\rm eq}={\rm 1.7325a.u.}$).

It takes $(1.80\pm 0.12)\times 10^{-3}$ seconds (averaged over 10 independent
trials) to predict the effective interaction for H2O in a particular geometry
using the VNet. This is a drastic saving over the relatively costly
downfolding procedure which has $\mathcal{O}(\mathscr{N}^{6})$ scaling with
the system size $\mathscr{N}$ for the underlying CCSD calculation and similar
scaling for the actual downfolding procedure.

Figure 10 illustrates the necessity to use the bare interaction to pretrain
VNet to obtain a good initial guess of the VNet parameters. Without using the
bare interaction to pretrain VNet, the MSE of the VNet loss function plateaus
around $10^{-5}$. In contrast, after using the bare interation to pretrain the
VNet first and then using the DUCC effective interaction to finetune to VNet
parameters, the MSE of the loss function reaches $10^{-7}$, which is two
orders of magnitude lower than the non-pretrained version.

In Figure 11, we compare the effectiveness of VNet and that of a generic fully
connected neural network with Fourier features [48] that takes the indices of
$(pq|rs)_{D}$ as well as the geometry $\mathbf{R}$ as the input and predicts
the numerical value of the corresponding tensor element as the output. The
generic NN with Fourier features is trained with the effective interaction
tensor at a few selected geometries. The details are provided in Appendix B.
As we can see from this figure, the training mean square error (MSE) can be as
low as $10^{-11}$. However, when this NN is used to infer effective
interaction tensor elements outside of the training geometries, the MSE can be
as large as $10^{-1}$. Clearly, the generic NN, which does not take advantage
of the special structure of the problem, is overfitted to the training data
and has trouble completing effective interaction tensors for untrained
geometries.

In addition to producing accurate approximations to effective interaction
tensors for different geometries, VNet also allows us to study the effect of
screening by comparing the bare and effective interaction kernel matrices
$W^{B}$ and $W^{D}$ which are independent of the geometry.

We computed the eigendcomposition of $W^{B}$, i.e.,

$W^{B}=Z\Lambda^{B}Z^{T},$ (18)

where the matrix $Z$ contains the eigenvectors of $W^{B}$, and the diagonal
matrix $\Lambda^{B}$ contains the corresponding eigenvalues
$\varepsilon_{1}^{B},\varepsilon_{2}^{B},...,\varepsilon_{\ell}^{B}$.

Using the eigenvector matrix $Z$ associated with $W^{B}$, we compute
$Z^{T}W^{D}Z$. Figures 8 and 9 show that this matrix is nearly diagonal for
both H2O and HF, which means that $W^{D}$ and $W^{B}$ share the almost same
set of eigenvectors. The diagonal elements of $Z^{T}W^{D}Z$ are close to the
eigenvalues of $W^{D}$. Therefore, it is meaningful to compare the eigenvalues
of $W^{B}$ and $W^{D}$ in order to understand the screening function of the
effective interaction kernel.

In Figures 8 and 9, we plotted the ratio
$\tau=(\varepsilon^{D}-\varepsilon^{B})/\varepsilon^{B}$ for all eigenvalues
of $W^{D}$ and the corresponding eigenvalues of $W^{B}$. These ratios are
sorted in increasing order before they are plotted. Note that the horizontal
axis does not correspond to a physical spatial variable. It is a latent
dimension created in the VNet representation of the molecular orbital. The
inset plots were obtained by excluding data points that correspond to nearly
zero $\varepsilon^{B}$ values. From these plots, we can see that the ratio
$\tau$ seems to follow a scaled tangent which suggests that $\varepsilon^{D}$
is related to $\varepsilon^{B}$ by
$\varepsilon^{D}_{i}=\varepsilon^{B}_{i}(1+\beta\tan(\alpha(i-i_{c})))$ for
some parameters $\alpha$, $\beta$ and $i_{c}$.

Figure 6: Maximum absolute error (MAE) and mean square error (MSE) between the
true two-body effective interaction tensor and the predicted tensor for the
H2O molecule.  Figure 7:  Maximum absolute error (MAE) and mean square error
(MSE) between the true two-body effective interaction tensor and the predicted
tensor for the HF molecule.

(a) $Z^{T}W^{D}Z$

(b) $(\varepsilon^{D}-\varepsilon^{B})/\varepsilon^{B}$

Figure 8: Analyzing the bare and effective interaction kernel matrices $W^{B}$
and $W^{D}$ for the H2O molecule. (a) $Z^{T}W^{D}Z$, where $Z$ denotes the
eigenvectors of $W^{B}$. (b) The ratio of all eigenvalues of $W^{D}$ to the
corresponding eigenvalues of $W^{B}$.

(a) $Z^{T}W^{D}Z$

(b) $(\varepsilon^{D}-\varepsilon^{B})/\varepsilon^{B}$

Figure 9: Analyzing the bare and effective interaction kernel matrices $W^{B}$
and $W^{D}$ for the HF molecule. (a) $Z^{T}W^{D}Z$, where $Z$ denotes the
eigenvectors of $W^{B}$. (b) The ratio of all eigenvalues of $W^{D}$ to the
corresponding eigenvalues of $W^{B}$.

(a) The HF molecule (b) The H2O molecule

Figure 10: Comparison of NN training on DUCC with and without bare pretrained
orbitals on the molecules (a) HF and (b) H2O.

(a) The HF molecule (b) The H2O molecule

Figure 11: Comparison of training and testing performance for regression using
Vnet and Coordinate-based MLP methods on the molecules (a) HF and (b) H2O. The
Fourier feature is a coordinate-based MLP approach that uses a special fully
connected network with input coordinates and outputs the corresponding value.
Although the training error of the Fourier feature method is nearly zero, the
testing error is relatively large, indicating a severe overfitting problem.

##  VI Conclusion

The CC downfolding theory opens the door to studying large molecules by
constructing an effective Hamiltonian in a small-dimensionality active space.
Although it is possible to construct such a Hamiltonian for a few selected
geometries, performing such calculations for a large number of geometries as
would required in a dynamics simulation is costly. In this work, we propose an
ML-based approach to predict the effective interactions of the downfolded
Hamiltonian by a special physics-informed neural network called VNet. The VNet
is trained with several bare interaction tensors as well as a few effective
interaction tensors associated with a few selected geometries. Using HF and
H2O molecules as examples, we have demonstrated that VNet can be effectively
trained to accurately predict the effective interactions for many geometries
not included in the training dataset. For these problems, the VNet model has
188,050 trainable parameters. For the HF case, training on the 215 Bare
tensors for 5,000 epochs takes 4.4 hours. Finetuning VNet using 4 DUCC tensors
for 500 epochs takes 10 minutes. Once VNet is trained, it only takes $(1.80\pm
0.12)\times 10^{-3}$ seconds (averaged over 10 independent trials) to predict
one effective interaction tensor. This presents a powerful alternative to
combat the $\mathcal{O}(\mathscr{N}^{6})$ scaling associated with the
downfolding procedure.

Not only can VNet allow us to accurately predict effective interactions
rapidly, but it can also be used to elucidate the nature of effective
interactions in terms of the screening effect of the effective interaction
kernel. In addition, the efficacy of the VNet framework was illustrated in a
complex scenario involving active orbitals, comprising a subset of valence
virtual orbitals and all occupied orbitals. The latter group of orbitals
exhibits diverse localization patterns and corresponds to different energy
regimes. For future applications, we intend to adopt a more uniform definition
of active orbitals, encompassing only valence occupied and virtual orbitals.

What we have presented here is the first attempt to use PINN-based ML
techniques to model and predict effective interactions for realistic
molecules. The VNet architecture we used in this work is by no means optimal
and can be further improved. More work is also required to investigate the
transferability of this approach to other molecular systems.

###### Acknowledgements.

This material is based upon work supported by the “Transferring exascale
computational chemistry to cloud computing environment and emerging hardware
technologies (TEC4)” project, which is funded by the U.S. Department of
Energy, Office of Science, Office of Basic Energy Sciences, the Division of
Chemical Sciences, Geosciences, and Biosciences (under FWP 82037). N.P.B. also
acknowledges support from the Laboratory Directed Research and Development
(LDRD) Program at Pacific Northwest National Laboratory. This work is also
supported by the U.S. Department of Energy, Office of Science, Office of
Advanced Scientific Computing Research and Office of Basic Energy Science,
Scientific Discovery through Advanced Computing (SciDAC) program under
Contract No. DE-AC02-05CH11231. (S.L. and C.Y.) This work used the
computational resources of the National Energy Research Scientific Computing
(NERSC) center under NERSC Award ASCR-ERCAP m1027 for 2024, which is supported
by the Office of Science of the U.S. Department of Energy under Contract No.
DE-AC02-05CH11231.

##  Appendix A Implementation details

This section will outline the implementation details, including dataset
preparation, neural network architecture, and neural network optimization.

###  A.1 Dataset

For the H2O molecule case, bare tensors with geometry ranging from relative
H-O bond lengths of 1.1 to 2.5 (a total of 66 different bond lengths) were
used to train a baseline VNet model. Then DUCC effective two-body interaction
tensors with relative H-O bond lengths 1.15, 1.45, 1.95 and 2.45 were used to
finetune the baseline VNet model. The finetuned model was then used to predict
the DUCC effective two-body interaction tensors associated with relative H-O
bond lengths ranging from 1.1 to 2.5 in equally spaced increments of 0.1.

For the HF molecule case, bare tensors with geometry ranging from relative H-F
bond lengths of 0.85 to 2.0 (a total of 215 different bond lengths) were used
to train the baseline VNet model. Then DUCC tensors with relative H-F bond
lengths 0.95, 1.35, 1.65, and 1.95 were used to finetune the VNet. The model
was then tested to predict the DUCC tensors with relative H-F bond lengths
ranging from 0.9 to 2.0 in equally spaced increments of 0.1.

###  A.2 Index set for reconstructing DUCC tensor

The size of the index set for a two-body interaction tensor (either DUCC or
Bare) is $(N_{\rm orb}^{\rm(act)})^{4}$. Since the DUCC tensor is 4-fold
symmetric, the total number of unique indices is approximately $(N_{\rm
orb}^{\rm(act)})^{4}/4$. We denote this unique index set by $\mathbb{S}$. A
considerable number of two-body interaction tensor elements are nearly zero
and these elements do not contribute significantly to the total energy
calculation. As a result, we only consider reconstructing the non-zero
components of the DUCC tensor. To identify the indices of the zero tensor
elements, we use the assumption that $(pq|rs)_{D}(\mathbf{R})=0$ if and only
if $(pq|rs)_{B}(\mathbf{R})=0$.

In practice, we set $(pq|rs)_{D}(\mathbf{R})$ to zero if
$\mathrm{mean}\\{(pq|rs)_{B}(\mathbf{R}_{i})\\}<\epsilon_{1}$ and
$\mathrm{std}\\{(pq|rs)_{B}(\mathbf{R}_{i})\\}<\epsilon_{2}$, where the mean
and the standard deviation (std) are taken over all geometries indexed by $i$.
We denote this set of indices $(p,q,r,s)$ associated with the zero DUCC tensor
elements by $\mathbb{P}$. We use $\epsilon_{1}=\epsilon_{2}=10^{-5}$ for all
tests performed in this work. The index set
$\mathbb{P}^{C}=\mathbb{S}\setminus\mathbb{P}$ represent all nonzero tensor
elements we construct. The mean square error (MSE) and maximum absolute error
(MAE) are calculated only on $\mathbb{P}^{C}$.

###  A.3 Network architecture

We approximate molecular orbital by a fully connected neural network, which
can be expressed as
$f(x)=W_{L}\sigma(\cdots(W_{2}\sigma(W_{1}x+b_{1})+b_{2}))+b_{L}$. Here, the
input $x$ is a vector of length $N_{\rm orb}^{\rm(act)}+2$, consisting of
three components. The first component is a scalar representing the geometry
(bond length) of the molecule. The second component is an integer that
describes the type of symmetry of the tensor. The third component is a one-hot
vector $h$ of length $N_{\rm orb}^{\rm(act)}$. It is an occupation
representation of the orbital index $i\in\\{1,2,\cdots,N_{\rm
orb}^{\rm(act)}\\}$, i.e., $h(i)=0$ and $h(j)=0$ if $j\neq i$. The weight
matrices are defined as $W_{1}\in\mathbb{R}^{(N_{\rm orb}^{\rm(act)}+2)\times
M}$, $W_{L}\in\mathbb{R}^{M\times k}$, and $W_{s}\in\mathbb{R}^{M\times M}$
for $s=2,\cdots,L-1$. The bias vectors are $b_{s}$ of length $M$ for
$s=1,\cdots,L-1$ and $b_{L}=k$. The function $\sigma$ is a non-linear
activation function applied element-wise to each entry of a vector.

We use the Sigmoid Linear Unit [49] as the activation function, i.e.,
$\sigma(x)=x/(1+\exp(-x))$. The values of $k$, $M$, and $L$ are chosen as 300,
200, and 4, respectively.

###  A.4 Optimization

During the training phase, the VNet parameters are optimized using the Adam
optimizer, a variant of the stochastic gradient descent method, for 5000
epochs. The learning rate is initialized at 0.001 and follows a cosine decay
schedule. In the finetuning phase, these parameters are further optimized for
500 epochs. The learning rate is set to 0.0002 for the H2O molecule case
(0.0005 for the HF case), which is smaller than the rate used in the initial
training phase. The learning rate is also adjusted according to a cosine decay
schedule.

##  Appendix B Comparison with generic NN method

To demonstrate the advantage of a PINN over a generic NN for completing the
effective two-body interaction, we compare the proposed VNet with a
coordinate-based multi-layer perceptron (MLP) [50, 51, 48]. A coordinate-based
MLP is a fully connected neural network architecture often designed to learn
continuous functions in several spatial coordinates. Unlike traditional MLPs
that operate on abstract data representations, coordinate-based MLPs take a
set of coordinates as the input and produce a scalar or vector as the output,
approximating the function value, e.g., color or density, defined at this set
of coordinates. The Fourier feature method (FFM) [48], a popular coordinate-
based MLP that adds a Fourier feature with random frequencies to the input
layer, is adopted for comparison in the experiments. The input to the FFM is
the index $(p,q,r,s)$ of a particular two-body interaction tensor element as
well as the geometry $\mathbf{R}$. The output of FFM approximates
$(pq|rs)_{D}(\mathbf{R})$. We used the same training setup as that used for
VNet except that the number of training epochs was set to 200. In Figure 11,
we compare the training and testing errors of VNet with those of the
coordinate-based MLP methods for both H2O and HF molecules. While the
coordinate-based MLP achieves almost zero loss function values on the training
data, the testing error on other geometries is relatively large. For VNet, the
loss function achieves nearly the same values on both the training and testing
datasets. Therefore, the physics-informed model provided by VNet appears to be
superior to a generic coordinate-based MLP NN.

## References

  * Jiménez-Hoyos and Scuseria [2015] C. A. Jiménez-Hoyos and G. E. Scuseria, Cluster-based mean-field and perturbative description of strongly correlated fermion systems: Application to the one-and two-dimensional hubbard model, Physical Review B 92, 085101 (2015).
  * Braunscheidel  _et al._ [2024] N. M. Braunscheidel, A. Bachhar, and N. J. Mayhall, Accurate and interpretable representation of correlated electronic structure via tensor product selected ci, Faraday Discussions (2024).
  * Gray and Chan [2024] J. Gray and G. K.-L. Chan, Hyperoptimized approximate contraction of tensor networks with arbitrary geometry, Physical Review X 14, 011009 (2024).
  * Hirata  _et al._ [2004] S. Hirata, R. Podeszwa, M. Tobita, and R. J. Bartlett, Coupled-cluster singles and doubles for extended systems, The Journal of chemical physics 120, 2581 (2004).
  * Booth  _et al._ [2013] G. H. Booth, A. Grüneis, G. Kresse, and A. Alavi, Towards an exact description of electronic wavefunctions in real solids, Nature 493, 365 (2013).
  * Zhang and Grüneis [2019] I. Y. Zhang and A. Grüneis, Coupled cluster theory in materials science, Frontiers in Materials 6, 123 (2019).
  * Liao  _et al._ [2021] K. Liao, T. Schraivogel, H. Luo, D. Kats, and A. Alavi, Towards efficient and accurate ab initio solutions to periodic systems via transcorrelation and coupled cluster theory, Physical Review Research 3, 033072 (2021).
  * Wang and Berkelbach [2020] X. Wang and T. C. Berkelbach, Excitons in solids from periodic equation-of-motion coupled-cluster theory, Journal of Chemical Theory and Computation 16, 3095 (2020).
  * Dean and Hjorth-Jensen [2004] D. Dean and M. Hjorth-Jensen, Coupled-cluster approach to nuclear physics, Physical Review C 69, 054320 (2004).
  * Kowalski  _et al._ [2004] K. Kowalski, D. Dean, M. Hjorth-Jensen, T. Papenbrock, and P. Piecuch, Coupled cluster calculations of ground and excited states of nuclei, Physical review letters 92, 132501 (2004).
  * Hagen  _et al._ [2014] G. Hagen, T. Papenbrock, M. Hjorth-Jensen, and D. J. Dean, Coupled-cluster computations of atomic nuclei, Reports on Progress in Physics 77, 096302 (2014).
  * Hu  _et al._ [2022] B. Hu, W. Jiang, T. Miyagi, Z. Sun, A. Ekström, C. Forssén, G. Hagen, J. D. Holt, T. Papenbrock, S. R. Stroberg, _et al._ , Ab initio predictions link the neutron skin of 208pb to nuclear forces, Nature Physics 18, 1196 (2022).
  * Riplinger and Neese [2013] C. Riplinger and F. Neese, An efficient and near linear scaling pair natural orbital based local coupled cluster method, The Journal of chemical physics 138 (2013).
  * Riplinger  _et al._ [2013] C. Riplinger, B. Sandhoefer, A. Hansen, and F. Neese, Natural triple excitations in local coupled cluster calculations with pair natural orbitals, The Journal of chemical physics 139 (2013).
  * Pinski  _et al._ [2015] P. Pinski, C. Riplinger, E. F. Valeev, and F. Neese, Sparse maps—a systematic infrastructure for reduced-scaling electronic structure methods. i. an efficient and simple linear scaling local mp2 method that uses an intermediate basis of pair natural orbitals, The Journal of chemical physics 143 (2015).
  * Riplinger  _et al._ [2016] C. Riplinger, P. Pinski, U. Becker, E. F. Valeev, and F. Neese, Sparse maps—a systematic infrastructure for reduced-scaling electronic structure methods. ii. linear scaling domain based pair natural orbital coupled cluster theory, The Journal of chemical physics 144 (2016).
  * Kowalski [2018] K. Kowalski, Properties of coupled-cluster equations originating in excitation sub-algebras, The Journal of Chemical Physics 148 (2018).
  * Shee  _et al._ [2024] A. Shee, F. M. Faulstich, B. Whaley, L. Lin, and M. Head-Gordon, A static quantum embedding scheme based on coupled cluster theory, arXiv preprint arXiv:2404.09078 (2024).
  * Huang  _et al._ [2023] R. Huang, C. Li, and F. A. Evangelista, Leveraging small-scale quantum computers with unitarily downfolded hamiltonians, PRX Quantum 4, 020313 (2023).
  * Kowalski and Bauman [2023] K. Kowalski and N. P. Bauman, Quantum flow algorithms for simulating many-body systems on quantum computers, Physical Review Letters 131, 200601 (2023).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Xia and Kais [2018] R. Xia and S. Kais, Quantum machine learning for electronic structure calculations, Nature communications 9, 4195 (2018).
  * Pfau  _et al._ [2020] D. Pfau, J. S. Spencer, A. G. Matthews, and W. M. C. Foulkes, Ab initio solution of the many-electron schrödinger equation with deep neural networks, Physical Review Research 2, 033429 (2020).
  * Hermann  _et al._ [2020] J. Hermann, Z. Schätzle, and F. Noé, Deep-neural-network solution of the electronic schrödinger equation, Nature Chemistry 12, 891 (2020).
  * Bassi  _et al._ [2024] H. Bassi, Y. Zhu, S. Liang, J. Yin, C. C. Reeves, V. Vlček, and C. Yang, Learning nonlinear integral operators via recurrent neural networks and its application in solving integro-differential equations, Machine Learning with Applications 15, 100524 (2024).
  * Karniadakis  _et al._ [2021] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Physics-informed machine learning, Nature Reviews Physics 3, 422 (2021).
  * Bauman  _et al._ [2019] N. P. Bauman, E. J. Bylaska, S. Krishnamoorthy, G. H. Low, N. Wiebe, C. E. Granade, M. Roetteler, M. Troyer, and K. Kowalski, Downfolding of many-body hamiltonians using active-space models: Extension of the sub-system embedding sub-algebras approach to unitary coupled cluster formalisms, The Journal of chemical physics 151 (2019).
  * Piecuch  _et al._ [1993] P. Piecuch, N. Oliphant, and L. Adamowicz, A state-selective multireference coupled-cluster theory employing the single-reference formalism, J. Chem. Phys. 99, 1875 (1993).
  * Piecuch and Adamowicz [1994] P. Piecuch and L. Adamowicz, State-selective multireference coupled-cluster theory employing the single-reference formalism: Implementation and application to the h8 model system, J. Chem. Phys. 100, 5792 (1994).
  * Adamowicz  _et al._ [1998] L. Adamowicz, P. Piecuch, and K. B. Ghose, The state-selective coupled cluster method for quasi-degenerate electronic states, Molecular Physics 94, 225 (1998).
  * Piecuch [2010] P. Piecuch, Active-space coupled-cluster methods, Mol. Phys. 108, 2987 (2010).
  * Kinoshita  _et al._ [2005] T. Kinoshita, O. Hino, and R. J. Bartlett, Coupled-cluster method tailored by configuration interaction, The Journal of chemical physics 123 (2005).
  * Veis  _et al._ [2016] L. Veis, A. Antalík, J. Brabec, F. Neese, O. Legeza, and J. Pittner, Coupled cluster method with single and double excitations tailored by matrix product state wave functions, The journal of physical chemistry letters 7, 4072 (2016).
  * Kowalski [2023] K. Kowalski, Sub-system self-consistency in coupled cluster theory, The Journal of Chemical Physics 158 (2023).
  * Purvis and Bartlett [1982] G. D. Purvis and R. J. Bartlett, A full coupled-cluster singles and doubles model: The inclusion of disconnected triples, J. Chem. Phys. 76, 1910 (1982).
  * Kowalski [2021] K. Kowalski, Dimensionality reduction of the many-body problem using coupled-cluster subsystem flow equations: Classical and quantum computing perspective, Physical Review A 104, 032804 (2021).
  * Kowalski and Bauman [2020] K. Kowalski and N. P. Bauman, Sub-system quantum dynamics using coupled cluster downfolding techniques, The Journal of chemical physics 152 (2020).
  * Evangelista  _et al._ [2019] F. A. Evangelista, G. K. Chan, and G. E. Scuseria, Exact parameterization of fermionic wave functions via unitary coupled cluster theory, The Journal of chemical physics 151 (2019).
  * Lee  _et al._ [2018] J. Lee, W. J. Huggins, M. Head-Gordon, and K. B. Whaley, Generalized unitary coupled cluster wave functions for quantum computation, Journal of chemical theory and computation 15, 311 (2018).
  * Bauman and Kowalski [2022] N. P. Bauman and K. Kowalski, Coupled cluster downfolding methods: The effect of double commutator terms on the accuracy of ground-state energies, The Journal of Chemical Physics 156 (2022).
  * Feldmann  _et al._ [2024] R. Feldmann, M. Mörchen, J. Lang, M. Lesiuk, and M. Reiher, Complete active space iterative coupled cluster theory, arXiv preprint arXiv:2404.06070 (2024).
  * Mutlu  _et al._ [2023] E. Mutlu, A. Panyala, N. Gawande, A. Bagusetty, J. Glabe, J. Kim, K. Kowalski, N. P. Bauman, B. Peng, H. Pathak, _et al._ , Tamm: Tensor algebra for many-body methods, The Journal of Chemical Physics 159 (2023).
  * Kissas  _et al._ [2020] G. Kissas, Y. Yang, E. Hwuang, W. R. Witschey, J. A. Detre, and P. Perdikaris, Machine learning in cardiovascular flows modeling: Predicting arterial blood pressure from non-invasive 4d flow mri data using physics-informed neural networks, Computer Methods in Applied Mechanics and Engineering 358, 112623 (2020).
  * Raissi  _et al._ [2019] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational physics 378, 686 (2019).
  * Raissi  _et al._ [2020] M. Raissi, A. Yazdani, and G. E. Karniadakis, Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations, Science 367, 1026 (2020).
  * Norambuena  _et al._ [2024] A. Norambuena, M. Mattheakis, F. J. González, and R. Coto, Physics-informed neural networks for quantum control, Physical Review Letters 132, 010801 (2024).
  * Kowalski  _et al._ [2024] K. Kowalski, B. Peng, and N. P. Bauman, The accuracies of effective interactions in downfolding coupled-cluster approaches for small-dimensionality active spaces, The Journal of Chemical Physics 160, 224107 (2024).
  * Tancik  _et al._ [2020] M. Tancik, P. Srinivasan, B. Mildenhall, S. Fridovich-Keil, N. Raghavan, U. Singhal, R. Ramamoorthi, J. Barron, and R. Ng, Fourier features let networks learn high frequency functions in low dimensional domains, Advances in neural information processing systems 33, 7537 (2020).
  * Hendrycks and Gimpel [2016] D. Hendrycks and K. Gimpel, Gaussian error linear units (gelus), arXiv preprint arXiv:1606.08415 (2016).
  * Sitzmann  _et al._ [2020] V. Sitzmann, J. Martel, A. Bergman, D. Lindell, and G. Wetzstein, Implicit neural representations with periodic activation functions, Advances in neural information processing systems 33, 7462 (2020).
  * Liang  _et al._ [2024] S. Liang, L. Lyu, C. Wang, and H. Yang, Reproducing activation function for deep learning, Communications in Mathematical Sciences 22, 285 – 314 (2024).