A general framework for active space embedding methods: applications in

quantum computing
Stefano Battaglia,1, a) Max Rossmannek,1, 2, b) Vladimir V. Rybkin,1, c) Ivano Tavernelli,2 and Jürg Hutter1
1)
Department of Chemistry, University of Zurich, Winterthurerstrasse 190, CH-8057 Zürich,
Switzerland
2)
IBM Quantum, IBM Research – Zurich, CH-8803 Rüschlikon, Switzerland
(Dated: 30 April 2024)
We developed a general framework for hybrid quantum-classical computing of molecular and periodic embed-
ding calculations based on an orbital space separation of the fragment and environment degrees of freedom.
arXiv:2404.18737v1 [physics.chem-ph] 29 Apr 2024

We show its potential by presenting a specific implementation of periodic range-separated DFT coupled to
a quantum circuit ansatz, whereby the variational quantum eigensolver and the quantum equation-of-motion
approach are used to obtain the low-lying spectrum of the embedded fragment Hamiltonian. Application of
this scheme to study strongly correlated molecular systems and localized electronic states in materials is show-
cased through the accurate prediction of the optical properties for the neutral oxygen vacancy in magnesium
oxide (MgO). Despite some discrepancies in absorption predictions, the method demonstrates competitive
performance with state-of-the-art ab initio approaches, particularly evidenced by the accurate prediction of
the photoluminescence emission peak.

I. INTRODUCTION which enables the efficient representation of an other-

wise inaccessible computational space. Many algorithms
Over the past decades, the electronic structure simu- have been developed to leverage this representation for
lation of molecular and solid-state systems has assumed the computation of ground and excited-state properties
an increasingly important role in the research and de- of chemical systems3 .However, these algorithms exceed
velopment of new materials. While the exact computa- the quantum computational capabilities of state-of-the-
tion of ground and excited state properties poses an ex- art devices when applied to large problems. On one hand,
ponentially scaling problem through the solution of the this is due to the deep quantum circuits that arise when
Schrödinger equation, many methods have been devel- encoding a fermionic WF in a set of qubits, while ensur-
oped to various degrees of approximation rendering their ing the physical nature of the generated superposition.
implementation feasible. In particular, density functional The long runtime of such quantum circuits exceeds the
theory (DFT) has established itself as a cheap yet effec- decoherence times of the available hardware, resulting
tive mean for the simulation of a wide range of systems in the accumulation of errors which require error cor-
of interest. However, due to its nature, it falls short in rection to be feasible before we can gain any significant
the description of problems that contain strongly cor- results3 .On the other hand, algorithms which trade the
related electrons. In general, the accurate description execution of a single long quantum circuit for a (varia-
of such systems requires methods that consider multiple tional) optimization problem involving many shorter cir-
configuration state functions (CSFs) in order to capture cuits quickly become infeasible; in particular when tar-
the complex entanglement between the electrons. These geting systems that are beyond the current capabilities
methods of greater accuracy come at increasingly higher of classical computational hardware4 .
computational costs approaching the theoretical expo- As a short-term solution to these limitations, combined
nential scaling of the exact solution, therefore limiting with the interest of the computational chemistry and ma-
their applicability to problems of relatively small sizes. terials science communities to leverage this new comput-
In the past decade, significant progress has been made ing platform, we have seen in recent years increasing ef-
in the development of quantum computers, which provide forts in the development of new hybrid quantum-classical
access to a new computational paradigm that promises algorithms. In particular, many embedding methods
to overcome this exponential barrier1,2 . This feat can have been adapted or newly developed to leverage the
be achieved by encoding the exponentially scaling wave quantum resources for the treatment of embedded frag-
function (WF) in a polynomial number of qubits, the ments, where in many cases the computation of a small
fundamental processing units of a quantum computer, region of a larger system is offloaded to the quantum
computer, while the rest is carried out on classical hard-
ware5–13 . These efforts have been very important for
a) Electronic exploring the current practical limitations of noisy de-
mail: stefano.battaglia@chem.uzh.ch; These authors
contributed equally to this work. vices and for benchmarking quantum computational ap-
b) Electronic mail: oss@zurich.ibm.com; These authors contributed proaches for the simulation of chemical and solid-state
equally to this work. systems. Crucially, embedding methods also provide the
c) now at: HQS Quantum Simulations GmbH, Rintheimer Strasse
means to systematically scale the problems of interest
23, DE-76131 Karlsruhe, Germany alongside the development of quantum computing hard-
 2

ware, such that the infrastructure and resources invested tronic Hamiltonian in the Born-Oppenheimer approxi-
now will eventually allow to reach system sizes and an mation. This can be written in atomic units as
accuracy beyond what is currently feasible. X 1X
Adding to these efforts, this work presents a general Ĥ = hpq â†p âq + gpqrs â†p â†r âs âq + V̂nn , (1)
framework for the implementation of active space embed- pq
2 pqrs
ding methods to couple classical and quantum computa-
tional codes. In particular, building up on an earlier work where V̂nn is the Coulomb repulsion between the nuclei,
by some of the authors5 , we extend the range-separated while hpq and gpqrs are one- and two-electron integrals
DFT embedding scheme to allow embedding into peri- given by
odic environments. The new and generally improved im- D E
plementation, coupling the CP2K package14 and Qiskit hpq = ψp (x) ĥ ψq (x) , (2)
Nature15,16 , is capable of simulating both molecular and
periodic systems alike. The communication between the gpqrs = ⟨ψp (x)ψr (x′ ) | ĝ | ψq (x)ψs (x′ )⟩ . (3)
two codes is handled through message passing permitting The variable x (x′ ) is a compound variable for the elec-
future extensions as well as providing a scalable path to- tron coordinates in space, r (r′ ), and its spin degree of
wards quantum-centric supercomputing.
freedom. The operators ĥ and ĝ in Eqs. (2) and (3) ac-
The rest of this paper is structured as follows. In
count for the kinetic energy of the electrons, the electron-
Section II we recap the theory of active space embed-
nuclear attraction, and the electron-electron repulsion,
ding routines and extend the previously published range-
and are defined as
separated DFT embedding to periodic systems. We also
review the most relevant concepts of quantum comput- 1 X −ZP
ing required for the simulation of electronic structure sys- ĥ(r) = − ∇2 + , (4)
2 |r − RP |
P
tems. Section III discusses the details of the integration
between CP2K and Qiskit Nature. We think that this 1
ĝ(r, r′ ) = , (5)
work will present useful and generally applicable guide- |r − r′ |
lines for other hybrid quantum-classical integrations in
where P labels the ion cores, while ZP and RP de-
the future. In Section IV we present and discuss the re-
note the corresponding nuclear charges and nuclear po-
sults obtained by applying the developed methodology
sitions, respectively. The indices p, q, r, s label general
to study the optimal properties of the neutral oxygen va-
one-particle functions (spin-orbitals), whose correspond-
cancy in magnesium oxide. Finally, Section V concludes
ing sums appearing in Eq. (1) run through the entire
this work.
basis set. The operator â†p (âp ) is a creation (annihila-
tion) operator adding (removing) an electron to (from)
spin-orbital ψp (x). In an embedding approach like the
II. THEORY
one presented here, the fragment of interest is defined in
terms of an active space consisting of a handful of ac-
We present a general framework for active space em- tive electrons and active orbitals. All the electrons that
bedding approaches, where a subset of electrons and or- are not part of the active space typically occupy the low-
bitals of a system — the fragment — are embedded in energy states of the system and are called inactive or-
an effective potential generated by the remaining elec- bitals; see Fig. 1 for an example active space selection in
trons of the systems and all the ion cores — the envi- molecules and materials. Once the active space has been
ronment. We do this in the framework of multiconfig- identified, a corresponding embedded fragment Hamilto-
urational range-separated density functional theory (rs- nian can be defined in a manner completely analogous to
DFT)17–23 and the Gaussian and plane waves (GPW) ap- Eq. (1), that is
proach24 , whereby the embedding potential is obtained
from a mean-field method, while the fragment Hamilto- X
emb † 1 X
Ĥ frag = Vuv âu âv + guvxy â†u â†x ây âv , (6)
nian is solved with a correlated wave function ansatz. 2 uvxy
uv
The approach and infrastructure we have developed is
completely general; it can treat both molecular and peri- with the only difference that the sums are limited to the
odic systems, it supports spin-polarized and unpolarized active orbitals, labeled by the indices u, v, x, y, and that
calculations, it describes the core electrons explicitly or the one-electron integrals, hpq , have been replaced by the
through pseuopotentials, and can be combined with both elements of an embedding potential, Vuv emb
. This poten-
classical WF and quantum circuit ansatzes alike. tial accounts for the interactions between the inactive and
active electrons in addition to the contributions from the
one-electron integrals. Notice that until this point this
A. Active space embedding formulation is completely general; we have not yet speci-
fied how to compute the embedding potential. In princi-
To introduce the framework for active space embed- ple, one could define an operator that explicitly accounts
ding methods, we start with the second-quantized elec- for all many-body interactions between the inactive and
 3

