# Scalable Variational Quantum Circuits for Autoencoder-based Drug Discovery

Junde Li and Swaroop Ghosh  Department of Computer Science and Engineering
The Pennsylvania State University
{jul1512<EMAIL_ADDRESS>
###### Abstract

The de novo design of drug molecules is recognized as a time-consuming and
costly process, and computational approaches have been applied in each stage
of the drug discovery pipeline. Variational autoencoder is one of the
computer-aided design methods which explores the chemical space based on
existing molecular dataset. Quantum machine learning has emerged as an
atypical learning method that may speed up some classical learning tasks
because of its strong expressive power. However, near-term quantum computers
suffer from limited number of qubits which hinders the representation learning
in high dimensional spaces. We present a scalable quantum generative
autoencoder (SQ-VAE) for simultaneously reconstructing and sampling drug
molecules, and a corresponding vanilla variant (SQ-AE) for better
reconstruction. The architectural strategies in hybrid quantum classical
networks such as, adjustable quantum layer depth, heterogeneous learning
rates, and patched quantum circuits are proposed to learn high dimensional
dataset such as, ligand-targeted drugs. Extensive experimental results are
reported for different dimensions including 8x8 and 32x32 after choosing
suitable architectural strategies. The performance of quantum generative
autoencoder is compared with the corresponding classical counterpart
throughout all experiments. The results show that quantum computing advantages
can be achieved for normalized low-dimension molecules, and that high-
dimension molecules generated from quantum generative autoencoders have better
drug properties within the same learning period.

###### Index Terms:

Quantum Machine Learning, Variational Autoencoder, Drug Discovery

##  I Introduction

Traditional drug discovery and development pipeline, from concept and market,
takes 10-18 years and costs billions of dollars [1]. As the productivity for
small-molecule drugs declines, the drug discovery trend shifts toward biologic
drugs which predominantly target disease-selective cellular ligands (i.e.,
molecules that bind to certain atoms to form biological complexes) [2].
Ligand-targeted drugs exhibiting little affinity for health cells but high
affinity for pathologic cells can achieve required therapeutic potency with
minimal toxicity [3]. Ligands would interact with protein pockets (associated
with a disease) in the binding sites if they are structurally complementary.
One can find such ligands by searching through all viable chemical compounds
and molecules adhering to a given set of construction principles and boundary
conditions called the _chemical space_. Navigating this impractically large
chemical space falls within the field of de novo drug discovery [4].

Figure 1: (a) The drug discovery pipeline in 4 general phases (not including
development, clinical trials and launch phases); (b) the lock and key model
showing the binding process between protein pocket and ligand.

Computational approaches have been broadly applied to nearly every drug
discovery phases, especially molecule design with desirable properties [5]
which corresponds to the lead generation phase shown in Fig. 1(a). Generative
models such as, variational autoencoders (VAEs) [6] and generative adversarial
networks (GANs) [7] are adopted to generate a plethora of drug candidates (in
string-based and graph-based [8, 9] representations) for further screening by
adopting the lock and key model as shown in Fig. 1(b). The set-level
adversarial loss of GAN is calculated by comparing with the set of training
samples; while generative VAEs add extra latent variables from vanilla
autoencoders (AEs) and thus support instance-level reconstruction. Such
reconstruction is helpful for finding a ligand match given a receptor protein
[10, 11], by going through the molecular docking process as illustrated in
Fig. 1(b). AEs support more accurate reconstruction for the lack of latent
variables but do not support sampling new ligand molecules. Therefore, both
AEs and VAEs are adopted for reconstructing graph-based molecules but only
VAEs for sampling new ligand molecules through learning the latent space
distribution as shown in Fig. 2(a).

Figure 2: (a) The classical VAE architecture for drug molecule discovery
includes a decoder and an encoder which outputs Gaussian parameters for
approximating the latent space distribution (the only part that AE does not
involve) from which random noises are sampled for generating new molecules
through the decoder network (the generation or sampling process is denoted in
red box and arrows); (b) a patched variational quantum circuit (each patch has
three qubits) serves as the quantum encoder/decoder network with corresponding
quantum input, hidden and output layers. A hybrid VAE/AE is realized by
connecting the final fully connected (FC) layer with the preceding measurement
outcome vector.

Quantum computing has shown unique advantages such as, strong expressive power
across a variety of learning and optimization applications e.g., detection,
classification and discovery [12, 13, 14]. Quantum variational autoencoders
(QVAEs) had been primarily used for quantum state compression on adiabatic
quantum computers [15, 16] but rarely on gate-model quantum computers. When
being implemented, quantum VAE/AE is more challenging than the usual quantum
neural network (QNN) since it contains two major components, namely encoder
and decoder. The scaling to high-dimensional dataset becomes even more
challenging given limited number of qubits and the input-output mapping
constraint imposed by measurements interpreted as the output. The output
dimension is constrained to be less than the qubit count when measurement
expectation value is taken as output. We aim to address these scaling
difficulties of quantum autoencoders by exploring various architectural
techniques for drug discovery.

The actual QVAE architecture primarily depends on the dimension and scale of
training dataset. A fully quantum autoencoder can achieve quantum learning
advantages for normalized low dimension dataset, whereas a hybrid quantum
classical model is necessary to implement the high-dimensional and large-scale
dataset. Appropriate choices of quantum embedding methods also impact the
quality of uploading classical data into quantum states in a Hilbert space.
Besides, the novel QNN architecture with patched quantum circuits, as shown in
Fig. 2(b), is proposed to alleviate such dilemma. The strategy of adopting
heterogeneous learning rates for quantum and classical layers is explored for
counteracting the parameter range imbalance between the two layer types.
Appropriate quantum layer depth is explored to improve learning quality as
well. A suitable combination among these architectural factors is identified
for learning different dimensional datasets accordingly.

We propose a new qubit-efficient scalable quantum generative autoencoder (SQ-
VAE) to learn latent molecular distribution with reconstruction and sampling
fidelity on par with classical ones. The architecture of classical and quantum
VAEs are shown in Fig. 2. Extensive experiments on the architectural
techniques are performed using molecular datasets with different dimensions
including 8x8 QM9 [17] and 32x32 PDBbind ligand [18]. Additionally, two image
datasets, 8x8 Digits and 32x32 grayscale CIFAR-10 [19], are used to better
visualize the learning qualities from classical and quantum autoencoders.

Contributions: We, (i) develop a scalable gate-model quantum generative
autoencoder for reconstructing and sampling real-size drug molecules with
better drug properties; (ii) demonstrate quantum learning advantages in terms
of training speed and representation quality experimentally on the low-
dimensional digit and molecule datasets; (iii) propose several novel
architectural techniques including the patched quantum circuit and
heterogeneous learning rates to learn high-dimensional data distributions.

##  II Preliminaries

###  II-A Computational Drug Discovery

A variety of classical computational approaches using GANs and VAEs have been
explored for learning molecular distribution and generating new molecules from
the learned chemical space [9, 8]. In the present study, VAE is primarily
selected because it learns both the inference network and generator network,
the former of which helps with the instance-level matching between receptor
protein and ligand. Quantum generative models have also been recently studied
for leveraging the strong expressive power and intrinsic probabilistic nature
of quantum circuits [11, 13]. However, only small molecule discovery was
targeted in these studies. One small molecular graph from QM9 and its matrix
representation are shown in Fig. 3. Each diagonal element in the molecule
matrix denotes the encoded atom type (1-C, 2-N, 3-O), and each off-diagonal
element represents the encoded bond type (0-NONE, 1-SINGLE, 2-DOUBLE,
4-AROMATIC) in the molecular graph. Only heavy atoms excluding Hydrogen are
displayed in the matrix.

###  II-B Variational Autoencoder

Autoencoder (AE) consists of an encoder network, which converts an input
vector to a code vector, and a decoder network which reconstructs the original
input with high fidelity [20]. Its generative version of VAE additionally
learns a variational model for the latent variables to capture the underlying
sample distribution at the expense of less accurate reconstruction.
Mathematically, VAE is represented by an inference network (a.k.a. encoder)
$q_{\boldsymbol{\phi}}(\boldsymbol{z}|\boldsymbol{x})$ and a generator network
(a.k.a. decoder) $p_{\boldsymbol{\theta}}(\boldsymbol{x}|\boldsymbol{z})$
where $\boldsymbol{z}$ denotes the latent variables [6]. The inference network
represented by a neural network outputs parameters for learning the Gaussian
distribution
$q_{\boldsymbol{\phi}}(\boldsymbol{z}|\boldsymbol{x})=\mathcal{N}(\mu_{\phi}(\boldsymbol{x}),\boldsymbol{\Sigma}_{\boldsymbol{\phi}}(\boldsymbol{x}))$.
The intractable true posterior
$p_{\boldsymbol{\theta}}(\boldsymbol{z}|\boldsymbol{x})$ is approximated by
the encoder network. The generative model depends on a continuous variable
$\boldsymbol{z}$ and the posterior. Assuming a training set
$S=\\{\boldsymbol{x}_{i}\\}_{i=1}^{N}$, the parameters $\boldsymbol{\phi}$ and
$\boldsymbol{\theta}$ are jointly learned by minimizing the negative of
evidence lower bound (ELBO) as follows:

$\displaystyle\min_{\boldsymbol{\theta},\boldsymbol{\phi}}\mathcal{L}_{\text{ELBO}}(\boldsymbol{\theta},\boldsymbol{\phi};\boldsymbol{x})=$
$\displaystyle-\mathbb{E}_{\boldsymbol{z}\sim
q_{\phi}(\boldsymbol{z}|\boldsymbol{x})}[\log{(p_{\boldsymbol{\theta}}(\boldsymbol{x}|\boldsymbol{z}))}]+D_{KL}[q_{\boldsymbol{\phi}}(\boldsymbol{z}|\boldsymbol{x})||p(\boldsymbol{z})]$

where $p(z)$ is the prior assumed to be $\mathcal{N}(0,I)$. The first term
$-\mathbb{E}_{\boldsymbol{z}\sim
q_{\phi}(\boldsymbol{z}|\boldsymbol{x})}[\log{(p_{\boldsymbol{\theta}}(\boldsymbol{x}|\boldsymbol{z}))}]$
in above equation is known as reconstruction loss, and the second term
$D_{KL}(q||p)$ is the Kullback-Leibler (KL) divergence between the
approximation density $q_{\boldsymbol{\phi}}(\boldsymbol{z}|\boldsymbol{x})$
and the Gaussian prior. Continuous random noise is drawn from
$p(\boldsymbol{z})$ for further approximating sampling new data points from
$p_{\text{data}}(\boldsymbol{x})$.

Figure 3: (a) A sample molecular graph where each circled number denotes the
atom sequence; (b) the corresponding molecule matrix with diagonal elements
representing atom types and off-diagonal ones representing bond types.

###  II-C Quantum Neural Networks

Quantum gates are quantum logic operations on a small number of qubits. A
single-qubit operation behaves like a matrix-vector multiplication, and two-
qubit operation are employed for creating quantum entanglement. Rotation gates
are adopted for variational learning, and each rotation gate contains three
trainable parameters. Moreover, non-parametric CNOT gate (see Fig. 2(b)) is
adopted for creating entanglement, and Pauli-Z gate for measuring the
expectation value for each qubit.

Quantum embedding uploads a given classical data sample $x$ into its
corresponding quantum state $\ket{\psi_{x}}$ in a Hilbert space. Embedding is
a crucial part of quantum circuit which affects the computational power. For
the vector $\boldsymbol{x}\in\mathcal{R}^{d}$, the amplitude encoded state is
defined as
$\ket{x}=\frac{1}{\|\boldsymbol{x}\|_{2}}\sum_{j=1}^{d}{x_{j}\ket{j}}$. Angle
embedding maps the feature vector $\boldsymbol{x}$ to the angles
$\phi(\boldsymbol{x})=(\phi_{1}(x_{1}),\phi_{2}(x_{2}),\ldots,\phi_{d}(x_{d}))$
with which $d$ qubits have to be rotated to reflect each feature value
$x_{i}$. Amplitude embedding has the stringent input-output constraint,
whereas angle embedding is not qubit efficient since it requires one qubit per
feature.

Measurement provides classical information in terms of expectation, variance,
or probabilities from quantum system. For expectations measured on the state
$\ket{\psi}=a\ket{0}+b\ket{1}$, the average of the observable $\textit{Z}=$
$\big{(}\begin{smallmatrix}1&0\\\ 0&-1\end{smallmatrix}\big{)}$ under the
state mathematically refers to the following:

$\displaystyle\braket{\textit{Z}}_{\ket{\psi}}\equiv\bra{\psi}\textit{Z}\ket{\psi}\equiv\text{Tr}[\ket{\psi}\bra{\psi}\textit{Z}]=|a|^{2}-|b|^{2}\in[-1,1]$

where the expectation bound $[-1,1]$ is obtained since the squared norm of
amplitude represents a probability. For probabilities measured on the state,
an array containing the probability $|\braket{i}{\psi}|^{2}$ for each
computational basis state $i$ is returned, satisfying
$\sum_{i}|\braket{i}{\psi}|^{2}=1$.

Figure 4: (a) Comparison of train MSE losses between VAE and BQ-VAE on
original Digits QM9 molecules; (b) comparison of train MSE losses on
normalized Digits and QM9 molecules by dividing the respective L1 norm for all
data points; (c) three sample digits, their respective BQ-VAE reconstructions,
and three random samples from the learned generator network; (d) one sample
QM9 molecule and its BQ-VAE reconstructions from input (original) and L1
normalized molecule, respectively.

##  III Variational Quantum Autoencoder

###  III-A Architectures

The general QNN architecture consists of quantum embedding, repeatable
variational entanglement layers, and measurement for constructing quantum
generative and vanilla autoencoders. However, the specific choices of quantum
embedding and measurement largely depend on the dimensions and scales of given
dataset. Fully quantum autoencoders only apply to normalized-scale
reconstruction since both expectation and probabilities have the upper bound
of 1. High-dimensional data reconstruction would require the patched
architecture to avoid the tradeoff between qubit-efficiency and input-to-
output mapping constraint. For simplicity, the repeatable hidden layer is
fixed with rotation gates $R(\psi,\theta,\omega)$ acting on each qubit,
followed by CNOT gates with a periodic layout, as shown in Fig. 2(b). We
demonstrate two quantum autoencoder architectures, i.e., baseline quantum
VAE/AE (BQ-VAE/AE) and scalable quantum VAE/AE (SQ-VAE/AE) to learn the low-
dimensional QM9 and high-dimensional PDBbind molecular datasets, respectively.

###  III-B Baseline Quantum Autoencoders

BQ-VAE/AE adopts amplitude embedding and expectation output for encoder, and
angle embedding and probability output for decoder. Low resolution 8x8 Digits
set is used to compare BQ-VAE/AE performance relative to the classical
counterpart. To demonstrate fully baseline quantum VAEs (F-BQ-VAEs), QM9
molecules and image digits are first normalized by directly dividing each non-
negative feature value by their sum. Latent space dimension (LSD) is
$log_{2}(64)=6$ by measuring expectation of each qubit. The LSD vector from
the latent space is reconstructed into $2^{6}=64$ dimension by computing the
probability of basis states. For original-scale (non-normalized) data points,
a final fully connected (FC) classical layer is used to reconstruct feature
values, resulting in the hybrid quantum variant (H-BQ-VAE/AE). Quantum hidden
layers are repeated for $L=3$ times in BQ-VAEs/AEs. Similarly, classical
encoder also takes 3 hidden linear layers followed by ReLU activation for
reducing the dimensions to 32, 16, and 6, respectively, while classical
decoder converts the dimensions in a reversed order. Table I compares of
number of trainable parameters between classical and baseline quantum
autoencoders.

TABLE I: Comparison of number of trainable parameters. Parameter Type | VAE(AE) | F-BQ-VAE(AE) | H-BQ-VAE(AE)
---|---|---|---
Quantum | 0 (0) | 108 (108) | 108 (108)
Classical | 5694 (5610) | 84 (0) | 4286 (4202)
Total | 5694 (5610) | 192 (108) | 4394 (4310)

Large number of classical parameters in H-BQ-VAE is due to the final linear
layer with a [64, 64] transformation. Apart from using fewer parameters, BQ-
VAE/AE even learns faster for normalized QM9 molecules in terms of the number
of training epochs as shown in Fig. 4(b). However, Fig. 4(a) does not show
quantum computational advantage for original-scale dataset due to the
probability output constraint. Three normalized input samples from Digits
dataset, their respective BQ-VAE reconstructions, and three randomly sampled
digits from the learned generator network are shown in Fig. 4(c). Note that,
normalized-scale digit images have the same quality as original-scale ones
when being visualized. However, molecule reconstructions from original and
normalized inputs differ so much that the latter hardly shares characteristics
with the input molecule, as shown in Fig. 4(d). This indicates that
normalization is not the feasible quantum autoencoder approach for learning
QM9 molecular representation and high-dimensional PDBbind molecules.

###  III-C Scalable Quantum Autoencoders

Scalable quantum autoencoders including SQ-VAEs and SQ-AEs are developed for
near-term quantum computers to learn high-dimensional dataset such as, PDBbind
ligands filtered with up to 32 heavy molecules. As illustrated in Fig. 5,
direct application of baseline quantum autoencoders to the 1024-dimensional
(32x32) data points do not work well due to the stringent input-to-output
mapping, i.e., latent space being too compact ($10=\log_{2}1024$ dimensional)
to faithfully represent PDBbind ligands. The fully quantum variant, F-BQ-AE,
hardly learns because probability values cannot be matched to the original-
scale ligand matrices. The extra classical layer in the hybrid quantum model
(H-BQ-AE) functions as the transformation to original scales. Fig. 5 (b) shows
lower test losses for classical AEs with larger latent space dimensions
whereas VAEs almost remain unchanged. Scalable quantum autoencoder relies on
the patched quantum circuits for scaling to high-dimensional dataset given
small size quantum computers. Angle embedding is used to encode the latent
variable $\boldsymbol{z}$ into the corresponding quantum state and measurement
expectation value is returned as the output since the probabilities from 1024
basis states are too miniscule to be reconstructed.

Patched quantum circuit (see Fig. 2(b) is derived from the patch quantum GAN
[21] which takes all features into each of the multiple circuits to obtain
richer representation of fake digits. However, we partition the entire feature
vector into multiple equal-sized sub-vectors, and each sub-vector is fed into
a quantum sub-circuit. Therefore, our revised method occupies fewer qubits in
each sub-circuit by embedding only a portion of the feature vector. Moreover,
patched quantum circuit increases the output dimension by taking more
expectation values as output. This alleviates the issue of having too compact
latent space from directly using baseline quantum autoencoder. It is worth
mentioning that the architecture advantage in terms of model complexity of SQ-
VAE/AE over its classical counterpart becomes more pronounced, relative to the
comparison aforementioned in Table I, because of the higher dimensional input.
The advantage here is not re-tabulated for brevity.

Figure 5: (a) Two BQ-AEs and the classical AE with three hidden layers applied
to high-dimensional PDBbind ligands; (b) AEs with increasingly large latent
space show lower test losses while VAEs remain almost unchanged.

Aside from the patched circuit, we also present three following architecture-
related factors that could affect SQ-VAE/AE performance. Firstly, the number
of quantum hidden layers $L$ determines the number of trainable weights which
relate to the trainability and expressiveness of variational quantum circuits.
A suitable circuit depth is selected after a series of sensitivity studies.
Secondly, similar to Fig. 5(b), the learning quality of SQ-VAE/AE also depends
on latent space dimension which in turn affects the number of patches in the
circuit. Thirdly, quantum parameters fall in the range $[-\pi,\pi]$, whereas
classical parameter space is much more vast given equal number of parameters.
The small-scale of quantum parameters may require a different quantum learning
rate schedule from classical one. This assumption is validated by adopting
heterogeneous learning rates for quantum and classical layers.

##  IV Experiments

This section discusses the experiment results of SQ-VAE/AE for learning high-
dimensional molecules since the learning result of baseline quantum
autoencoders has already been demonstrated on QM9 molecules in Section III-B.

###  IV-A Dataset and Metrics

The refined PDBbind 2019 dataset [18] contains 4852 protein-ligand complexes,
whereas only ligand molecules are used to demonstrate the SQ-VAE/AE. We filter
out molecules with more than 32 heavy atoms since 32x32 (a power of 2)
dimension can directly leverage quantum information processing. Similar to QM9
dataset (see the 9x9 example in Fig. 3), PDBbind ligands are also represented
by molecule matrices where diagonal elements indicate the encoded atom types
(1-C, 2-N, 3-O, 4-F, 5-S) and off-diagonal elements indicate the encoded bond
types. Each molecule matrix is symmetric, therefore, half of the off-diagonal
elements could be removed. Ligands consisting of other atom types are removed.
Therefore, only 2492 ligands are finally taken for learning SQ-VAE/AE. The
dataset is split into 85% and 15% for train and test sets, respectively. SQ-
VAE/AE performance on PDBbind ligands is evaluated using the metric of MSE
loss together with three molecule property metrics of quantitative estimate of
druglikeness (QED), log octanol-water partition coefficient (logP) and
synthetic accessibility (SA) retrieved using RDKit [22]. Additionally, 32x32
resolution gray-scaled CIFAR-10 [19] images are adopted for visualizing SQ-AE
reconstruction.

###  IV-B Implementation Details

Patched quantum circuit is included to scale the SQ-VAE/AE network to high-
dimensional ligands. Other architecture-related factors i.e., adjustable
quantum layer depth and heterogeneous learning rates are finalized through
sensitivity studies. Both quantum encoder and decoder are connected to a
classical layer to map the measurement to original ligand features.

All scalable quantum autoencoder variants are trained with a mini-batch of 32
molecules using the Adam optimizer with default $\beta_{1}=0.9$ and
$\beta_{2}=0.999$ on a single RTX 2080 Ti GPU and the PennyLane simulation
platform for the quantum stage. The learning rate is set to 0.001 throughout
20 training epochs for tuning quantum layer depth. A set of data points is
randomly selected from the test set to evaluate SQ-VAE/AE reconstructions.
Sampling results from SQ-VAEs are examined by sampling Gaussian noise from
learned latent spaces.

Figure 6: Training and test reconstruction MSE loss comparison among SQ-VAEs
with 1 to 9 strongly entangling layer(s).

###  IV-C Ablation Study

The number of spurious local minima depends exponentially on the number of
quantum parameters [23]. Therefore, the suitable quantum layer depth needs to
be identified for better performance. Adopting heterogeneous learning rates
for quantum and classical layers can speed up the overall learning.

The suitable quantum layer depth is identified by sweeping only the SQ-AE
layers from 1 to 9. Note that SQ-VAE has the same architecture as SQ-AE except
the extra latent variables. The losses after 5 and 10 epochs are shown in Fig.
6 where the configuration with 5 entangling layers produce the lowest test
losses. The result indicates that SQ-AE with too few quantum layers hurts its
expressive power, whereas too many layers possibly create unwanted number of
spurious local minima. Thus the following SQ-VAE/AE experiments are all
configured with 5 layers.

A single loss term is used to iteratively update the parameters in both
quantum and classical layers. However, the quantum rotation angle parameter
space is quite different from the classical parameter space. Therefore,
studies of heterogenous learning rates are conducted only for SQ-AEs, due to
architecture similarity with SQ-VAEs, to potentially balance the learning
qualities in the hybrid algorithms. Five learning rates i.e., [0.001, 0.003,
0.01, 0.03, 0.1] are examined for both quantum and classical layers. Fig. 7
shows the sensitivity results of total 25 combinations. Classical learning
rate of 0.01 and quantum learning rate of 0.03 produce the most favorable
training loss. All the following experiments are configured with this
combination.

Figure 7: Training loss comparison among SQ-VAEs configured with different
combinations of quantum and classical learning rates.

###  IV-D Evaluation Results

Figure 8: (a) Comparison of train losses between VAE and SQ-VAE/AE with
different LSDs on PDBbind ligands; (b) comparison of train losses between
VAE/AE and SQ-VAE/AE on grayscale CIFAR-10 images; (c) three sample CIFAR-10
images and their respective reconstructions from AE and SQ-AE.

As shown in Fig. 5(a), baseline quantum autoencoders could not learn the high-
dimensional PDBbind ligands well due to the input-output mapping constraint of
the 10 LSD. Larger intermediate dimensions up to 64 improves the learning
result for classical AEs (Fig. 5(b)). However, input dimension and amplitude
embedding determine the latent space dimension for quantum circuits. To
increase LSD, patched quantum circuit in the scalable quantum autoencoder
splits the holistic circuit into multiple patches. The performances of SQ-AEs
with 4 different patches are evaluated. The variant with $p$ patches has
$p\log_{2}(1024/p)$ dimensions for latent space. Therefore, the LSD increases
from 10 in baseline quantum autoencoders to 18, 32, 56 and 96 in scalable
quantum autoencoders with 2, 4, 8 and 16 patches, respectively. Classical AEs
with corresponding LSDs are created as well for comparison. To evaluate the
sampling results using drug property metrics, SQ-VAE variants with
corresponding patches are also conducted even though the negative MSE results
are observed for SQ-VAEs with larger LSDs as shown in Fig. 5(b).

Table II compares the drug properties of 1000 new molecules generated from SQ-
VAEs and VAEs with different LSDs after 20-epoch learning on PDBbind ligands.
SQ-VAE with 18 LSD generates molecules with better logP and SA scores. The
best QED score is achieved at SQ-VAE with 56 LSD using 8 circuit patches in
the hybrid algorithm. VAEs with larger LSDs show slightly worse MSE losses as
shown in Fig. 8(a), however, logP and SA scores increase considerably with
increasing LSDs. SQ-VAE drug properties do not vary much with LSD.

TABLE II: Drug properties of sampled ligands from SQ-VAEs and VAEs with different LSDs.  Metrics | LSD-18 | LSD-32 | LSD-56 | LSD-96
---|---|---|---|---
VAE-QED | 0.138 | 0.179 | 0.139 | 0.142
SQ-VAE-QED | 0.153 | 0.177 | 0.204 | 0.167
VAE-logP | 0.357 | 0.472 | 0.496 | 0.761
SQ-VAE-logP | 0.780 | 0.616 | 0.709 | 0.740
VAE-SA | 0.192 | 0.292 | 0.307 | 0.599
SQ-VAE-SA | 0.626 | 0.479 | 0.534 | 0.547

Sampling results can only be evaluated using generative autoencoders, while
reconstruction can be evaluated using both generative and vanilla
autoencoders. Reconstruction results from SQ-VAE/AE are approximately
reflected by loss bars displayed in Fig. 8(a). Besides, reconstructions from
autoencoders with 18 LSD are also evaluated using CIFAR images as shown in
panels (b-c). SQ-VAE/AE shows reconstruction results on par with classical
counterparts. Reconstructed images in panel (c) only show the sketches of
original inputs for both classical and quantum AEs after learning in 20
epochs. The simulation inefficiency of quantum circuits precludes a longer
learning period which leads to more accurate reconstruction.

##  V Conclusion

This paper introduces a baseline quantum autoencoder that shows computational
advantages for only low-dimensional molecules. The scalable quantum variants
leverage the patched circuits and heterogeneous learning rates for scaling to
high-dimensional ligands and show better sampling results in terms of drug
properties and reconstruction results comparable to classical counterparts.
The proposed scalable quantum autoencoder also applies to other tasks such as
image generation.

## References

  * [1] S. Myers and A. Baker, “Drug discovery—an operating model for a new era,” _Nature biotechnology_ , vol. 19, no. 8, pp. 727–730, 2001.
  * [2] W. Li, I.-H. Pang, M. T. F. Pacheco, and H. Tian, “Ligandomics: a paradigm shift in biological drug discovery,” _Drug discovery today_ , vol. 23, no. 3, pp. 636–643, 2018.
  * [3] M. Srinivasarao and P. S. Low, “Ligand-targeted drug delivery,” _Chemical reviews_ , vol. 117, no. 19, pp. 12 133–12 164, 2017.
  * [4] J.-L. Reymond and M. Awale, “Exploring chemical space for drug discovery using the chemical universe database,” _ACS chemical neuroscience_ , vol. 3, no. 9, pp. 649–657, 2012.
  * [5] S. Ekins, A. C. Puhl, K. M. Zorn, T. R. Lane, D. P. Russo, J. J. Klein, A. J. Hickey, and A. M. Clark, “Exploiting machine learning for end-to-end drug discovery and development,” _Nature materials_ , vol. 18, no. 5, p. 435, 2019\.
  * [6] D. P. Kingma and M. Welling, “Auto-encoding variational bayes.” in _Proceedings of the International Conference on Learning Representations (ICLR)_ , 2014.
  * [7] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” in _Advances in neural information processing systems_ , 2014.
  * [8] M. J. Kusner, B. Paige, and J. M. Hernández-Lobato, “Grammar variational autoencoder,” in _International Conference on Machine Learning_. PMLR, 2017, pp. 1945–1954.
  * [9] N. De Cao and T. Kipf, “Molgan: An implicit generative model for small molecular graphs,” _arXiv preprint arXiv:1805.11973_ , 2018.
  * [10] G. Parmar, D. Li, K. Lee, and Z. Tu, “Dual contradistinctive generative autoencoder,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 823–832.
  * [11] J. Li, M. Alam, C. M. Sha, J. Wang, N. V. Dokholyan, and S. Ghosh, “Drug discovery approaches using quantum machine learning,” _arXiv preprint arXiv:2104.00746_ , 2021.
  * [12] Y. Du, M.-H. Hsieh, T. Liu, D. Tao, and N. Liu, “Quantum noise protects quantum classifiers against adversaries,” _Physical Review Research_ , vol. 3, no. 2, p. 023153, 2021.
  * [13] J. Li, R. Topaloglu, and S. Ghosh, “Quantum generative models for small molecule drug discovery,” _arXiv preprint arXiv:2101.03438_ , 2021.
  * [14] A. Abbas, D. Sutter, C. Zoufal, A. Lucchi, A. Figalli, and S. Woerner, “The power of quantum neural networks,” _Nature Computational Science_ , vol. 1, no. 6, pp. 403–409, 2021.
  * [15] J. Romero, J. P. Olson, and A. Aspuru-Guzik, “Quantum autoencoders for efficient compression of quantum data,” _Quantum Science and Technology_ , vol. 2, no. 4, p. 045001, 2017.
  * [16] A. Khoshaman, W. Vinci, B. Denis, E. Andriyash, H. Sadeghi, and M. H. Amin, “Quantum variational autoencoder,” _Quantum Science and Technology_ , vol. 4, no. 1, p. 014001, 2018.
  * [17] R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A. von Lilienfeld, “Quantum chemistry structures and properties of 134 kilo molecules,” _Scientific Data_ , vol. 1, 2014.
  * [18] Z. Liu, Y. Li, L. Han, J. Li, J. Liu, Z. Zhao, W. Nie, Y. Liu, and R. Wang, “Pdb-wide collection of binding data: current status of the pdbbind database,” _Bioinformatics_ , vol. 31, no. 3, pp. 405–412, 2015.
  * [19] A. Krizhevsky, G. Hinton _et al._ , “Learning multiple layers of features from tiny images,” 2009.
  * [20] G. E. Hinton and R. S. Zemel, “Autoencoders, minimum description length, and helmholtz free energy,” _Advances in neural information processing systems_ , vol. 6, pp. 3–10, 1994.
  * [21] H.-L. Huang, Y. Du, M. Gong, Y. Zhao, Y. Wu, C. Wang, S. Li, F. Liang, J. Lin, Y. Xu _et al._ , “Experimental quantum generative adversarial networks for image generation,” _arXiv preprint arXiv:2010.06201_ , 2020.
  * [22] “Rdkit: Open-source cheminformatics software,” _https://www.rdkit.org/_.
  * [23] X. You and X. Wu, “Exponentially many local minima in quantum neural networks,” in _Proceedings of the 38th International Conference on Machine Learning_ , 2021.


$\displaystyle\frac{1}{\sqrt{2}}\sum_{\tilde{i}}\tilde{c}_{\tilde{i}}\bigl{[}{\bm{\vartheta}}_{\tilde{i}}({\bf
r}_{1})\otimes{\bm{\vartheta}}_{\tilde{\tilde{i}}}({\bf r}_{2})$
$\displaystyle{}-{\bm{\vartheta}}_{\tilde{i}}({\bf
r}_{2})\otimes{\bm{\vartheta}}_{\tilde{\tilde{i}}}({\bf r}_{1})\bigr{]},$
(123)

where the index ${\tilde{i}}$ runs only over half of the PS, and we have
introduced the ${\bm{\vartheta}}$ spinors that form Schmidt pairs
($\tilde{i},\tilde{\tilde{i}}$). Note that these Schmidt orbitals are only
defined up to a unitary transformation within the degenerate subspace.

The 1-RDM of the wavefunction given by 123 reads as

${\bf n}_{1}^{+,2e-}({\bf r},{\bf
r}^{\prime})=\sum_{\tilde{i}}\lvert\tilde{c}_{i}\rvert^{2}_{\tilde{i}}\left[{\bm{\vartheta}}_{\tilde{i}}({\bf
r}){\bm{\vartheta}}^{\dagger}_{\tilde{i}}({\bf
r}^{\prime})+{\bm{\vartheta}}_{\tilde{\tilde{i}}}({\bf
r}){\bm{\vartheta}}^{\dagger}_{\tilde{\tilde{i}}}({\bf r}^{\prime})\right],$
(124)

which is already given in its diagonal representation and from which we find
$\lvert\tilde{c}_{i}\rvert^{2}=n_{\tilde{i}}$

When the Hamiltonian operator preserves time-reversal symmetry (i.e. it is
time-independent and does not contain external magnetic fields) it commutes
with the Kramers’ operator 181818${\widehat{\mathcal{K}}}$ in the many
electron case leads to the complex conjugate of the CI coefficients and
transforms all spinors in the tensor product [141, 55]. This implies that we
can choose all our eigenstates to satisfy [141]

${\Psi}^{2e-}_{+}({\bf r}_{1},{\bf
r}_{2})=\widehat{\mathcal{K}}\Psi^{2e-}_{+}({\bf r}_{1},{\bf r}_{2}).$ (125)

which has the right symmetry for real CI coefficients. This means that the
degenerate configurations in (123) transform into each other. Since the
Schmidt orbitals are only defined up to a unitary transformation with the
degenerate subspace, we can use this freedom to transform the Schmidt pairs to
Kramers pairs. So the Schmidt pairs can be replaced with Kramers pairs

$\displaystyle\Psi_{+}^{2e-}({\bf r}_{1},{\bf r}_{2})={}$
$\displaystyle\frac{1}{\sqrt{2}}\sum_{i}c_{i}\bigl{[}\widetilde{\bm{\chi}}_{i}({\bf
r}_{1})\otimes\widetilde{\bm{\chi}}_{\bar{i}}({\bf r}_{2})$
$\displaystyle{}-\widetilde{\bm{\chi}}_{i}({\bf
r}_{2})\otimes\widetilde{\bm{\chi}}_{\bar{i}}({\bf r}_{1})\bigr{]},$ (126)

where

$c_{i}=\sqrt{n_{i}}\textrm{e}^{-\textrm{i}\zeta_{i}}.$ (127)

The phases $\zeta_{i}=k\pi$ with $k\in\mathbb{Z}$ to preserve the Kramers’
symmetry.

The 1-RDM is given by

${\bf n}_{1}^{+,2e-}({\bf r},{\bf
r}^{\prime})=\sum_{i}n_{i}\left[\widetilde{\bm{\chi}}_{i}({\bf
r})\widetilde{\bm{\chi}}^{\dagger}_{i}({\bf
r}^{\prime})+\widetilde{\bm{\chi}}_{\bar{i}}({\bf
r})\widetilde{\bm{\chi}}^{\dagger}_{\bar{i}}({\bf r}^{\prime})\right].$ (128)

Evaluating the electron-electron interaction using Eq. (126) we arrive to

$\displaystyle\langle\Psi_{+}^{2e-}|\widehat{\widetilde{W}}|\Psi_{+}^{2e-}\rangle$
$\displaystyle={}\sum_{i,j}\frac{\sqrt{n_{i}n_{j}}}{2}\textrm{e}^{\textrm{i}(\zeta_{j}-\zeta_{i})}\int
d{\bf r}_{1}d{\bf
r}_{2}\textrm{Tr}\left[\frac{1}{r_{12}}\left(\mathbb{I}_{16\times
16}-{\bm{\alpha}}\odot{\bm{\alpha}}\right)\right.$
$\displaystyle{}\times\left((\widetilde{{\bm{\chi}}}_{i}({\bf
r}_{1})\otimes\widetilde{{\bm{\chi}}}_{\bar{i}}({\bf
r}_{2})(\widetilde{{\bm{\chi}}}^{\dagger}_{j}({\bf
r}_{1})\otimes\widetilde{{\bm{\chi}}}^{\dagger}_{\bar{j}}({\bf
r}_{2}))-(\widetilde{{\bm{\chi}}}_{i}({\bf
r}_{1})\otimes\widetilde{{\bm{\chi}}}_{\bar{i}}({\bf
r}_{2}))(\widetilde{{\bm{\chi}}}^{\dagger}_{j}({\bf
r}_{2})\otimes\widetilde{{\bm{\chi}}}^{\dagger}_{\bar{j}}({\bf
r}_{1}))\right.$
$\displaystyle{}-\left.\left.(\widetilde{{\bm{\chi}}}_{i}({\bf
r}_{2})\otimes\widetilde{{\bm{\chi}}}_{\bar{i}}({\bf
r}_{1}))(\widetilde{{\bm{\chi}}}^{\dagger}_{j}({\bf
r}_{1})\otimes\widetilde{{\bm{\chi}}}^{\dagger}_{\bar{j}}({\bf
r}_{2}))+(\widetilde{{\bm{\chi}}}_{i}({\bf
r}_{2})\otimes\widetilde{{\bm{\chi}}}_{\bar{i}}({\bf
r}_{1}))(\widetilde{{\bm{\chi}}}^{\dagger}_{j}({\bf
r}_{2})\otimes\widetilde{{\bm{\chi}}}^{\dagger}_{\bar{j}}({\bf
r}_{1}))\right)\right],$ (129)

where $\mathbb{I}_{16\times 16}$ is the 16$\times$16 identity matrix. In this
form, we notice that the energy minimization procedure also implies the
optimization w.r.t. the phases ($\\{\zeta_{i}\\}$). Consequently, fixing the
phases as in the nonrelativistic case and exploding Kramers’ symmetry on the
integrals we arrive to the energy expression given by Eq. (III.6). Therefore,
the rel-GNOF/PNOF$x$ functionals are equivalent to the relativistic version of
the Fixed-Phases functional for ‘spin-compensated’ two-electron systems.

Finally, let us show that the relativistic Fixed-Phases functional contains
its nonrelativistic counterpart. To that end, let us first focus on the NOs
forming the Kramers’ pairs, which are built as four component spinors, i.e.

$\widetilde{\bm{\chi}}_{i}({\bf r})=\begin{pmatrix}\widetilde{\phi}_{i,1}({\bf
r})\\\ \widetilde{\phi}_{i,2}({\bf r})\\\ \widetilde{\phi}_{i,3}({\bf r})\\\
\widetilde{\phi}_{i,4}({\bf r})\end{pmatrix}$ (130)

and

$\widehat{\mathcal{K}}\widetilde{\bm{\chi}}_{i}({\bf
r})=\begin{pmatrix}[r]-\widetilde{\phi}^{*}_{i,2}({\bf r})\\\
\widetilde{\phi}^{*}_{i,1}({\bf r})\\\ -\widetilde{\phi}^{*}_{i,4}({\bf r})\\\
\widetilde{\phi}^{*}_{i,3}({\bf
r})\end{pmatrix}=\widetilde{\bm{\chi}}_{\bar{i}}({\bf r}).$ (131)

It is easy to prove that the tensor product of 4-component NOs contains
singlet and triplet contributions, but we can introduce some constraints to
remove the triplet contribution. To that end, let us assume that the so-called
small component terms of the NOs are negligible (i.e.
$\widetilde{\phi}_{i,3}({\bf r})$ and $\widetilde{\phi}_{i,4}({\bf
r})\rightarrow 0$) 191919This is normally the case in the nonrelativistic
limit [55]. Then, adding the following constraint $\widetilde{\phi}_{i,2}({\bf
r})=a\widetilde{\phi}_{i,1}({\bf r})$ with $a\in\mathbb{C}$ we may write the
components as $\widetilde{\phi}_{i,1}({\bf r})=\Re\widetilde{\phi}_{i,1}({\bf
r})+\textrm{i}\Im\widetilde{\phi}_{i,1}({\bf r})$ and
$\widetilde{\phi}_{i,2}({\bf r})=a{}\Re\widetilde{\phi}_{i,1}({\bf
r})+a\textrm{i}\Im\widetilde{\phi}_{i,1}({\bf r})$, where
$\Re\widetilde{\phi}_{i,1}({\bf
r})=\textrm{Real}\left[\widetilde{\phi}_{i,1}({\bf r})\right]$ and
$\Im\widetilde{\phi}_{i,1}({\bf
r})=\textrm{Imaginary}\left[\widetilde{\phi}_{i,1}({\bf r})\right]$. Thence,
we arrive to the following expression for the NOs

$\widetilde{\bm{\chi}}_{i}({\bf
r})=\begin{pmatrix}\Re\widetilde{\phi}_{i,1}({\bf
r})+\textrm{i}\Im\widetilde{\phi}_{i,1}({\bf r})\\\
a\Re\widetilde{\phi}_{i,1}({\bf r})+a\textrm{i}\Im\widetilde{\phi}_{i,1}({\bf
r})\\\ \end{pmatrix}$ (132)

and

$\widetilde{\bm{\chi}}_{\bar{i}}({\bf
r})=\begin{pmatrix}-a^{*}\Re\widetilde{\phi}_{i,1}({\bf
r})+a^{*}\textrm{i}\Im\widetilde{\phi}_{i,1}({\bf r})\\\
\Re\widetilde{\phi}_{i,1}({\bf r})-\textrm{i}\Im\widetilde{\phi}_{i,1}({\bf
r})\\\ \end{pmatrix},$ (133)

where $a^{*}$ is the complex-conjugate of $a$ and we have already omitted the
small component terms.

Assuming that the scalar spinors are real (i.e.
$\Im\widetilde{\phi}_{i,1}({\bf r})=0$) and inserting the NOs in Eq. (126) we
arrive to 202020Retaining only the imaginary part of the scalar spinors
(setting ${\Re\widetilde{\phi}_{i,1}({\bf r})=0}$) we would arrive to similar
conclusions. Nevertheless, to retain the real and the imaginary parts we need
to impose ${\Re\widetilde{\phi}_{i,1}({\bf r})=\Im\widetilde{\phi}_{i,1}({\bf
r})}$ to remove the triplet contribution.

$\displaystyle\Psi_{+}^{2e-}({\bf r}_{1},{\bf r}_{2})$
$\displaystyle=\frac{1}{\sqrt{2}}\sum_{i}c_{i}\left[\Re\widetilde{\phi}_{i,1}({\bf
r}_{1})\begin{pmatrix}1\\\ a\end{pmatrix}\Re\widetilde{\phi}_{i,1}({\bf
r}_{2})\begin{pmatrix}-a^{*}\\\ 1\end{pmatrix}-\Re\widetilde{\phi}_{i,1}({\bf
r}_{2})\begin{pmatrix}1\\\ a\end{pmatrix}\Re\widetilde{\phi}_{i,1}({\bf
r}_{1})\begin{pmatrix}-a^{*}\\\ 1\end{pmatrix}\right],$ (134)

which in the $a=0$ case allows us to introduce the usual nonrelativistic spin
functions $\alpha=\begin{pmatrix}1\\\ 0\end{pmatrix}$ and
$\beta=\begin{pmatrix}0\\\ 1\end{pmatrix}$ 212121In the $a=1$ case the spin
functions are those of the ${\widehat{S}_{x}}$ operator. On the other hand,
taking $a=\textrm{i}$ we arrive to spin functions of the ${\widehat{S}_{y}}$
operator.. Inserting the spin functions may rewrite the wavefunction in a
simplified notation as

$\Psi_{+}^{2e-}({\bf r}_{1},{\bf
r}_{2})=\frac{1}{\sqrt{2}}\sum_{i}c_{i}\Re\widetilde{\phi}_{i,1}({\bf
r}_{1})\Re\widetilde{\phi}_{i,1}({\bf r}_{2})(\alpha\beta-\beta\alpha)$ (135)

that corresponds to the singlet wavefunction in the nonrelativistic limit with
the $\\{\Re\widetilde{\phi}_{i,1}\\}$ being the scalar NOs. Indeed, from Eq.
(135) the nonrelativistic Fixed-Phases functional can be readily obtained;
therefore, the relativistic Fixed-Phases functional presented in this work
contains its nonrelativistic counterpart.

## References

  * Laerdahl  _et al._ [1998] J. Laerdahl, K. Faegri Jr, L. Visscher, and T. Saue, J. Chem. Phys. 109, 10806 (1998).
  * Tecmer  _et al._ [2011] P. Tecmer, A. S. P. Gomes, U. Ekström, and L. Visscher, Phys. Chem. Chem. Phys. 13, 6249 (2011).
  * Samultsev  _et al._ [2015] D. Samultsev, Y. Y. Rusakov, and L. Krivdin, Russ. Chem. Bull. 64, 2756 (2015).
  * Castro  _et al._ [2019] A. C. Castro, H. Fliegl, M. Cascella, T. Helgaker, M. Repisky, S. Komorovsky, M. Á. Medrano, A. G. Quiroga, and M. Swart, Dalton Trans. 48, 8076 (2019).
  * Tiesinga  _et al._ [2021] E. Tiesinga, J. Kłos, M. Li, A. Petrov, and S. Kotochigova, New J. Phys. 23, 085007 (2021).
  * Valderrama  _et al._ [1997] E. Valderrama, E. V. Ludeña, and J. Hinze, J. Chem. Phys. 106, 9227 (1997).
  * Valderrama  _et al._ [1999] E. Valderrama, E. V. Ludeña, and J. Hinze, J. Chem. Phys. 110, 2343 (1999).
  * Valderrama  _et al._ [2001] E. Valderrama, J. Mercero, and J. Ugalde, J. Phys. B: Atom. Molec. Phys. 34, 275 (2001).
  * Cioslowski [1991] J. Cioslowski, Phys. Rev. A 43, 1223 (1991).
  * Ramos-Cordoba  _et al._ [2016] E. Ramos-Cordoba, P. Salvador, and E. Matito, Phys. Chem. Chem. Phys. 18, 24015 (2016).
  * Benavides-Riveros  _et al._ [2017] C. L. Benavides-Riveros, N. N. Lathiotakis, and M. A. Marques, Phys. Chem. Chem. Phys. 19, 12655 (2017).
  * Ramos-Cordoba and Matito [2017] E. Ramos-Cordoba and E. Matito, J. Chem. Theory Comput. 13, 2705 (2017).
  * Via-Nadal  _et al._ [2017] M. Via-Nadal, M. Rodríguez-Mayorga, and E. Matito, Phys. Rev. A 96, 050501 (2017).
  * Via-Nadal  _et al._ [2019] M. Via-Nadal, M. Rodríguez-Mayorga, E. Ramos-Cordoba, and E. Matito, J. Phys. Chem. Lett. 10, 4032 (2019).
  * Gagliardi and Roos [2000] L. Gagliardi and B. O. Roos, Chem. Phys. Lett. 331, 229 (2000).
  * Gagliardi  _et al._ [2001] L. Gagliardi, B. Schimmelpfennig, L. Maron, U. Wahlgren, and A. Willetts, Chem. Phys. Lett. 344, 207 (2001).
  * Gagliardi and Roos [2005] L. Gagliardi and B. O. Roos, Nature 433, 848 (2005).
  * Gagliardi and Roos [2007] L. Gagliardi and B. O. Roos, Chem. Soc. Rev. 36, 893 (2007).
  * Knecht  _et al._ [2014] S. Knecht, Ö. Legeza, and M. Reiher, J. Chem. Phys. 140, 041101 (2014).
  * Battaglia  _et al._ [2018] S. Battaglia, S. Keller, and S. Knecht, J. Chem. Theory Comput. 14, 2353 (2018).
  * Fleig  _et al._ [2007] T. Fleig, L. K. Sørensen, and J. Olsen, Theor. Chim. Acta (Berlin) 118, 347 (2007).
  * Ghosh  _et al._ [2016] A. Ghosh, R. K. Chaudhuri, and S. Chattopadhyay, J. Chem. Phys. 145, 124303 (2016).
  * Eliav and Kaldor [2010] E. Eliav and U. Kaldor, in  _Recent Progress in Coupled Cluster Methods_ (Springer, 2010) pp. 113–144.
  * Oleynichenko  _et al._ [2020] A. V. Oleynichenko, A. Zaitsevskii, L. V. Skripnikov, and E. Eliav, Symmetry 12, 1101 (2020).
  * Tang  _et al._ [2017] Y.-B. Tang, B.-Q. Lou, and T.-Y. Shi, Phys. Rev. A 96, 022513 (2017).
  * Shee  _et al._ [2018] A. Shee, T. Saue, L. Visscher, and A. Severo Pereira Gomes, J. Chem. Phys. 149, 174113 (2018).
  * Lew-Yee  _et al._ [2021] J. F. H. Lew-Yee, M. Piris, and J. M. del Campo, J. Chem. Phys. 154, 064102 (2021).
  * Wang and Baerends [2022] J. Wang and E. J. Baerends, Phys. Rev. Lett. 128, 013001 (2022).
  * López-Sandoval and Pastor [2002] R. López-Sandoval and G. Pastor, Phys. Rev. B 66, 155118 (2002).
  * Lathiotakis and Marques [2008] N. Lathiotakis and M. A. Marques, J. Chem. Phys. 128, 184103 (2008).
  * Sharma  _et al._ [2013] S. Sharma, J. Dewhurst, S. Shallcross, and E. Gross, Phys. Rev. Lett. 110, 116403 (2013).
  * Pernal [2015] K. Pernal, New J. Phys. 17, 111001 (2015).
  * Mitxelena  _et al._ [2017] I. Mitxelena, M. Piris, and M. Rodríguez-Mayorga, J. Phys. Condens. Matter 29, 425602 (2017).
  * Mitxelena and Piris [2020a] I. Mitxelena and M. Piris, J. Phys. Condens. Matter 32, 17LT01 (2020a).
  * Mitxelena and Piris [2020b] I. Mitxelena and M. Piris, J. Chem. Phys. 152, 064108 (2020b).
  * Cohen  _et al._ [2008] A. J. Cohen, P. Mori-Sánchez, and W. Yang, Science 321, 792 (2008).
  * Su  _et al._ [2018] N. Q. Su, C. Li, and W. Yang, PNAS 115, 9678 (2018).
  * Seidl  _et al._ [2007] M. Seidl, P. Gori-Giorgi, and A. Savin, Phys. Rev. A 75, 042511 (2007).
  * Becke [2013a] A. D. Becke, J. Chem. Phys. 138, 074109 (2013a).
  * Becke [2013b] A. D. Becke, J. Chem. Phys. 138, 161101 (2013b).
  * Schmidt  _et al._ [2019] J. Schmidt, C. L. Benavides-Riveros, and M. A. Marques, Phys. Rev. B 99, 224502 (2019).
  * Benavides-Riveros  _et al._ [2020] C. L. Benavides-Riveros, J. Wolff, M. A. Marques, and C. Schilling, Phys. Rev. Lett. 124, 180603 (2020).
  * Ohsaku  _et al._ [2002] T. Ohsaku, S. Yamanaka, D. Yamaki, and K. Yamaguchi, Int. J. Quant. Chem. 90, 273 (2002).
  * Gilbert [1975] T. L. Gilbert, Phys. Rev. B 12, 2111 (1975).
  * Toulouse [2021] J. Toulouse, SciPost Chem. 1, 002 (2021).
  * Chaix and Iracane [1989] P. Chaix and D. Iracane, Journal of Physics B: Atomic, Molecular and Optical Physics 22, 3791 (1989).
  * Chaix  _et al._ [1989] P. Chaix, D. Iracane, and P.-L. Lions, Journal of Physics B: Atomic, Molecular and Optical Physics 22, 3815 (1989).
  * Mittleman [1981] M. H. Mittleman, Phys. Rev. A 24, 1167 (1981).
  * Saue and Visscher [2003] T. Saue and L. Visscher, in  _Theoretical chemistry and physics of heavy and superheavy elements_ (Springer, 2003) pp. 211–267.
  * Kramers [1930] H. A. Kramers, Proc. Acad. Amst 33 (1930).
  * Kutzelnigg [2000] W. Kutzelnigg, Theor. Chem. Acc. 103, 182 (2000).
  * Kutzelnigg  _et al._ [2010] W. Kutzelnigg, K. Shamasundarb, and D. Mukherjee, Mol. Phys. 108, 433 (2010).
  * Note [1] The spinor-components annihilation field operators are defined as ${\mathaccent 866{\psi}^{\vphantom{\dagger}}}_{\mu}({\bf r})=\sum@\slimits@_{I}{\mathaccent 866{b}^{\vphantom{\dagger}}}_{I}\phi_{I,\mu}({\bf r})+\sum@\slimits@_{R}{\mathaccent 866{d}^{\tmspace+{.2777em}\dagger}}_{R}\phi_{R,\mu}({\bf r})$; the creation field operators are written as ${\mathaccent 866{\psi}^{\tmspace+{.2777em}\dagger}}_{\mu}({\bf r})=\sum@\slimits@_{I}{\mathaccent 866{b}^{\tmspace+{.2777em}\dagger}}_{I}\phi_{I,\mu}({\bf r})+\sum@\slimits@_{R}{\mathaccent 866{d}^{\vphantom{\dagger}}}_{R}\phi_{R,\mu}({\bf r})$.
  * Dirac [1930] P. A. M. Dirac, Proc. R. Soc. Lond. A. 126, 360 (1930).
  * Dyall and Fægri Jr [2007] K. G. Dyall and K. Fægri Jr,  _Introduction to relativistic quantum chemistry_ (Oxford University Press, 2007).
  * Note [2] We exploit the fact that any unitary matrix can be written in terms of an auxiliary antihermitian matrix ${\bm{\kappa}}$.
  * Valatin [1958] J. Valatin, Il Nuovo Cimento (1955-1965) 7, 843 (1958).
  * Bogoljubov [1958] N. Bogoljubov, Il Nuovo Cimento (1955-1965) 7, 794 (1958).
  * Note [3] Also known as the single Slater determinant wavefunction.
  * Helgaker  _et al._ [2000] T. Helgaker, P. Jorgensen, and J. Olsen,  _Molecular Electronic Structure Theory_ (Wiley, Chichester, 2000).
  * Visscher  _et al._ [1991] L. Visscher, P. Aerts, O. Visser, and W. Nieuwpoort, Int. J. Quant. Chem. 40, 131 (1991).
  * Shabaev  _et al._ [2004] V. Shabaev, I. Tupitsyn, V. Yerokhin, G. Plunien, and G. Soff, Physical review letters 93, 130405 (2004).
  * Dyall [2012] K. G. Dyall, Chem. Phys. 395, 35 (2012).
  * Note [4] Assuming that the opposite charge operator ($\widehat{Q}$) and the interacting Hamiltonian operator ($\widehat{H}$) commute, which is only valid for an appropriate definition of the fermion-fermion interaction and with the adequate definition for the normal ordering.
  * Note [5] Could also be call as the natural spinor representation.
  * Levy [1979] M. Levy, Proc. Natl. Acad. Sci. USA 76, 6062 (1979).
  * Lieb [1983] E. H. Lieb, Int. J. Quant. Chem. 24, 243 (1983).
  * Valone [1991] S. M. Valone, Physical Review B 44, 1509 (1991).
  * Note [6] We currently have no easy characterization of the set of $N$-representable 1-RDMs $\mathcal{D}_{Q}$. We at least know that the trace of these 1-RDMs should be equal to the charge $Q$. Other necessary conditions will be investigated in future work.
  * Note [7] Whether we have a true minimum, is an open question, but we suspect it to be the case like in the nonrelativistic limit [67].
  * Note [8] Recall that this vacuum state is given using some orthonormal basis ${\bm{\psi}}_{A}({\bf r})$ that is the initial guess. Also notice that the operators ${\mathaccent 866{b}^{\tmspace+{.2777em}\dagger}}_{I}$, ${\mathaccent 866{d}^{\tmspace+{.2777em}\dagger}}_{R}$, ${\mathaccent 866{b}^{\vphantom{\dagger}}}_{I}$, and ${\mathaccent 866{d}^{\vphantom{\dagger}}}_{R}$ are the initial ones.
  * Lieb and Siedentop [2000] E. H. Lieb and H. Siedentop, Commun. Math. Phys. 213, 673 (2000).
  * Hainzl  _et al._ [2005a] C. Hainzl, M. Lewin, and É. Séré, Commun. Math. Phys. 257, 515 (2005a).
  * Hainzl  _et al._ [2005b] C. Hainzl, M. Lewin, and E. Séré, J. Phys. A Math. Theor. 38, 4483 (2005b).
  * Hainzl  _et al._ [2007a] C. Hainzl, M. Lewin, and J. P. Solovej, Comm. Pure Appl. Math. 60, 546 (2007a).
  * Hainzl  _et al._ [2007b] C. Hainzl, M. Lewin, E. Séré, and J. P. Solovej, Phys. Rev. A 76, 052104 (2007b).
  * Lewin [2011] M. Lewin, in  _Mathematical Results In Quantum Physics: (With DVD-ROM)_ (World Scientific, 2011) pp. 45–59.
  * Gravejat  _et al._ [2011] P. Gravejat, M. Lewin, and É. Séré, Commun. Math. Phys. 306, 1 (2011).
  * Brouder [2002] C. Brouder, EPJ direct 4, 1 (2002).
  * Sucher [1980] J. Sucher, Phys. Rev. A 22, 348 (1980).
  * Hubbell [2006] J. H. Hubbell, Radiat. Phys. Chem. 75, 614 (2006).
  * Note [9] Let us mention that NS are not populated within the np approximation, but they cannot be neglected when optimizing the PS (i.e. NS play a similar role as the secondary space in nonrelativistic complete active space calculations).
  * Garrod and Percus [1964] C. Garrod and J. K. Percus, J. Math. Phys. 5, 1756 (1964).
  * Herbert and Harriman [2003] J. M. Herbert and J. E. Harriman, J. Chem. Phys. 118, 10835 (2003).
  * Visscher  _et al._ [1994] L. Visscher, O. Visser, P. J. Aerts, H. Merenga, and W. Nieuwpoort, Comput. Phys. Commun. 81, 120 (1994).
  * Gaunt [1929] J. A. Gaunt, Proc. R. Soc. A 112, 513 (1929).
  * Löwdin and Shull [1956] P.-O. Löwdin and H. Shull, Phys. Rev. 101, 1730 (1956).
  * Piris [2013] M. Piris, Int. J. Quant. Chem. 113, 620 (2013).
  * Piris [2014] M. Piris, J. Chem. Phys. 141, 044107 (2014).
  * Piris [2018] M. Piris, Phys. Rev. A 98, 022504 (2018).
  * Piris [1999] M. Piris, J. Math. Chem. 25, 47 (1999).
  * Note [10] Let us comment that the formation of pairs of degenerated spinor is not unique, but we have adopted the one given by Kramers’ pairs in this work.
  * Note [11] The same kind of partitioning of the spinor states can be applied to the NS.
  * Note [12] The uppercase indices that enter in the definition of ${\mathaccent 866{x}^{s_{1}s_{2}}_{IK,JL}}$ are readily used in the creation and annihilation operators, but they are used as (lowercase) unbarred indices in $T_{ij}$ operator. Actually, the $T_{ij}$ is only defined with unbarred indices.
  * Piris [2017] M. Piris, Phys. Rev. Lett. 119, 063002 (2017).
  * Rodríguez-Mayorga  _et al._ [2021] M. Rodríguez-Mayorga, I. Mitxelena, F. Bruneval, and M. Piris, J. Chem. Theory Comput. 17, 7562 (2021).
  * Mentel  _et al._ [2014] L. M. Mentel, R. van Meer, O. V. Gritsenko, and E. J. Baerends, J. Chem. Phys. 140, 214105 (2014).
  * Rodríguez-Mayorga  _et al._ [2017] M. Rodríguez-Mayorga, E. Ramos-Cordoba, M. Via-Nadal, M. Piris, and E. Matito, Phys. Chem. Chem. Phys. 19, 24029 (2017).
  * Sun  _et al._ [2021] S. Sun, T. F. Stetina, T. Zhang, H. Hu, E. F. Valeev, Q. Sun, and X. Li, J. Chem. Theory Comput. 17, 3388 (2021).
  * Hafner [1980] P. Hafner, J. Phys. B 13, 3297 (1980).
  * Müller [1984] A. M. K. Müller, Phys. Lett. 105A, 446 (1984).
  * Buijse [1991] M. A. Buijse,  _Thesis: Electron Correlation. Fermi and Coulomb holes, dynamical and nondynamical correlation_ , Ph.D. thesis, Vrije Universiteit, Amsterdam, The Netherlands (1991).
  * Buijse and Baerends [2002] M. A. Buijse and E. J. Baerends, Mol. Phys. 100, 401 (2002).
  * Cioslowski and Pernal [1999] J. Cioslowski and K. Pernal, J. Chem. Phys. 111, 3396 (1999).
  * Sharma  _et al._ [2008] S. Sharma, J. K. Dewhurst, N. N. Lathiotakis, and E. K. U. Gross, Phys. Rev. B 78, 201103 (2008).
  * Cioslowski and Pernal [2000] J. Cioslowski and K. Pernal, Phys. Rev. A 61, 034503 (2000).
  * Goedecker and Umrigar [1998] S. Goedecker and C. J. Umrigar, Phys. Rev. Lett. 81, 866 (1998).
  * Marques and Lathiotakis [2008] M. A. L. Marques and N. N. Lathiotakis, Phys. Rev. A 77, 032509 (2008).
  * Lieb and Oxford [1981] E. H. Lieb and S. Oxford, Int. J. Quant. Chem. 19, 427 (1981).
  * Kamil  _et al._ [2016] E. Kamil, R. Schade, T. Pruschke, and P. E. Blöchl, Phys. Rev. B 93, 085141 (2016).
  * Löwdin [1955] P.-O. Löwdin, Phys. Rev. 97, 1474 (1955).
  * Piris  _et al._ [2013] M. Piris, J. Matxain, and X. Lopez, J. Chem. Phys. 139, 234109 (2013).
  * Mitxelena  _et al._ [2018] I. Mitxelena, M. Rodríguez-Mayorga, and M. Piris, Eur. Phys. J. B 91, 109 (2018).
  * Note [13] A variant proposed for PNOF7 functional that is used to define NOF-MP2 method [90].
  * Note [14] The relativistic version of PNOF6 functional is not included because one of us proved that beyond the perfect-pairing approach this functional does not obey the sum rule [98].
  * Piris [2021] M. Piris, Phys. Rev. Lett. 127, 233001 (2021).
  * Shull and Löwdin [1959] H. Shull and P.-O. Löwdin, J. Chem. Phys. 30, 617 (1959).
  * Pernal [2013] K. Pernal, Comput. Theor. Chem. 1003, 127 (2013).
  * Coleman [1963] A. J. Coleman, Rev. Mod. Phys. 35, 668 (1963).
  * Smith [1966] D. W. Smith, Phys. Rev. 147, 896 (1966).
  * Borland and Dennis [1972] R. Borland and K. Dennis, J. Phys. B 5, 7 (1972).
  * Coleman [1978] A. Coleman, Int. J. Quant. Chem. 13, 67 (1978).
  * Klyachko [2006] A. A. Klyachko, in  _Journal of Physics: Conference Series_ , Vol. 36 (IOP Publishing, 2006) p. 72.
  * Altunbulak and Klyachko [2008] M. Altunbulak and A. Klyachko, Commun. Math. Phys. 282, 287 (2008).
  * Mazziotti [2000] D. A. Mazziotti, in  _Many-Electron Densities and Reduced Density Matrices_ (Springer, 2000) pp. 139–163.
  * Schilling  _et al._ [2018] C. Schilling, M. Altunbulak, S. Knecht, A. Lopes, J. D. Whitfield, M. Christandl, D. Gross, and M. Reiher, Phys. Rev. A 97, 052503 (2018).
  * Theophilou  _et al._ [2015] I. Theophilou, N. N. Lathiotakis, M. A. L. Marques, and N. Helbig, J. Chem. Phys. 142, 154108 (2015).
  * Chakraborty and Mazziotti [2014] R. Chakraborty and D. A. Mazziotti, Phys. Rev. A 89, 042505 (2014).
  * Chakraborty and Mazziotti [2016] R. Chakraborty and D. A. Mazziotti, Int. J. Quant. Chem. 116, 784 (2016).
  * Mazziotti [2016] D. A. Mazziotti, Phys. Rev. A 94, 032516 (2016).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Pino  _et al._ [2007] R. Pino, O. Bokanowski, E. V. Ludeña, and R. L. Boada, Theor. Chim. Acta (Berlin) 118, 557 (2007).
  * Penz and van Leeuwen [2021] M. Penz and R. van Leeuwen, arXiv preprint arXiv:2106.15370 (2021).
  * Note [15] Notice that the $\lambda_{\bar{i}\bar{j},\bar{k}\bar{l}}$ terms are obtained by replacing the unbarred indices by barred ones.
  * Note [16] For $i=j$, the elements are 0 by definition and they fulfill the condition.
  * Piris  _et al._ [2011] M. Piris, X. Lopez, F. Ruipérez, J. M. Matxain, and J. M. Ugalde, J. Chem. Phys. 134, 164102 (2011).
  * Note [17] Recall that the ‘spin’ information is contained in the 4-component spinors.
  * Carlson and Keller [1961] B. Carlson and J. M. Keller, Phys. Rev. 121, 659 (1961).
  * Schmidt [1907] E. Schmidt, Math. Ann. 63, 433 (1907).
  * Note [18] ${\mathaccent 866{\mathcal{K}}}$ in the many electron case leads to the complex conjugate of the CI coefficients and transforms all spinors in the tensor product [141, 55].
  * Jørgen A. Jensen  _et al._ [1996] H. Jørgen A. Jensen, K. G. Dyall, T. Saue, and K. Fægri Jr, J. Chem. Phys. 104, 4083 (1996).
  * Note [19] This is normally the case in the nonrelativistic limit [55].
  * Note [20] Retaining only the imaginary part of the scalar spinors (setting ${\Re\mathaccent 869{\phi}_{i,1}({\bf r})=0}$) we would arrive to similar conclusions. Nevertheless, to retain the real and the imaginary parts we need to impose ${\Re\mathaccent 869{\phi}_{i,1}({\bf r})=\Im\mathaccent 869{\phi}_{i,1}({\bf r})}$ to remove the triplet contribution.
  * Note [21] In the $a=1$ case the spin functions are those of the ${\mathaccent 866{S}_{x}}$ operator. On the other hand, taking $a=\textrm{i}$ we arrive to spin functions of the ${\mathaccent 866{S}_{y}}$ operator.

Figure 2: Graphical abstract.