are the elements of the one- and two-particle reduced

virtual conduction
subspace band states density matrices (RDMs), D and d, respectively. By
separating the one-particle RDM (1-RDM) into inactive
and active components, D = DI + DA , and factorizing
fragment

the elements of the inactive two-particle RDM (2-RDM)

active into a product of 1-RDMs (see Appendix A.1 and A.2
subspace

fragment
of Rossmannek et al. 5 for a detailed derivation), we can
express the total energy as a sum of inactive and active
parts, E = E I + E A , with
X
EI = (hij + Vijemb )Dij
I
+ Vnn , (10)
environment

ij

environment
inactive
subspace and
valence
band states X 1 X
EA = emb A
Vuv Duv + guvxy dA
uvxy . (11)
uv
2 uvxy

FIG. 1. Example of possible selections of active and inac- In Eqs. (10) and (11), the indices i, j label inactive or-
tive spaces for the water molecule (left) and for a positively bitals and the superscripts I and A on the density ma-
charged boron vacancy in hexagonal boron nitride (right). In trix elements emphasize to which subspace they belong
both cases, only a small number of orbitals is included in the (even though the indices and sums implicitly account for
active space: for water, it is spanned by the HOMO-LUMO that information). At last, notice that the choice of one-
pair, while for boron nitride, by the localized defect orbitals. particle functions is completely general: one can choose
Notice that also the number of active electrons must be spec- localized molecular orbitals in case of molecules, crys-
ified, though, once the active orbitals have been selected, the talline orbitals or Wannier functions in solid-state sys-
number of active electrons is easily identified. tems, or a combination thereof, e.g. to describe point
defects in materials.
active subsystems, though, probably this would result in
a methodology that is computationally as expensive as
solving directly for the entire problem with such an ap- B. Periodic range-separated DFT embedding
proach. Therefore, in practice, the embedding potential
introduces approximations to describe the low-energy de- As a starting point for introducing the periodic range-
grees of freedom and the interactions between active and separated DFT embedding, we can rely on the working
inactive subsystems. For example, one such option would equations of Rossmannek et al. 5 that we will briefly recap
be to use the Hartree–Fock (HF) approximation for the in the following. The first ingredient is the definition of
inactive electrons, such that the active electrons only in- the one-particle embedding potential, with elements
teract with the active ones in a mean-field manner. In
emb I,LR SR
emb
this case, Vuv would simply correspond to the elements Vpq = Fpq + VH,pq [ρI ] + VH,pq
SR
[ρA ] + Vxc,pq
SR
[ρ] , (12)
of the Fock matrix. Similarly, as we will discuss in more
detail in the next section, describing the environment us- where the elements of the inactive long-range Fock oper-
ing DFT translates into an embedding potential similar ator are defined as
to the Kohn–Sham (KS) one. It is important to real- I,LR LR
emb
ize that in general, Vuv always depends on the inactive Fpq = hpq + VH,pq [ρI ] + VHFX,pq
LR
[ρI ] , (13)
electronic degrees of freedom, but possibly also on the
active subsystem, in which case the resulting embedding along with the classical Hartree potential, VH [ρ], the
scheme has to be solved self-consistently (see Fig. 2). exact Hartree-Fock exchange potential, VHFX [ρ], and
To compute the total energy of the system, we can the DFT exchange-correlation potential, Vxc [ρ], evalu-
start from the expectation value of Eq. (1) with respect ated over the indicated electron densities, ρI , ρA and
to the total WF of the system, that is ρ = ρI + ρA (see the Supplementary Information (SI) for
the explicit definition of these operators in a one-particle
D E X 1X basis). The two-electron integrals over the Coulomb op-
E = Ψ Ĥ Ψ = hpq Dpq + gpqrs dpqrs + Vnn ,
pq
2 pqrs erator are split into long-range (LR) and short-range
(7) (SR) components,
where
ĝ = ĝ ω,LR + ĝ ω,SR (14)
Dpq = ⟨Ψ|â†p âq |Ψ⟩ , (8) erf(ω|r − r′ |) erfc(ω|r − r′ |)
= + (15)
dpqrs = ⟨Ψ|â†p â†r âs âq |Ψ⟩ , (9) |r − r′ | |r − r′ |
 4

which give rise to the superscripts LR and SR in Eqs. (12) Ĥ frag should also provide this quantity (more generally,
and (13). The range separation is obtained with the error it should provide the 1-RDM). Crucially, this dependence
function and its complement (as indicated by Eq. (15)), of V emb on ρA implies that our embedding approach re-
where ω is the range-separation (RS) parameter of units quires a self-consistent solution, whereby an updated ac-
a−1
0 . tive density is obtained at each iteration, which is used to
In practice, two issues arise for the direct computation build a refined embedding potential and updated inactive
of the inactive energy and embedding potential accord- energy that account for the feedback of the active subsys-
ing to Eqs. (10) and (12). First, it is computationally tem on the environment degrees of freedom. A scheme
disadvantageous when the inactive subsystem becomes depicting this self-consistent loop is shown in Fig. 2, when
very large. Second, for periodic calculations, the sums discussing the implementation details in Section III.
concerning the electron-electron, electron-nuclear, and Finally, it should be emphasized that Eqs. (18)
nuclear-nuclear interactions are conditionally convergent, and (19) generalize to the molecular and periodic em-
and cannot be easily separated into inactive and active bedding settings, since only the computation of the total
components. Hence, we express the inactive terms in- Fock operator and energy are affected by this change,
directly as the difference between the total system and at least when invoking the Γ point approximation. Fur-
the (localized) active subsystem. We can achieve this thermore, we point out the limiting cases provided by
by defining the inactive 1-RDM and electron density as the RS scheme: they allow us to recover the common
DI = D − DA and ρI = ρ − ρA , respectively. Reformu- HF embedding scheme as ω approaches infinity as well
lating Eq. (13) by replacing ρI = ρ − ρA , yields as KS DFT as ω approaches zero. This can be seen in
Eqs. (18) and (19), where the standard KS case is ev-
I,LR tot SR SR
Fpq = Fpq − VH,pq [ρ] − Vxc,pq [ρ] ident as all LR terms simply disappear. The common
LR
− VH,pq [ρA ] − VHFX,pq
LR
[ρA ] , (16) HF embedding is also evident from Eq. (17), in which
the only DFT-specific term for the exchange-correlation
SR
where the total rsDFT Fock operator is defined as interaction, Vxc , vanishes.

tot LR SR
Fpq = hpq + VH,pq [ρ] + VHFX,pq [ρ] + Vxc,pq [ρ] . (17)
C. Quantum computing
Inserting the same relation for the inactive electron den-
sity as well as Eq. (16) into Eq. (12) results in We leverage quantum computing for finding the ground
emb tot LR and excited state solutions of the embedded fragment
Vpq = Fpq − VH,pq [ρA ] − VHFX,pq
LR
[ρA ] . (18) Hamiltonian (cf.Eq. (6)). We do so through the means
of the newly developed integration between CP2K14 and
We can proceed analogously for the expression of the
Qiskit Nature15,16 which we discuss in more detail in Sec-
inactive energy, obtaining
tion III. In recent years, quantum computing has made
X significant progress in emerging as a new computational
E I = E tot − tot A
Fuv LR A
Duv + EH LR
[ρ ] + EHFX [ρA ] . (19) paradigm that promises great advances for the simulation
uv
of chemistry and material science problems3 .In particu-
The active energy component that is needed to compute lar for the latter, the ability to treat localized fragments
the total energy, E = E I + E A , simply corresponds to embedded into periodic systems using quantum comput-
the ground state of the fragment Hamiltonian, Eq. (6). ing platforms is particularly appealing. While we rely
Owing to the similar structure of Eqs. (1) and (6), es- on state-of-the-art hybrid quantum-classical algorithms
sentially any electronic structure method can be used in such as the variational quantum eigensolver (VQE)4 and
combination with our embedding scheme. In practice, quantum equation of motion (qEOM)27 , the presented
because the space spanned by the active orbitals and elec- embedding framework is not coupled to these choices and
trons is relatively small, exact diagonalization or a good can leverage any advancements in the field of quantum
approximation thereof is the method of choice. Electron- computing that are yet to come. This also holds for the
ically excited states can also be targeted by the embed- modular integration of the two computational codes as
ding method, either by directly calculating the spectrum we will show later. In this section, we briefly review the
of Ĥ frag or by linear response. However, one has to be theoretical foundations of the quantum computational
careful that the inactive subspace is normally optimized tools used throughout this work.
for the ground state, unless some form of state-averaging
or orbital optimization similar to classical multiconfig-
urational quantum chemical methods is introduced25,26 . 1. The fermion-to-qubit mapping
As will be discussed in the next section, we have used
quantum circuit ansatzes to obtain the ground and ex- In order to simulate a fermionic system on a quan-
cited states energies of Eq. (6). Owing to the dependence tum computer one must first map the second-quantized
of the embedding potential to the active space electron Hamiltonian (cf. Eq. (6)) into a form that the quantum
density, ρA , the method chosen to get the spectrum of computer can work with. Since the fundamental opera-
 5

tional unit of a quantum computer is a two-level system, sampling of many shorter duration circuits. The funda-
a so-called qubit, the translation routines are referred mental principle of the VQE is based on the expectation
to as fermion-to-qubit mappings. Many such mappings value computation of an observable, Ô, with respect to
exist28–33 , the most common of which is arguably the some reference state, Ψ, as
Jordan-Wigner mapping28 . It maps the fermionic cre-
ation, â†p , and annihilation, âp , operators acting on spin ⟨Ψ|Ô|Ψ⟩
⟨Ô⟩ = . (23)
orbital, p, to a tensor product of identities, I, and Pauli ⟨Ψ|Ψ⟩
matrices, {σ̂pX , σ̂pY , σ̂pZ }, which correspond to the princi-
pal single-qubit rotations along the Cartesian axes. The Finding the ground state of a Hamiltonian, Ĥ, amounts
mapping can be written as to using a parameterized ansatz for the WF, Ψ(θ), and
variationally optimizing the parameters, θ, with respect
p−1
! M
!
O O to the expectation value, ⟨Ĥ⟩. In recent years, many vari-
âp → σ̂qZ σ̂p− Iq , (20) ants of the VQE algorithm have been developed (see for
q=1 q=p+1
instance34–37 ), which oftentimes leverage specific struc-
p−1
! M
!
O O tures found in particular ansatzes for the WF. Since we
â†p → σ̂qZ σ̂p+ Iq , (21) are only simulating the execution of quantum circuits
q=1 q=p+1 on classical computers, we do not require improved cir-
cuit depths or other benefits brought about by these al-
where σ̂ ± = (σ̂ X ± iσ̂ Y )/2 are the so-called ladder op- gorithm variants. Thus, we employ the unaltered VQE
erators. This mapping is the most straightforward since algorithm.
it encodes the occupation of one spin orbital into the
occupation of a single qubit. Therefore, it will encode
a fermionic Hamiltonian acting on M spin orbitals into 3. The wave function ansatz
a qubit Hamiltonian acting on M qubits. The anti-
commutation relations of the fermionic operators are pre- Choosing an ansatz for the VQE algorithm is a diffi-
served by the chain of σ̂ Z rotations on lower-indexed cult task. On one hand, it has to be expressive enough to
qubits. The resulting qubit Hamiltonian is a weighted contain the true ground state in its parameterized sub-
sum of the form1 space of the entire Hilbert space. On the other hand,
X it should be limited in its circuit depth implementation
Ĥ qu = cp P̂p , (22) to ensure that it can be executed on near-term quantum
p
computers. Both of these properties can be achieved by
means of hardware efficient ansatzes (HEAs), which can
where each Pauli string, P̂p , is a tensor product of identi- even be tailored to respect limited qubit connectivity of
ties and Pauli matrices. Crucially, the number of unique the quantum computing hardware. However, optimiza-
terms in Ĥ qu scales as O(M 4 ), just like the number of tion of such ansatzes can pose to be challenging due to
two-body interactions in the original second-quantized the large number of variational parameters38 . Moreover,
Hamiltonian, Eq. (6). chemical problems, such as the ones we are interested in
Another very common mapping is the parity map- here, are constrained to the Fock space that often ex-
ping30 which can be thought of as the dual to the Jordan- hibit additional symmetries, further restricting the size
Wigner mapping. It encodes the parity information lo- of the subspace containing the true physical ground state.
cally on the qubits and delocalizes the occupation infor- Therefore, chemistry-inspired ansatzes have been devel-
mation of the fermionic spin orbitals. This has an added oped which are designed to explore only this physical
benefit that for particle-number conserving Hamiltoni- subspace, a famous example of which is the unitary cou-
ans, two qubits carry redundant information and may be pled cluster (UCC) ansatz39 . Its major drawback is the
omitted (which is referred to as two-qubit reduction 30 ). significant circuit depth overhead associated with the im-
This is the mapping which was used for all simulations in plementation of these constraints. Nonetheless, it is still
this work, but any other fermion-to-qubit mapping could orders of magnitude cheaper than an implementation of
be used in its place. the QPE algorithm39 .

2. The variational quantum eigensolver 4. Computing excited states

The VQE4 is a hybrid quantum-classical algorithm to Many different (hybrid) quantum algorithms exist for
find the ground state of any Hamiltonian. It has gained the computation of excited states 27,40,41 ; in this work
a widespread interest as an alternative to quantum phase we rely on the qEOM method27 . It has the advantage
estimation (QPE) since it is more amenable to near-term that, once a ground state solution has been found, one
quantum computers. It does so by replacing the execu- only needs to perform additional measurements on this
tion of a single long-running quantum circuit with the optimized WF. Other algorithms, however, may require
 6

an entire new optimization procedure to be completed CP2K Qiskit

for each targeted excited state40 . quantum
The fundamental idea of qEOM relies on a classical SCF user preparation
computer solving the equations of motion while a quan- Socket
tum computer is used for the measurement of the matrix active space computation
elements that go into the system of equations. When in- connect handshake connect
cluding only single and double excitations in the opera-
tor basis, the measurement cost of these matrix elements initialization LR fetch
scales just like the measurement cost of the Hamiltonian guvxy (15)
i = 0 transferred
expectation value with its number of terms O(M 4 ). This data &

loop
can be seen from the expectation values that give rise to update
emb A,(i)
Vuv [ρ ] (18) construct
the excitation energies that we are after i += 1 Ĥ frag (6)
E I [ρA,(i) ] (19)
No
⟨0|[Ôn , Ĥ, Ôn† ]|0⟩ Duv
A,(i+1)
E0n = , (24) done? solve (11)
⟨0|[Ôn , Ôn† ]|0⟩ E A,(i+1)

Yes
where |0⟩ denotes the ground state and Ôn† = |n⟩ ⟨0| is shutdown close shutdown
the excitation operator from the ground state to the n-
th excited state (and Ôn is the matching de-excitation
operator). For more details we refer the interested reader FIG. 2. Workflow diagram depicting the interaction of CP2K
to the original paper by Ollitrault et al. 27 . and Qiskit Nature. The user configures the two classical pro-
Once again, our choice of using qEOM is not constrain- cesses and the socket for the IPC. Each process then follows
ing the future applicability of any other (hybrid) quan- the computational steps (rectangular boxes) outlined inside
tum algorithm to be used for the computation of excited of their respective frames. The data that gets computed
states. and transferred is indicated by the rounded boxes. Num-
bers in parentheses refer to the respective equations in this
manuscript. The self-consistent embedding requires a loop
which is highlighted by the gray box. This loop is terminated
III. IMPLEMENTATION based on the decision (diamond shape) taken by the CP2K
process.
In this section, we present the implementation details
of the integration of CP2K14 and Qiskit Nature15,16 . The
developments of this work have been released as part of in order to exchange data between the two codes. Partic-
CP2K v2024.1, Qiskit Nature v0.7.0, as well as a sep- ularly, for this initial implementation the messages and
arate module handling more specific parts of the inte- data are sent over a socket file. This is inspired by a sim-
gration called qiskit-nature-cp2k42 . In the first part, ilar architecture used by the i–Pi project44 . Additional
we discuss the technical aspects and challenges. Later, technical details are available in the SI.
we review the future scalability and extensibility of this A socket is an application programming interface
design. (API) used for inter-process communication (IPC). Using
this protocol, it is possible for the communicating pro-
cesses to run on the same physical machine or different
A. Technical realization ones connected via the internet. The calculation proceeds
identically in both scenarios with the only difference be-
Interfacing CP2K with Qiskit Nature for the imple- ing the latency of the communication. However, this is
mentation of an iterative embedding scheme poses a num- not of a concern to us, since the rate-limiting factor of
ber of challenges. While CP2K is primarily written in the communication is (in any case) the bandwidth of the
Fortran and provides the means to efficiently run highly connection to the quantum hardware.
parallelized simulations in a variety of computational se- Fig. 2 summarizes the computational workflow of our
tups, Qiskit Nature (and the underlying Qiskit software integration between CP2K and Qiskit Nature. The di-
development kit (SDK)) is mostly developed in Python agram depicts a user who has to configure three parts
and has not yet43 reached a computational maturity com- of their calculation; the CP2K and Qiskit (Nature) pro-
parable to CP2K. When Qiskit Nature was coupled to cesses depicted on the left and right, respectively, as well
other Python-based computational programs in the past, as the socket itself via which the messages are passed
they could easily share the same Python runtime execu- between the two codes. Both computational codes will
tion environment and, thus, share all data directly in start in parallel. While CP2K starts out by finding a
memory5,13 . For the integration discussed here, this was low-level solution to the entire system (SCF), Qiskit can
not possible in such a straightforward manner. Instead, use this time to perform certain preparational tasks that
our implementation relies on a message passing protocol are unique to the execution on quantum computing hard-
 7

ware and do not require problem-specific data. For both and design as scalable. Furthermore, future improve-
programs, the user has full flexibility to leverage their ments to aid in the transfer of data to the QPU that will
respective capabilities during these initial steps. Upon be implemented into the quantum stack will be accessi-
completion of their respective steps, both codes will syn- ble directly to the end-users of our integration because it
chronize by performing a handshake through the socket. only serves as a middle-man between the two codes.
If either process reaches this point before the other, it
awaits the other one. CP2K reaches this point inside its
active space module which was released as part of CP2K IV. RESULTS AND DISCUSSION
v2024.1. The input to this module configures the active
fragment to be embedded into its environment and com-
putes the one- and two-body terms of the active space To test our implementation we have studied the opti-
Hamiltonian (i.e. the fragment Hamiltonian). It al- cal properties of the F 0 -center (neutral oxygen vacancy)
lows both single-shot and iterative embedding routines in magnesium oxide, whose nature remains unclear de-
to be performed. In Fig. 2 we have depicted the rsDFT spite the many experimental45–48 and computational49–56
embedding protocol presented in Section II B. However, studies carried out in the past decades. This data will
some of the components do not change throughout the allow us to compare the accuracy of the periodic rsDFT
course of the self-consistent embedding. These can be embedding with respect to both experimental spectra as
pre-computed only once during the initialization proce- well as state-of-the-art ab initio methods.
dure. A key example are the LR-electron-repulsion inte-
grals (ERIs) (cf. Eqs. (3) and (15)). Therefore, these only
have to be transferred to the Qiskit process once. Com- A. Computational details
ponents which depend on the active density of the cur-
rent iteration, ρA,(i) , have to be updated and exchanged
We have considered four different supercell sizes: the
during every iteration of the loop indicated by the gray
2x2x2, 3x3x3 and 4x4x4 supercells constructed from the
frame in Fig. 2. During every such iteration, Qiskit Na-
primitive unit cell, and the 2x2x2 supercell constructed
ture constructs the Hamiltonian of the active fragment
from the conventional unit cell; see Fig. 3. For each sys-
(cf. Eq. (6)) using the LR-ERIs and embedding potential,
emb A,(i) tem, we have first optimized the supercell parameters and
Vuv [ρ ], (cf. Eq. (18)). It then proceeds with finding
geometry of the pristine system with spin-unpolarized
the ground-state solution to this Hamiltonian using the
KS-DFT and the Perdew, Burke, and Ernzerhof (PBE)
quantum circuit ansatz specified by the user. Upon com-
functional57 , together with the geometrical response va-
pletion, it will return the active energy, E A,(i+1) , and
A,(i+1) lence triple-ζ basis set58 and correlation-consistent polar-
active 1-RDM, Duv , to the CP2K process. CP2K ization functions (ccGRB-T basis set in CP2K). The core
will then perform a convergence check based on the total electrons were described by the Goedecker-Teter-Hutter
energy, E = E I + E A,(i+1) . If convergence has not been pseudopotential optimized for the PBE functional59–61 .
reached, CP2K and Qiskit will return back to their re- The DFT calculations were performed within the GPW
spective steps of the embedding protocol to proceed with approach24,62 , using plane-wave absolute and relative
another iteration. If this check succeeds, both processes cutoffs of 1000 and 50 Rydberg, respectively, and a 4-
will be signalled to terminate. During their shutdown layer grid for the numerical integration.
procedures, both processes may perform additional post- To create the vacancy, we removed a single oxygen
processing steps. For example, Qiskit Nature may com- atom from each optimized pristine supercell and relaxed
pute additional properties using the final ground-state again the geometry of all systems with the same set-
WF including the computation of excited state energies. tings as just discussed, while keeping the cell parame-
ters fixed. Note that the ground state of the F 0 -center
in MgO is closed-shell, and therefore spin-unpolarized
B. Future scalability DFT describes it well. To perform the rsDFT embed-
ding calculations we reduced the basis set to a double-ζ
At this point one might wonder how scalable this de- plus polarization (ccGRB-D) for all atoms but the six
sign is for the future. Indeed, the transfer of the two-body magnesium ones surrounding the vacancy, for which we
integrals is the limiting factor here. If CP2K and Qiskit kept the triple-ζ one. In addition, we have also added
Nature were able to leverage a shared memory, this would the triple-ζ basis functions of oxygen at the vacancy site,
alleviate the need for transfer completely. However, only caculated as the center of mass of the six coordinated
up to the point where the mapped qubit Hamiltonian magnesium atoms. Furthermore, we changed the pseu-
needs to be transferred to the QPU. Until we have a di- dopotential to the one optimized for hybrid functionals.
rect high-bandwidth connection between the CPU and We used the local density approximation (LDA) func-
QPU, this transfer of data (also known as the data load- tional in its SR form63–65 in combination with a trun-
ing problem in quantum computing) will remain the rate- cated Coulomb potential66 for the LR HF component of
limiting factor. Therefore, within the scope of this more the functional, using a truncation radius of 4.25 Å. The
general problem, we deem the current implementation RS parameter was set to 0.14 a−1 0 , which was optimally
 8

Primitive (2x2x2) Primitive (3x3x3) Primitive (4x4x4) Conventional (2x2x2)

15 atoms 53 atoms 127 atoms 63 atoms

FIG. 3. The four supercells used in the calculations. Magnesium atoms are in green, oxygen atoms are in red, the oxygen
vacancy is colored in black.

proach27 . The energies for the photoluminescence spec-

trum were computed at the relaxed geometry of the 3T1u
excited state, optimized using time-dependent density
t1u t1u t1u
functional theory (TDDFT) within the Tamm–Dancoff
CB
8 approximation68,69 on top of the 1A1g ground state. The
Orbital energy (eV)

CBM
other computational settings were the same as those used
6 for the ground state geometry relaxation. The AS embed-
a1g ding calculations for the singlet and triplet states were
4 based on spin-unpolarized and spin-polarized DFT, re-
spectively, and were converged to an energy change of
2 1 × 10−8 E h . At last, notice that the q-UCCSD ansatz
for 2 electrons is equivalent to an exact diagonalization,
0
which we confirm by performing all calculations using
VB the classical full configuration interaction (FCI) solver of
PySCF70 instead of VQE plus qEOM. Additional infor-
FIG. 4. Band diagram of the F 0 -center in magnesium ox- mation on the computational details is provided in the
ide. The ground state is a closed-shell A1g singlet, with two SI.
electrons occupying a defect orbital within the gap. Singlet
and triplet excitons occur when an electron from the mid-gap
orbital is excited either to the fully delocalized CBM state, or B. Absorption and photoluminescence spectra
to one of the three t1u defect orbitals within the conduction
band. The 4 defect orbitals, along with the CBM one are
included in the active space. A neutral oxygen vacancy in MgO, typically called F 0 -
center or V0O vacancy, introduces a 1 s-type localized de-
fect orbital at mid-gap that is doubly occupied in the
1
tuned by matching the bandgap of the largest pristine su- A1g electronic ground state. We shall call this orbital
percell to the experimental bandgap value of 7.77 eV67 . the mid-gap orbital. Three more defect-localized degen-
The embedding calculations were performed by includ- erate one-particle states of t1u symmetry (p-like orbitals)
ing 2 electrons and 5 orbitals (10 spin-orbitals in total) in appear within the conduction band. These three orbitals
the active space (AS), which are shown within the band are energetically slightly above the CBM, which corre-
structure diagram in Fig. 4. Four of the five active or- sponds to the delocalized s-band. The four orbitals local-
bitals are localized at the vacancy and are labeled accord- ized at the oxygen vacancy and the CBM one are shown
ing to the (localized) octahedral symmetry around the in Fig. 4 and are believed to be responsible for the optical
defect (Oh point group). The remaining orbital included properties of defective MgO; for this reason we included
corresponds to the conduction band minimum (CBM) them in the active space of our embedding calculation,
and is a fully delocalized conduction s-band, which we along with the 2 electrons occupying the mid-gap defect
label as CBM in the following. The ground state energy state.
of the embedded fragment Hamiltonian was obtained us- Absorption. Experimentally, it is well established
ing the VQE algorithm with a quantum UCC singles and that the absorption peak of the F 0 -center is at 5.03 eV,
doubles (q-UCCSD) ansatz, and the resulting quantum which is extremely close to that of the F + -center (that
circuits were simulated without the addition of artificial is, a positively charged oxygen vacancy) at 4.96 eV45,46 .
noise using Qiskit (version 0.45). We identified the vertical excitation energies correspond-
The absorption and emission spectra were obtained by ing to the transition of one electron from the mid-gap
calculating the lowest ten excitonic states of singlet and orbital to either the CBM or the t1u orbitals, denoted
triplet spin symmetry, respectively, using the qEOM ap- as 1A1g → 1CBM and 1A1g → 1T1u , respectively. In Fig. 5
 9

fact that they likely target the CBM state rather than
12 the defect one. This may suggest that the transition to
10 the 1T1u obtained with those method is possibly over-
Energy (eV)

8 estimating the experimental absorption maximum. Our

6 calculations consistently predict a blue-shifted absorp-
4 tion by ∼0.5 eV compared to the embedded-BSE@DDH
of Vorwerk and Galli 55 . One possible reason could be
2 1
A1g → 1CBM 3
CBM→ 1A1g
1 the small active space of 2 electrons in 5 orbitals for
0 A1g → 1T1u 3
T1u → 1A1g
both the ground state VQE-q-UCCSD calculation and
−2 the subsequent qEOM one, which does not allow suffi-
0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07
cient relaxation of the t1u orbitals, as these are empty
−1
Natoms in the ground state calculation and thus not optimized.
Possible solutions for this problem are discussed in Sec-
FIG. 5. Vertical absorption and emission energies as a func- tion V. Interestingly, density-matrix embedding theory56
tion of the (inverse) number of atoms in the supercell. Filled predicts the absorption to the CBM state to be higher
green squares and blue circles correspond to calculated singlet than the one to the localized state, perhaps owing to
absorption energies, while filled upside orange and downside the missing extrapolation to the TDL. The only method
yellow triangles correspond to calculated triplet emission ener-
that underestimates the absorption is TDDFT, with an
gies. Dashed lines are linear extrapolation curves to the TDL,
whose value is marked with a corresponding empty symbol.
absorption peak centered at 4.85 eV52 .
Photoluminescence. The photoluminescence (PL) of
defective MgO is considerably more complicated than
−1 1
cell Natoms A1g → 1CBM 1A1g → 1T1u the absorption, and several interpretations have been
p 2x2x2 0.067 10.27 12.59 brought forward throughout the years. The main source
p 3x3x3 0.019 6.24 7.74 of ambiguity is the vicinity of the absorption peaks of the
c 2x2x2 0.016 5.59 7.06 F 0 , F + and F 2+ centers, which are all likely to be ex-
p 4x4x4 0.008 5.52 6.83 cited by incoming irradiation at 5 eV, significantly com-
→∞ → 0.0 4.60 5.78 plicating the assignment of the emission bands to the
correct point defect and electronic state. Experimentally,
TABLE I. Optical singlet absorption energies (in eV) of the there are two very distinct peaks visible in the PL spec-
F 0 -center in MgO. The bottom row contains the values ex- trum, one at 2.3 eV and one at 3.2 eV47,48 . The former
trapolated to the TDL. has been associated to an emission from the F 0 -center,
while the latter to an emission from the F + -center47,48 .
The initial explanation for the long-lived nature of the
we show in blue and green the energies obtained for these 2.3 eV band was based on temperature-dependent exper-
states as a function of the number of atoms in the super- iments carried out on samples prepared in different ways
cell. Notice that the energies calculated with our em- and containing different concentrations of F 0 -centers and
bedding method are obtained directly from many-body hydrogen impurities. The proposed mechanism involves
wave functions, thus they account for electron-hole inter- the escape of an electron in the conduction band upon
actions and the exciton binding energy. We performed excitation of the F 0 -center, leaving behind a positively
an extrapolation to the TDL assuming a Natoms −1
conver- charged oxygen vacancy, an F + -center. Hydrogen im-
gence to correct for finite-size effects, which is shown as purities in the sample then act as traps for the mobile
dashed lines in Fig. 5. Our best estimates are reported electrons, which may be thermally released back at a
at the bottom of Table I, along with all the energies later time into the conduction band. Two processes can
obtained for the different supercell sizes. The transi- happen with the released electrons: they may encounter
tion to the 1CBM state is predicted to be 4.60 eV, while F + -centers left behind after the absorption process, in
for the (triply degenerate) one to the 1T1u state we get which case they recombine in an excited 1T1u state of the
5.78 eV. The experimental absorption peak at 5.03 eV F 0 -center that quickly emits light, or they may be recap-
corresponds to the 1A1g → 1T1u , as the excitation to the tured by H – traps, slowing down the overall emission pro-
CBM is dipole-forbidden; we therefore overestimate the cess. An alternative, more straightforward interpretation
experimental value by 0.75 eV. This result can be com- is that the emission band is simply due to triplet phos-
pared to a number of other computational studies on phorescence from a localized 3T1u state at the defect,
the F 0 -center in magnesium oxide, which have been per- accessed via inter-system crossing from the excited 1T1u
formed with a large variety of methods from both the state. This second interpretation is the most accepted ex-
solid-state physics and quantum chemistry communities. planation in recent computational studies, corroborated
From Table II we can see that most of the quantum chem- by calculations based on advanced ab initio methodolo-
istry methods overestimate the absorption to the local- gies52,53,55,56 (in contrast to the first interpretation based
ized 1T1u state, while FN-DMC51 and G0 W0 -BSE50 re- on older semi-empirical methods and experiments47,48 ).
produce the experimental absorption peak, despite the In light of this analysis, we also investigate the second
 10

method structure 1A1g → 1CBM 1A1g → 1T1u emission

qEOM-srLDA (this work) periodic 4.60 5.78 2.44
NEVPT2-DMET@ROHF56 periodic 5.67a 5.24 2.89
embedded-BSE@DDH55 periodic 4.13 5.23 2.93b
EOM-CCSD (GTOs)54 periodic - 5.31 -
EOM-CCSD (PWs)53 periodic - 5.28 3.66
TDDFT@PBE052 periodic - 4.85 2.90
FN-DMC@PBE51 periodic 5.00 - 3.80b
G0 W0 @LDA0-BSE50 periodic 4.95 - 3.40b
CASPT2(2,2)49 cluster - 5.44 4.09
experiment - dark 5.0345,46 2.3047,48
a Value corresponding to the p 4x4x4 system.
b Value obtained for the emission from a different state than 3T1u .

TABLE II. Comparison of predicted optical absorption and photoluminescence energies (in eV) obtained with different com-
putational methods. Different computational studies computed the emission energies from different states, see the footnotes
and the main text for more information.

cell −1
Natoms 3
T1u → 1A1g 3
CBM→ 1A1g band to the F 0 -center, though, originating from a differ-
ent mechanism than the originally proposed one.
p 2x2x2 0.067 6.59 10.05
p 3x3x3 0.019 3.54 5.74
c 2x2x2 0.016 3.38 5.38
p 4x4x4 0.008 3.02 5.02 V. CONCLUSIONS
→∞ → 0.0 2.44 4.14
We developed a general framework for hybrid
TABLE III. Photoluminescence energies (in eV) of the F 0 - quantum-classical molecular and periodic embedding cal-
center in MgO. The bottom row contains the values extrapo- culations based on an orbital space separation of the sys-
lated to the TDL. tem into fragment and environment. This framework has
been implemented in the CP2K package, leveraging many
of its functionalities and taking advantage of its high par-
pathway as the leading emission process and compute the allel efficiency. The modular nature of the implementa-
photoluminescence from the relaxed 3T1u structure and tion allows to easily develop several types of embedding
state. Our results are shown by the orange and yellow schemes and to couple different solvers to obtain ground
curves in Fig. 5 and listed in Table III. The predicted and excited states energies and properties of the embed-
value for the emission from the localized triplet state ded fragment Hamiltonian. It supports both classical
is 2.44 eV and is in very good agreement with the ex- wave function and quantum circuits ansatzes, and the
perimental value of 2.3 eV. In particular, we are much communication between CP2K and the solver is handled
closer to this value than other methodologies, which con- by sockets, which seamlessly integrates within current
sistently predict higher energies for this emission band supercomputing facilities, but is also ready for a more
as shown in the last column of Table II. Nevertheless, quantum-centric high-performance computing vision.
one has to be careful when comparing different theo- To demonstrate the potential of the new framework in
retical works, since these focused on different states or practice, we have implemented a range-separated DFT
mechanisms. The work by Vorwerk and Galli 55 reported embedding scheme that enables the study of both finite
2.93 eV as the PL from the F + -center, hence to be com- and extended systems. This approach is essentially an
pared with the experimental value of 3.2 eV. The ear- extension of multiconfigurational range-separated DFT
lier work based on FN-DMC51 and G0 W0 -BSE50 re- to periodic boundary conditions and relies on the range-
ported the emission from the singlet state (whether the separation of the two-electron integrals in long- and
CBM or T1u state is unclear), and have concluded that short-range components. Within this approach, any cor-
the assignment from the experimental studies should be related wave function method can in principle be coupled
re-evaluated, with the 3.2 eV peak assigned to the F 0 - with DFT in a self-consistent scheme, whereby the for-
center rather than the F + center. The methodologies mer is used to obtain the spectrum of a fragment Hamil-
based on quantum chemistry methods, that is NEVPT2- tonian, and the latter to construct an embedding poten-
DMET, EOM-CCSD, TDDFT and CASPT2, calculate tial generated by the environment degrees of freedom. In
the transition 3T1u → 1A1g like us and compare their re- particular, as part of this work, we have implemented
sults against the 2.3 eV band49,52,53,56 . Here we do the an interface to Qiskit Nature42 that allows to map the
same and we get the best theoretical result so far, which fermionic fragment Hamiltonian to a qubit Hamiltonian,
corroborates the experimental assignment of the 2.3 eV whose ground and excited states can be obtained with
 11

the quantum algorithm of choice. can be interfaced to other classical active space solvers,
The developed rsDFT embedding scheme has a wide such as those based on selected configuration interaction,
application scope, allowing the investigation of both which would also allow to study larger fragment Hamil-
strongly correlated molecular systems, as well as localized tonians solely on classical hardware.
electronic states in materials, such as those arising from
vacancies and impurities. To this end, we have demon-
strated its accuracy and applicability by studying the op- ACKNOWLEDGEMENTS
tical properties of the neutral oxygen vacancy in MgO,
whereby both defect-localized and delocalized states have The authors thank Valery Weber for insightful discus-
been treated on equal footing, and the low-lying spec- sions and guidance during the This research was sup-
trum of the embedded fragment Hamiltonian has been ported by the NCCR MARVEL, a National Centre of
calculated by VQE and qEOM, in combination with the Competence in Research, funded by the Swiss National
q-UCCSD ansatz. Our calculations for the absorption Science Foundation (grant number 205602). This work
spectrum predict a peak at 5.78 eV, a value that over- has been also supported by the Swiss National Sci-
estimates the experimental result by 0.75 eV, but which ence Foundation in the form of Ambizione grant No.
lies in the same ballpark as other sophisticated compu- PZ00P2 174227 (VVR). early stages of the code develop-
tational approaches. On the other hand, the predicted ment. IBM, the IBM logo, and ibm.com are trademarks
PL emission of 2.44 eV from the 3T1u state almost per- of International Business Machines Corp., registered in
fectly matches with the experimentally measured signal many jurisdictions worldwide. Other product and ser-
at 2.3 eV, and provides new insight on a system that has vice names might be trademarks of IBM or other compa-
eluded state-of-the-art ab initio approaches for the last nies. The current list of IBM trademarks is available at
decade. While the accuracy of the method for the ab- https://www.ibm.com/legal/copytrade.
sorption leaves room for improvement, and the excellent
agreement for the emission is certainly helped by favor-
able error compensation, the present study shows that SUPPLEMENTARY INFORMATION
current hybrid quantum-classical algorithms can compete
with classical state-of-the-art ab initio methodologies for The Supplementary Information contains the deriva-
problems beyond simple model systems. tion of the HF approach within the GPW formalism,
Many possible future directions are envisioned based showing the expression of the potentials expressed in the
on this work. On one hand, the periodic rsDFT embed- molecular orbital basis. It also contains the technical de-
ding scheme can be extended in several ways. For in- tails about the message passing interface used to couple
stance, introducing orbital optimization would allow to CP2K and Qiskit nature, and additional computational
incorporate the feedback from the correlated active space details regarding the calculations performed. All input
wave function onto the inactive long-range component; and output files, structure files, tabulated raw data and
this would be particularly important for accounting the scripts to perform the simulations are available on the
changes in the environment when targeting states other Materials Cloud platform at 10.24435/materialscloud:47-
than the ground state. Furthermore, it would make the 6g.
embedding scheme truly variational, significantly simpli- 1 P. K. Barkoutsos, J. F. Gonthier, I. Sokolov, N. Moll, G. Salis,
fying the calculation of analytical forces. State-averaging A. Fuhrer, M. Ganzhorn, D. J. Egger, M. Troyer, A. Mezzacapo,
would also be a useful extension, allowing a more bal- S. Filipp, and I. Tavernelli, Phys. Rev. A 98, 022322 (2018).
2 M. Motta and J. E. Rice, WIREs Computational Molecular Sci-
anced description of several states simultaneously. This is
fundamental in cases where near-degeneracies are promi- ence 12, e1580 (2022).
3 Y. Alexeev, M. Amsler, P. Baity, M. A. Barroca, S. Bassini,
nent, such as in molecules and materials containing open- T. Battelle, D. Camps, D. Casanova, Y. jai Choi, F. T. Chong,
shell transition metals. While optimally-tuned short- C. Chung, C. Codella, A. D. Corcoles, J. Cruise, A. D. Meglio,
range LDA performed well in this study, implementing J. Dubois, I. Duran, T. Eckl, S. Economou, S. Eidenbenz,
more SR functionals will provide alternatives to cases B. Elmegreen, C. Fare, I. Faro, C. S. Fernández, R. N. B. Fer-
reira, K. Fuji, B. Fuller, L. Gagliardi, G. Galli, J. R. Glick,
where LDA is not sufficiently accurate.
I. Gobbi, P. Gokhale, S. P. de la Gonzalez, J. Greiner, B. Gropp,
Future directions that are not strictly tied to the rs- M. Grossi, E. Gull, B. Healy, B. Huang, T. S. Humble, N. Ito,
DFT embedding scheme are also envisaged. For in- A. F. Izmaylov, A. Javadi-Abhari, D. Jennewein, S. Jha, L. Jiang,
stance, one possibility is to implement orbital localization B. Jones, W. A. de Jong, P. Jurcevic, W. Kirby, S. Kister,
M. Kitagawa, J. Klassen, K. Klymko, K. Koh, M. Kondo, D. M.
schemes based on maximally localized Wannier functions Kurkcuoglu, K. Kurowski, T. Laino, R. Landfield, M. Leininger,
and pair natural orbitals, which would allow to study V. Leyton-Ortega, A. Li, M. Lin, J. Liu, N. Lorente, A. Luckow,
pristine solid-state materials, where localized states do S. Martiel, F. Martin-Fernandez, M. Martonosi, C. Marvinney,
not arise naturally due to symmetry-breaking of the su- A. C. Medina, D. Merten, A. Mezzacapo, K. Michielsen, A. Mi-
percell. Owing to the local nature of electron correlation, tra, T. Mittal, K. Moon, J. Moore, M. Motta, Y.-H. Na, Y. Nam,
P. Narang, Y.-y. Ohnishi, D. Ottaviani, M. Otten, S. Pakin,
orbital localization could also simplify the construction of V. R. Pascuzzi, E. Penault, T. Piontek, J. Pitera, P. Rall, G. S.
hardware efficient ansatzes, thereby potentially increas- Ravi, N. Robertson, M. Rossi, P. Rydlichowski, H. Ryu, G. Sam-
ing the maximum size of the active space. Finally, CP2K sonidze, M. Sato, N. Saurabh, V. Sharma, K. Sharma, S. Shin,
 12

G. Slessman, M. Steiner, I. Sitdikov, I.-S. Suh, E. Switzer, 31 Z. Jiang, A. Kalev, W. Mruczkiewicz, and H. Neven, Quantum
W. Tang, J. Thompson, S. Todo, M. Tran, D. Trenev, C. Trott, 4, 276 (2020).
H.-H. Tseng, E. Tureci, D. G. Valinas, S. Vallecorsa, C. Wever, 32 Y. Shee, P.-K. Tsai, C.-L. Hong, H.-C. Cheng, and H.-S. Goan,

K. Wojciechowski, X. Wu, S. Yoo, N. Yoshioka, V. W.-z. Yu, Physical Review Research 4, 023154 (2022).
S. Yunoki, S. Zhuk, and D. Zubarev, arXiv:2312.09733v1 (2023). 33 A. Miller, Z. Zimborás, S. Knecht, S. Maniscalco, and G. Garcı́a-
4 A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, Pérez, PRX Quantum 4, 030314 (2023).
P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nature Commu- 34 H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall,

nications 5, 4213 (2014). Nature Communications 10, 3007 (2019).

5 M. Rossmannek, P. K. Barkoutsos, P. J. Ollitrault, and I. Tav- 35 H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J.

ernelli, The Journal of Chemical Physics 154, 114105 (2021). Mayhall, E. Barnes, and S. E. Economou, PRX Quantum 2,
6 B. Huang, M. Govoni, and G. Galli, PRX Quantum 3, 010339 020310 (2021).
(2022). 36 Z.-J. Zhang, T. H. Kyaw, J. S. Kottmann, M. Degroote, and
7 W. Li, Z. Huang, C. Cao, Y. Huang, Z. Shuai, X. Sun, J. Sun, A. Aspuru-Guzik, Quantum Science and Technology 6, 35001
X. Yuan, and D. Lv, Chemical Science 13, 8953 (2022). (2021).
8 M. Otten, M. R. Hermes, R. Pandharkar, Y. Alexeev, S. K. Gray, 37 P. G. Anastasiou, Y. Chen, N. J. Mayhall, E. Barnes, and S. E.

and L. Gagliardi, Journal of Chemical Theory and Computation Economou, Physical Review Research 6, 013254 (2024).
18, 7205 (2022). 38 J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and
9 C. Vorwerk, N. Sheng, M. Govoni, B. Huang, and G. Galli, Na- H. Neven, Nature Communications 9, 4812 (2018).
ture Computational Science 2, 424 (2022). 39 J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love,
10 B. Huang, N. Sheng, M. Govoni, and G. Galli, Journal of Chem- and A. Aspuru-Guzik, Quantum Science and Technology 4, 14008
ical Theory and Computation 19, 1487 (2023). (2018).
11 R. Izsák, C. Riplinger, N. S. Blunt, B. de Souza, N. Holzmann, 40 O. Higgott, D. Wang, and S. Brierley, Quantum 3, 156 (2019).

O. Crawford, J. Camps, F. Neese, and P. Schopf, Journal of Com- 41 M. Motta, W. Kirby, I. Liepuoniute, K. J. Sung, J. Cohn, A. Mez-

putational Chemistry 44, 406 (2023). zacapo, K. Klymko, N. Nguyen, N. Yoshioka, and J. E. Rice,
12 Y. Liu, O. R. Meitei, Z. E. Chin, A. Dutt, M. Tao, I. L. Chuang, arXiv:2312.00178v1 (2023).
and T. Van Voorhis, Journal of Chemical Theory and Computa- 42 Rossmannek, Max and Battaglia, Stefano, Qiskit Na-
tion 19, 2230–2247 (2023). ture and CP2K embedding, https://github.com/mrossinek/
13 M. Rossmannek, F. Pavošević, A. Rubio, and I. Tavernelli, The qiskit-nature-cp2k (2024), accessed on: 2024-4-26.
Journal of Physical Chemistry Letters 14, 3491 (2023). 43 At the time of writing.
14 T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. See- 44 V. Kapil, M. Rossi, O. Marsalek, R. Petraglia, Y. Litman,

wald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiff- T. Spura, B. Cheng, A. Cuzzocrea, R. H. Meißner, D. M. Wilkins,
mann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, B. A. Helfrecht, P. Juda, S. P. Bienvenue, W. Fang, J. Kessler,
V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Laz- I. Poltavsky, S. Vandenbrande, J. Wieme, C. Corminboeuf, T. D.
zaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Ander- Kühne, D. E. Manolopoulos, T. E. Markland, J. O. Richardson,
matt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belle- A. Tkatchenko, G. A. Tribello, V. Van Speybroeck, and M. Ce-
flamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C. J. Mundy, riotti, Computer Physics Communications 236, 214 (2019).
C. Plessl, M. Watkins, J. VandeVondele, M. Krack, and J. Hut- 45 Y. Chen, R. T. Williams, and W. A. Sibley, Physical Review 182,

ter, The Journal of Chemical Physics 152, 194103 (2020). 960 (1969).
15 Qiskit contributors, Qiskit: An open-source framework for quan- 46 L. A. Kappers, R. L. Kroes, and E. B. Hensley, Physical Review

tum computing, Zenodo, doi: 10.5281/zenodo.2573505 (2024). B 1, 4151 (1970).

16 The Qiskit Nature developers and contributors, Qiskit Nature, 47 G. P. Summers, T. M. Wilson, B. T. Jeffries, H. T. Tohver,

Zenodo, doi: 10.5281/zenodo.7828768 (2024). Y. Chen, and M. M. Abraham, Physical Review B 27, 1283
17 A. Savin, in Theoretical and Computational Chemistry, Vol. 4 (1983).
(Elsevier, 1996) pp. 327–357. 48 G. H. Rosenblatt, M. W. Rowe, G. P. Williams, R. T. Williams,
18 T. Leininger, H. Stoll, H.-J. Werner, and A. Savin, Chemical and Y. Chen, Physical Review B 39, 10309 (1989).
Physics Letters 275, 151 (1997). 49 C. Sousa and F. Illas, The Journal of Chemical Physics 115, 1435
19 R. Pollet, A. Savin, T. Leininger, and H. Stoll, The Journal of (2001).
Chemical Physics 116, 1250 (2002). 50 P. Rinke, A. Schleife, E. Kioupakis, A. Janotti, C. Rödl, F. Bech-
20 J. Toulouse, P. Gori-Giorgi, and A. Savin, Theoretical Chemistry stedt, M. Scheffler, and C. G. Van De Walle, Physical Review
Accounts 114, 305 (2005). Letters 108, 126404 (2012).
21 E. Fromager, J. Toulouse, and H. J. A. Jensen, The Journal of 51 E. Ertekin, L. K. Wagner, and J. C. Grossman, Physical Review

Chemical Physics 126, 074111 (2007). B 87, 155210 (2013).

22 E. D. Hedegård, J. Toulouse, and H. J. A. Jensen, The Journal 52 J. Strand, S. K. Chulkov, M. B. Watkins, and A. L. Shluger, The

of Chemical Physics 148, 214103 (2018). Journal of Chemical Physics 150, 044702 (2019).
23 K. Pernal and M. Hapka, WIREs Computational Molecular Sci- 53 A. Gallo, F. Hummel, A. Irmler, and A. Grüneis, The Journal of

ence 12, e1566 (2022). Chemical Physics 154, 064106 (2021).

24 G. Lippert, J. Hutter, and M. Parrinelloo, Molecular Physics 92, 54 B. T. G. Lau, B. Busemeyer, and T. C. Berkelbach, Optical prop-

477 (1997). erties of defects in solids via quantum embedding with good ac-
25 B. O. Roos, P. R. Taylor, and P. E. M. Siegbahn, Chemical tive space orbitals (2023), arXiv:2301.09668.
Physics 48, 157 (1980). 55 C. Vorwerk and G. Galli, Physical Review Materials 7, 033801
26 H.-J. Werner, Journal of Chemical Physics 74, 5794 (1981). (2023).
27 P. J. Ollitrault, A. Kandala, C.-F. Chen, P. K. Barkoutsos, 56 S. Verma, A. Mitra, Y. Jin, S. Haldar, C. Vorwerk, M. R. Hermes,

A. Mezzacapo, M. Pistoia, S. Sheldon, S. Woerner, J. M. Gam- G. Galli, and L. Gagliardi, The Journal of Physical Chemistry
betta, and I. Tavernelli, Physical Review Research 2, 43140 Letters 14, 7703 (2023).
(2020). 57 J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Let-
28 P. Jordan and E. Wigner, Zeitschrift für Physik 47, 631 (1928). ters 77, 3865 (1996).
29 S. B. Bravyi and A. Y. Kitaev, Annals of Physics 298, 210 (2002). 58 J. Hutter, BASIS, https://github.com/juerghutter/BASIS
30 S. Bravyi, J. M. Gambetta, A. Mezzacapo, and K. Temme, (2024), accessed on: 2024-04-17.
Tapering off qubits to simulate fermionic hamiltonians (2017), 59 S. Goedecker, M. Teter, and J. Hutter, Physical Review B 54,

arXiv:1701.08213. 1703 (1996).

 13

60 C. Hartwigsen, S. Goedecker, and J. Hutter, Physical Review B 68 S. Hirata and M. Head-Gordon, Chemical Physics Letters 314,
58, 3641 (1998). 291 (1999).
61 M. Krack, Theoretical Chemistry Accounts 114, 145 (2005). 69 A.-S. Hehn, B. Sertcan, F. Belleflamme, S. K. Chulkov, M. B.
62 J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, Watkins, and J. Hutter, Journal of Chemical Theory and Com-
T. Chassaing, and J. Hutter, Computer Physics Communications putation 18, 4186 (2022).
167, 103 (2005). 70 Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt,
63 J. P. Perdew and Y. Wang, Physical Review B 45, 13244 (1992). N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen,
64 P. M. W. Gill, R. D. Adamson, and J. A. Pople, Molecular Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Ko-
Physics 88, 1005 (1996). val, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain,
65 S. Paziani, S. Moroni, P. Gori-Giorgi, and G. B. Bachelet, Phys- M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto,
ical Review B 73, 155111 (2006). P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer,
66 M. Guidon, J. Hutter, and J. VandeVondele, Journal of Chemical H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upad-
Theory and Computation 5, 3010 (2009). hyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J.
67 D. M. Roessler and W. C. Walker, Physical Review 159, 733 Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkel-
(1967). bach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, The Journal
of Chemical Physics 153, 024109 (2020).