11institutetext: University of California, San Diego
11email<EMAIL_ADDRESS>
# Towards Explainable Automated Neuroanatomy

Kui Qian1,∗ Anonymous Organizations
11email<EMAIL_ADDRESS>Litao Qiao1 Anonymous Organizations
11email<EMAIL_ADDRESS>Beth Friedman2 Anonymous Organizations
11email<EMAIL_ADDRESS>Edward O’Donnell3 Anonymous Organizations
11email<EMAIL_ADDRESS>
David Kleinfeld3,4 and Yoav Freund2,5,∗ Anonymous Organizations
11email<EMAIL_ADDRESS>
###### Abstract

We present a novel method for quantifying the microscopic structure of brain
tissue. It is based on the automated recognition of interpretable features
obtained by analyzing the shapes of cells. This contrasts with prevailing
methods of brain anatomical analysis in two ways. First, contemporary methods
use gray-scale values derived from smoothed version of the anatomical images,
which dissipated valuable information from the texture of the images. Second,
contemporary analysis uses the output of black-box Convolutional Neural
Networks, while our system makes decisions based on interpretable features
obtained by analyzing the shapes of individual cells. An important benefit of
this open-box approach is that the anatomist can understand and correct the
decisions made by the computer. Our proposed system can accurately localize
and identify existing brain structures. This can be used to align and
coregistar brains and will facilitate connectomic studies for reverse
engineering of brain circuitry.

###### Keywords:

Explainable ML Brain texture Computational anatomy

1 Department of Electrical and Computer Engineering, University of California,
San Diego, La Jolla, CA 92093, USA
2 Department of Computer Science and Engineering, University of California,
San Diego, La Jolla, CA 92093, USA
3 Department of Physics, University of California, San Diego, La Jolla, CA
92093, USA
4 Department of Neurobiology, University of California, San Diego, La Jolla,
CA 92093, USA
5 Halıcıoğlu Data Science Institute, University of California, San Diego, La
Jolla, CA 92093, USA
∗ Correspondence<EMAIL_ADDRESS><EMAIL_ADDRESS>

##  1 Introduction

One of the first steps in brain analysis is to answer the “where” question. To
answer this question the anatomist typically relies on brain cytoarchitecture,
namely, the spatial organization of neural elements. However, manual labeling
of the brain structures is a labor-intensive task. Typically, identifying and
marking the boundaries of 40 standard landmarks in a single brain takes a
trained anatomist many weeks of work.

Based on [6], our new work introduces a machine learning-based approach to
automate the identification of structures within high-resolution mouse brain
images (Figure 1). The method in [6] introduced the use of high resolution
texture of brain tissue to identify different brain regions. The process of
identification relied on a pre-trained convolutional neural network (CNN) that
processes pixel-level features. This “black box” approach leads to inherently
uninterpretable results. Our new system utilizes interpretable cell shape
features for structure detection. Our underlying assumption is that the
distribution of the cell shapes inside and outside the structures should be
significantly different, which is analogous to the criteria used by anatomists
for structure identification [4]. By focusing on individual cells as the
primary unit of analysis, our method not only makes self-explainable decisions
for the anatomists but also maintains robustness across different staining
procedures and imaging techniques. We demonstrate that the features we compute
for cell shapes can be used as inputs to structure detectors of different
brain regions.

Of interest, a recent independent study [9] also uses cell shapes to build an
interpretable machine learning method for cortical cytoarchitecture analysis.
This study solved a classification challenge of detection of human cortex
laminae using a neuron-centric approach. Here we apply a neuron-centric
approach to detect the more general case of non-laminated geometrically
diverse spatial distributions of neurons. Such distributions are particularly
prevalent in the mouse brainstem, a region that is typically very challenging
to map. Our method relies on two innovative approaches. Firstly, on top of
manually designed features, our approach utilizes unsupervised learning to
extract nuanced and interpretable features that more accurately represent the
diversity of cell shapes. By not relying solely on predefined features, our
system can uncover patterns that might otherwise be overlooked. Secondly, we
incorporate regional features that summarize the statistical attributes of
cell populations, rather than individual cells alone. These features capture
the collective properties of cell groups, such as density and orientation
distributions, to provide a robust framework for our detector. This approach
enables the system to recognize and classify structures even when individual
cell shapes are ambiguous or when cells exhibit subtle differences that are
only discernible as a population. The integration of these advanced regional
features significantly enhances the detector’s robustness, making it more
resilient to variations in staining methods and imaging modalities.

##  2 Methods

Figure 1:  The pipelines of our proposed system, where the figures indicate
the outputs of the stages and Roman letters enumerate the procedures. I:
OpenCV for cell segmentations. II: K-means, Diffusion Mapping, and feature
space alignment for extracting cell features. III: Cumulative CDFs for
generating regional features. IV: XGBoost for the classification of the
region.

###  2.0.1 Quantifying cell shapes

Since our method is based on identifying the shapes of the cells inside and
outside the interested structure, the first part of our method is to generate
cell features to quantify the shapes of all cells in a single brain. Our
method starts with cell segmentations from all brain images using OpenCV [3],
a library of computer vision tools that includes adaptive thresholding and
connected components to isolate cell images from the noisy images. Then all
cell images are zero-padded to create pre-defined uniformly sized cell patches
for efficient handling of the subsequent processing steps. To efficiently
represent each cell patch, we employ a dimensionality reduction technique via
Diffusion Mapping (DM) [2, 7]. This creates a continuous non-linear mapping
from single-cell patches into $m$-dimensional feature vectors that represent
the $m$ most significant eigenvectors of the Laplace-Beltrami operator. In our
case, we choose $m=10$ for all brains, which significantly reduces the
dimensionality from the vast number of original pixel values to a more
manageable feature space that best quantifies the cell shapes. However, as we
typically extract tens of millions of cells from each brain, the dataset used
to train DM is too large to fit in the computer memory. We thus used a
streaming implementation of the K-means algorithm [1] as an efficient way to
create a small number of representative cell patches as the training set for
the DM algorithm (supplementary materials).

Conceptually, DM transforms cell patches into a 10-dimensional feature space,
which allows the visualization of these patches by treating their features as
coordinates within this space. By projecting cell patches in the same brain
onto a 2D plane from this space, we create a “patch cloud” that represents
cell shape diversity along the chosen two dimensions. The visualization of the
patch clouds shows that the cell features learned through DM capture the cell
shapes in a visually explainable way. Furthermore, after visually analyzing
the patch clouds generated from different brains, we found that these patch
clouds share similar shapes, despite the fact that the brain images are
obtained using different stains (thionin vs. NeuroTrace blue) and imaging
techniques (brightfield vs. fluorescent). In particular, the clouds from the
different brains only differ by the orders of the axes or the orientations of
the shapes when they are visualized in the 2D planes. This motivates the idea
that the cell features from any brain can be aligned to a fixed set of
features using a simple affine transformation. By adopting the cell features
from the selected brain images as the reference, we formulated a root mean
squared (RMS) optimization procedure (supplementary materials) that has a
closed-form solution to formulate the affine transformation matrix. For
example, Figure 2 shows the two patch clouds projected onto the subspace
formed by the 1st and 4th eigenvectors. Each cloud consists of the cell
patches obtained using the K-means algorithm that are representatives of the
cells in two different brains that are stained using thionin and NeuroTrace
blue, respectively. The patch clouds are visualized by projection onto the
subspace formed by the 1st and 4th eigenvectors. The clouds overlap almost
perfectly (Figure 2a), which confirms the effectiveness of the affine
transformation for aligning cell features from different brains. Further, the
cell patches with similar shapes are grouped together in different regions
(Figure 2b–d), which implies that the 4th eigenvector on the x-axis describes
the aspect ratio of the cell shape and the 1st eigenvector on the y-axis gives
the sizes of the cells.

Figure 2:  (a) Visualization of patch clouds from two brains with different
modalities on the 1st and 4th fixed features. Each “dot” in the image is a
cell patch. The red patches are from a brain that is stained with thionin and
imaged using brightfield, and the blue patches are from a brain that is
stained with NeuroTrace blue and imaged using fluorescence. (b) The region of
very small cells. (c) The region of large and round cells. (d) The region of
thin cells.

In addition to the 10 cell features identified through the DM technique, we
incorporate an additional 10 manually designed features (supplementary
materials) that explicitly describe the shape of the cells. These
supplementary features are directly extracted from the images, without
undergoing a learning process. Finally, our feature vector for each cell
encompasses a total of 20 attributes. To facilitate the extraction of cell
features for groups of cells in subsequent steps of our method, we will
establish a database that stores the feature vectors for all cells within the
brain. It is important to highlight that this feature extraction process is
based on unsupervised learning algorithms. This can be performed even without
manual annotations of the structures by experts in anatomy.

###  2.0.2 Structure detection

One challenge of using the neuron-centric approach for the task of structure
detection lies in the fact that the single-cell shape does not possess the
distinctive properties to be classified as either inside or outside a target
structure. Since cytoarchitecture is a property of collections of cells, a
well-trained anatomist usually determines the boundary of the structure based
on the distribution of a group of cells in a small region. This insight
motivates us to develop our structure detection system to classify image
regions rather than individual cells, which are characterized using regional
features that represent the cell shape distribution in a specific region.
Also, analyzing regions instead of individual cells can greatly reduce the
sensitivity of the analysis to segmentation errors. By leveraging the 20 cell
features designed for characterizing individual cell shapes, we can describe a
region that contains a collection of cells through the distributions of these
features, which are represented by their empirical cumulative distribution
functions (CDFs). Therefore, the two regions have different characteristic
cell shapes if their CDFs are different (Figure 3).

Figure 3:  Cell shape distributions for two regions near the Abducens Nucleus
(6N). The left image highlights a region that roughly corresponds to the 6N
structure (circled in green), while the right image shows a region outside the
structure (circled in red). The CDF graph represents the rotation feature of
cells, with the green curve depicting the CDF for cells within region 1 and
the red curve for those within region 2.

To compute the region feature for a region containing a group of cells, we
first query the cell feature database to retrieve features for all cells
within the region, followed by generating cumulative CDFs for these cell
features. Each region’s feature CDF is then discretized into 99 points,
resulting in a comprehensive region feature vector of $20\times 99=1980$
dimensions. Then we append two additional features: cell density per unit area
and the area ratio covered by cells, which results in a final region feature
vector with 1982 elements. Finally, our structure detection model takes the
vector of the regional features and predicts the likelihood of the region
belonging to a specific structure. Our model employs XGBoost [5], which is a
supervised learning algorithm that is particularly suitable for our case.
XGBoost not only inherits Adaboost’s resistance to overfitting but also is
interpretable in the sense that feature importance can be readily derived from
the trained model which is critical for anatomists to understand the decision
of the model. Given that all our brains have 26 annotated structures, we opted
for 26 binary classifiers rather than a single multi-class classifier, as this
allows for more precise and focused detection of each specific structure.

##  3 Results

###  3.0.1 Classifier accuracy

We evaluate our method by comparing the predictive performance of our
structure detector with the CNN texture classifier introduced in [6], both
quantitively and qualitatively. The quantitive evaluation of our method was
conducted using the same three human-annotated brain image stacks captured
from sagittal views, stained with thionin, and scanned at a resolution of
$0.5\mu m$ using brightfield imaging [6]. For the CNN texture classifiers, we
utilized the pre-existing models from [6] without further training. Unlike the
CNN approach, whose input image dimension is limited to $224\times 224$
pixels, our method can process regions of any shape. However, to ensure a fair
comparison, we adopted the same preprocessing procedures as [6] to prepare the
training set for our method and the test set for both methods. In particular,
all images from the three brains were split into small image tiles of size
$224\times 224$ pixels using the sliding window technique, which will be used
as the input regions for our method. Given that we are training 26 binary
classifiers corresponding to each of the 26 structures, a training image tile
is labeled as positive for a specific classifier if more than half of the tile
is within the anatomists’ defined structure boundary and is labeled as
negative for that classifier otherwise. The image tiles from two of the three
brain images were used as the training set, while the third one was used as
the test set for performance evaluation.

The performance of both methods was assessed using the Area Under the Receiver
Operating Characteristic Curve (ROC AUC) metric. The ROC curves for six chosen
structures using our method are shown in Figure 4a, while Figure 4b displays
the scores of both methods across all 26 structures. The ROC AUC scores for
our method consistently demonstrate high predictive performance, with the
lowest scores surpassing 0.75 and the highest nearing 1.0. Thus our method
reliably differentiates between the various structures. The average ROC AUC
score for our method is 0.892, which, while slightly lower than the CNN’s
average of 0.924, still represents a robust performance.

Figure 4:  Comparison of ROC AUC scores for detecting different structures.
Please see [6] for the list of all 26 structures and their abbreviations.

###  3.0.2 Probability map

The unique advantage of our method is that our model can be directly applied
to images that have different textures from the training images without
retraining. In the following experiment, we used the same classifiers used in
the first experiment, which were trained on thionin-stained brain images, to
produce probability maps for brain images derived from a brain stained with
NeuroTrace blue. This staining method yields brain images with pixel
intensities markedly different from those stained with thionin. We followed
the same step used in [6] to generate the probability maps for both methods,
wherein the output of the classifier determines the probability of each image
tile. Because of the dearth of annotations of structures for this brain, we
will evaluate the probability maps from both methods qualitatively through
visual inspection. In Figure 5a, we show 7 different image sections of the
left locus coeruleus (LC), where we can clearly see the structure in each
image section defined by the black and grey cells. Our method yields high-
probability regions that are closely aligned with the actual distribution of
cells, as shown in the original image patches (Figure 5b). This alignment is
visually more precise in the probability maps of our method than in the CNN
probability maps (Figure 5c). Note the regions of high probability derived
using the CNN method do not accurately overlay the cell-marked areas. This
observation underscores the robustness of our method in generating structure-
specific probability maps across different staining modalities, ensuring its
utility without the necessity of retraining the model for each new staining
technique.

Figure 5:  Comparison of structure identification across different sections
using probability maps. (a) Original image patches of the LC with cells
delineated in black. (b) Probability maps for each of the sections of our
method. (c) Probability maps of the CNN method. Bright and dark colors
represent high and low probabilities of being inside the structure,
respectively.

###  3.0.3 Understanding the detector

Utilizing XGBoost as our classification algorithm offers the distinct
advantage of providing feature importance metrics for each input variable.
This is particularly useful given that our cell and regional features are
intrinsically interpretable, which allows us to give anatomists a visual
explanation for the model’s decisions. Specifically, we can trace back to the
original images and highlight cells that satisfy a feature deemed critical by
XGBoost in identifying a particular structure.

We now consider explanations of our automated method. In Figure 6a the green
outline labels the superior colliculus (SC) as identified by anatomists. In
contrast, the areas our detector identifies as SC are shaded in red, which is
suboptimal as there is a large area of the unrecognized region (blue color)
inside the structure. Intriguingly, this behavior is understandable as the
discernible boundary marked by the black line made by our detector aligns with
known biological layers within SC: the superficial (SG) and deep (DG) grey [8]
(Figure 6b). We can understand why our method makes such decisions by
highlighting the cells that satisfy the features deemed important by the
XGBoost model. In this case, the XGBoost model put the highest importance on
the feature “rotation–11.3”, which represents the cumulative probability that
a cell’s rotation angle is less than $11.3\degree$. By highlighting the cells
whose rotation angles are in the range between $-65\degree$ and $11.3\degree$
in dark red, we find that the highlighted cells predominantly populate the red
regions, while being scarcely present in the blue regions (Figure 6c). Thus, a
biological rationale for our detector’s partitioning is primarily based on the
density of the cells oriented from the lower left to the upper right. The
disparities in the empirical CDFs for the rotation angle of cells illustrate a
clear statistical distinction between the cell orientations in the regions
identified with high ($>1$) and low ($<-1$) XGBoost scores (Figure 6e). It is
worth noting that we deliberately omit the cells with angles less than
$-65\degree$ for visual clarity, since they contribute very little to the
XGboost decision process because of the similar densities of such cells in
both highlighted regions; see CDFs in Figure 6d.

Figure 6:  Explanation of the SC structure detection. (a) The result of our
detector (green contour), where there is a clear boundary between the
recognized regions (red) and unrecognized regions (blue). (b) The SG and DG
gray layers within the SC, which roughly correspond to the red and blue
regions predicted by our detector. (c) and (d) A closer view of the cellular
orientation that our detector uses for its classification shows a
concentration of cells oriented from the lower left to the upper right in the
red region, in contrast to the sparsity of such cells in the blue region. (e)
The CDFs of the cell rotation angles for cells in high-score and low-score
regions, where the cells with angles between $-65\degree$ and $11.3\degree$
are highlighted in (c) for a visual explanation.

##  4 Conclusion

Our method introduces significant contributions to computational neuroanatomy:
(1) an unsupervised learning procedure for efficiently extracting quantifiable
cell shape features, (2) region features that encapsulate the statistical
properties of cell populations, and (3) the deployment of supervised learning,
specifically XGBoost, for robust structure detection. By focusing on cell
shapes and statistical properties, our approach achieves high accuracy and
interpretability, aligning closely with the criteria used by anatomists.
Future directions include refining segmentation techniques, exploring advanced
unsupervised learning models for feature extraction, extending the method’s
applicability to other species or brain regions, and integrating with
neuroanatomical databases to uncover new brain structures.

## References

  * [1] Arthur, D., Vassilvitskii, S., et al.: k-means++: The advantages of careful seeding. In: Soda. vol. 7, pp. 1027–1035 (2007)
  * [2] Belkin, M., Niyogi, P.: Laplacian eigenmaps for dimensionality reduction and data representation. Neural Computation 15(6), 1373–1396 (2003)
  * [3] Bradski, G.: The OpenCV Library. Dr. Dobb’s Journal of Software Tools (2000)
  * [4] Braitenberg, P.D.V.: On the texture of brains. In: Heidelberg Science Library (1977), https://api.semanticscholar.org/CorpusID:10281267
  * [5] Chen, T., Guestrin, C.: Xgboost: a scalable tree boosting system. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. pp. 785–794 (2016), https://github.com/dmlc/xgboost
  * [6] Chen, Y., McElvain, L.E., Tolpygo, A.S., Ferrante, D., Friedman, B., Mitra, P.P., Karten, H.J., Freund, Y., Kleinfeld, D.: An active texture-based digital atlas enables automated mapping of structures and markers across brains. Nature Methods 16(4), 341–350 (2019)
  * [7] Coifman, R.R., Lafon, S., Lee, A.B., Maggioni, M., Nadler, B., Warner, F., Zucker, S.W.: Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps. Proceedings of the National Academy of Sciences 102(21), 7426–7431 (2005)
  * [8] Franklin, K.B., Paxinos, G.: Paxinos and Franklin’s the Mouse Brain in Stereotaxic Coordinates, Compact: The Coronal Plates and Diagrams. Academic Press (2019)
  * [9] Štajduhar, A., Lipić, T., Lončarić, S., Judaš, M., Sedmak, G.: Interpretable machine learning approach for neuron-centric analysis of human cortical cytoarchitecture. Scientific Reports 13 (2023), https://api.semanticscholar.org/CorpusID:257928559

Supplementary Materials for “Towards Explainable Automated Neuroanatomy”
Anonymous

###  4.0.1 Details about OpenCV procedure

For cell segmentation from brain section images, we used the adaptive
thresholding technique provided by OpenCV. This process, specifically through
the cv.adaptiveThreshold function, computes the threshold for a pixel based on
a small region around it. The parameters used in this function were as
follows:

        thresholdType: cv2.THRESH_BINARY,
        adaptiveMethod: cv.ADAPTIVE_THRESH_GAUSSIAN_C,
        blockSize: 101 (chosen based on the typical cell size),
        C: -12 (selected manually based on segmentation quality).

###  4.0.2 Details about KMeans

We applied the Kmeans clustering algorithm on an extensive dataset of ten
million cell patches extracted from a single brain’s segmentation to find
around one thousand representative cell patches. The Kmeans procedure included
two key steps: initialization with Kmeans++ and subsequent refinement.
Kmeans++ was used to select 2,000 initial cluster centroids. Then these
centroids were used to aggregate similar cell patches into clusters. The mean
of each cluster was computed as the final set of representative cell patches.
Clusters containing fewer than 5 samples were excluded. Figure 7 shows samples
of the original cell images and of representative cell patches selected using
K-means after normalizing rotation.

###  4.0.3 Details about Diffusion map

We used the public implementation of DM
111https://github.com/DiffusionMapsAcademics/pyDiffMap. The parameters used in
this function were as follows:

        n_evecs  Number of diffusion map eigenvectors: 100,
        k  Number of nearest neighbors to construct the kernel: 100,
        epsilon: 5000 (chosen based on the typical cell size),
        alpha: 1.0, neighbor_params: {’n_jobs’: -1, ’algorithm’: ’ball_tree’}.

###  4.0.4 Details about RMS-based optimization

The mathematical formulation of finding a linear transformation between two
different diffusion maps can be described as follows: Given a set of vector
pairs $(a_{1},b_{1}),\ldots,(a_{n},b_{n})$ where each of the vectors
$a_{i},b_{i}$ are in a $d$ dimensional space $R^{d}$, find an offset vector
$\mu\in r^{d}$ and a linear transformation $M$ which is a $d\times d$ matrix
so that the following cost function is minimized:

$cost(\mu,M)=\frac{1}{n}\sum_{i=1}^{n}||\mu+Ma_{i}-b_{i}||_{2}^{2}$ (1)

We first compute the hessian matrix of 1 with respect to $\mu$ and $M$ to see
if the problem can be directly solved.

$\frac{\partial^{2}cost(\mu,M)}{\partial\mu\partial\mu^{T}}=2I\rightarrow p.d$
(2) $\frac{\partial^{2}cost(\mu,M)}{\partial M\partial
M^{T}}=\frac{2}{n}\sum_{i=1}^{n}a_{i}a_{i}^{T}\rightarrow p.s.d$ (3)

The results show that both hessian matrices are Positive Semi-definite
matrices and thus we can find the solution by setting derivatives to 0. We
have:

$(\frac{1}{n}\sum_{i=1}^{n}a_{i}a_{i}^{T}-\frac{1}{n}\sum_{i=1}^{n}a_{i}\times\frac{1}{n}\sum_{i=1}^{n}a_{i}^{T})M^{T}=\frac{1}{n}\sum_{i=1}^{n}a_{i}b_{i}^{T}-\frac{1}{n}\sum_{i=1}^{n}a_{i}\times\frac{1}{n}\sum_{i=1}^{n}b_{i}^{T}$
(4)

###  4.0.5 Details about 10 manually designed features

The 10 manually designed features encompass the following: width, height, and
area of a cell; rotation angle along with its confidence level; mean and
standard deviation of pixel intensities within the cell image; the size of the
cell patch, and standard deviation of horizontal and vertical coordinates of
all cell pixels in the cell patch.

###  4.0.6 Details about XGBoost classifiers

The parameters used to train our XGBoost classifiers were as follows:

        params: {’max_depth’:3, ’eta’: 0.2,
                ’objective’:’binary:logistic’, ’num_class’:1},
        num_boost_round: 100.

Figure 7:  (a) Original cell patches, (b) Cell patches selected using K-means
after normalizing rotation.


11affiliationtext: Department of Computer Science and Engineering, Ghent
University, Ghent, Belgium22affiliationtext: QBee B.V. Quantum Accelerated,
Leuven, Belgium

# References

Kourosh Sayar Dogahe  Tamara Sarac  Delphine De Smedt  Koen Bertels

# "Toward" Metal-Organic Framework Design by Quantum Computing

Kourosh Sayar Dogahe  Tamara Sarac  Delphine De Smedt  Koen Bertels

Abstract

The article summarizes the study performed in the context of the Deloitte
Quantum Climate Challenge in 2023 111Deloitte hosts an annual Quantum Climate
Challenge aiming to encourage climate-relevant collaborations between
sustainability and quantum computing experts. The QBee team won first place in
the 2023 challenge among 118 registrations from 33 countries. We want to
express our gratitude to Deloitte and IBM Quantum companies for their
assistance in providing the necessary materials. Link: [ Delloite Challenge
’23]. 222We would also like to thank QBee Company for their generous support
and for helping us with our Quantum Computing research and development needs.
To learn more about QBee Company, please visit their website at: [ QBee.eu]..
We present a hybrid quantum-classical method for calculating Potential Energy
Surface scans, which are essential for designing Metal-Organic Frameworks for
Direct Air Capture applications. The primary objective of this challenge was
to highlight the potential advantages of employing quantum computing. To
evaluate the performance of the model, we conducted total energy calculations
using various computing frameworks and methods. The results demonstrate, at a
small scale, the potential advantage of quantum computing-based models. We
aimed to define relevant classical computing model references for method
benchmarking. The most important benefits of using the PISQ approach for
hybrid quantum-classical computational model development and assessment are
demonstrated.

Keywords: quantum computing; quantum computational chemistry; PISQ; molecular
simulation; VQE.

##  1 Introduction

With the anticipated progress in quantum technology, quantum computing could
bring about significant advantages over classical computing in the future.
These include the ability to tackle previously intractable problems and the
potential for substantial computational speedup. Computational chemistry is a
domain anticipated as among the first ones to potentially benefit from quantum
computing. It tackles the behavior of electrons and atoms at the quantum
level, and classical computers struggle to accurately simulate those systems.
Therefore, it is crucial to develop quantum technology-based computational
methods, also in the domain of computational chemistry.

Quantum computing is a multidisciplinary domain, that combines hardware and
software development. Extensive work on quantum computing in the last decade
led to the establishment of approaches like the NISQ (Noisy Intermediate-Scale
Quantum), which refers to a category of near-term quantum computers handling
noisy qubits, resulting from different noise sources such as decoherence, gate
errors, and measurement errors [1]. Going further, there is the PISQ (Perfect
Intermediate-Scale Quantum) approach that enables various scientific and
industrial communities to step into the quantum computing field by letting the
hardware concerns aside and focusing on the quantum computing logic for their
specific expert domains [2]. The "perfect qubit" is a theoretical concept,
defined after the famous talk by Richard Feynman [3], that is based on perfect
qubit behavior simulation on a classical computer. Perfect qubits are often
used to explore the full potential of quantum computing algorithms and
protocols. PISQ and NISQ represent complementary strategies that are expected
to converge to an operational field within the next 10-15 years. [1, 2]

In the framework of Deloitte’s Quantum Climate Challenge 2023, an interesting
use case was proposed. The primary goal was to investigate how quantum
computers may help to improve materials used in Direct Air Capture (DAC) of
$\mathrm{CO_{2}}$. We approached the challenge from the chemistry domain
following the PISQ approach, with a goal to assess the structure of metal-
organic frameworks (MOFs), with particular attention to MOF74 [4].

The main observation indicates that a potential advantage of quantum
computing-based models over classical methods can be demonstrated at a small
scale. Namely, quantum computation based results reach the accuracy of the
methods that incorporate an analytical approach, going beyond the accuracy
obtained by conventional classical methods employment. Our attempt to scale up
the quantum-based computation revealed some drawbacks of present quantum
models and tools.

The paper is structured as follows. We first present the goal of the challenge
and our approach to it in the Problem Definition Section. In the Method
section, we provide detailed descriptions of all the methods and
implementations. This is followed by a Results and Discussion section where
the methods are benchmarked and evaluated. The article concludes by
summarizing our impressions of the quantum model’s application and announcing
future steps.

##  2 Problem definition

MOFs consist of two key components: an inorganic metal cluster and an organic
molecule. The choice of metal and organic components influences the structure
and, consequently, the properties of the MOF. There is a broad range of
possibilities for combining metallic cations and organic ligands to construct
the MOF74 [5]. A well-designed MOF is capable of efficiently capturing
$\mathrm{CO_{2}}$ in DAC filters while hindering the clogging of the filter by
other gases. Screening for potentially suitable MOF74 structures involves
calculating the Dissociation Energy (DE) of a complex formed by MOF74 and gas
molecules. To calculate the DE of a MOF-gas complex, we require Potential
Energy Surface (PES) scans of both the complex and individual molecules, as
depicted in Equation 1.

$\mathrm{DE=PES_{gas-ion}-(PES_{gas}+PES_{ion})}$ (1)

The method for computing PES scans was developed using Variational Quantum
Eigensolver (VQE), which is defined as a hybrid quantum-classical, variation
principle-based algorithm[6]. The model is exemplified using the interaction
between a magnesium ion ($\mathrm{Mg^{2+}}$) and carbon dioxide gas
($\mathrm{CO_{2}}$).

Due to the structural complexity of the MOF-gas system, the latter must be
down-scaled for the analysis of quantum-based results. Fig. 1 presents an
example of system simplification, from a network of MOF unit cells, passing
through a singular unit cell with six open reaction sites (metallic ions), and
culminating in one preferential binding site (metal ion) - gas interaction,
neglecting the organic part [4].

Figure 1: The MOF structure simplification scheme

The electronic structure of $\mathrm{CO_{2}-Mg^{2+}}$, modeled using the
minimal basis set STO-3G, encompasses a total of 48 spin orbitals and
approximately 700,000 fermionic terms. During the construction of the
Electronic Hamiltonian, these spin orbitals require encoding using roughly 48
qubits, leading to the generation of over 12 million Pauli terms. To this day,
quantum computing resources are constrained by a limited number of qubits and
low circuit depth for both physical quantum computers and classical emulation
of quantum computing processes. Consequently, the number of spin orbitals was
reduced within the scope of this study using the AS Transformation, as
described in the Method section.

The challenge was to demonstrate the potential benefits of a quantum
computing-based method. To evaluate this, we benchmarked total energy
calculations using various computing approaches and methods, including:

  1. (a)

VQE algorithm as a hybrid quantum-classical computation implemented by perfect
qubits;

  2. (b)

VQE algorithm as a hybrid quantum-classical computation implemented by
superconducting physical qubits, i.e., using quantum hardware;

  3. (c)

Ab-initio classical methods: the low-cost Restricted Hartree-Fock (RHF) and
the higher-cost Coupled Cluster single and Double excitation (CCSD) [7];

  4. (d)

Complete Active Space (AS) Configuration Interaction (CASCI) analysis as the
classical reference method [8].

A 3D potential energy surface scan in quantum mechanics offers a comprehensive
exploration of molecular behavior by systematically altering atomic positions
and calculating corresponding potential energies [9, 10]. Such scans are
valuable for interpreting chemical reactions, reaction pathways, and
mechanisms, as well as pinpointing transition states and energy barriers. The
potential of VQE to produce PES scans that consider multiple degrees of
freedom has been also examined through the creation of
$\mathrm{CO_{2}-Mg^{2+}}$ 3D PES scan.

Finally, a draft approach to compute the energy of a full unit cell was
proposed. It is based on what is known as the classical hybrid approach,
involving the combination of low-cost computational methods on a larger scale;
and local high-cost computational methods focused on the open reaction sites
of a molecule [11]. Following a similar logic, we applied a low-cost classical
method on the unit cell - CO2 interaction and quantum computing-based
calculations on open reaction sites, which are metal ions in this case. The
obtained results are well-converged and they follow the expected trend.

##  3 Method

Geometry Optimization. The Avogadro software was used to define the optimized
geometry of $\mathrm{CO_{2}-Mg^{2+}}$. To evaluate different degrees of
freedom within the molecule, we selected the Z-matrix coordinate system. In
the optimization process using the Universal Force Field (UFF) method, we
executed the Steepest Descent method with 500 steps.

Active Space Transformation. The AS Transformation involves the construction
of an AS configuration by selecting a certain number of orbitals from the full
electronic structure. It is challenging to find the AS configuration that
contributes the most to the total energy since it relies on varying factors
such as the Highest Occupied Molecular Orbitals (HOMOs) and Lowest Unoccupied
Molecular Orbitals (LUMOs) within the electronic structure [8]. Quantum
computation is performed only for the defined AS. The number of HOMO and LUMO
orbitals chosen to construct the AS configuration was varied. By reducing the
number of contributor spin orbitals in the model, the computational resources
required for the calculations were also reduced (see Table 1).

Table 1: The overview of resource requirements for $\mathrm{CO_{2}-Mg^{2+}}$ at different levels of approximation | System
---
Size
| System
---
Abv.
| No. Spin Orbitals
---
= No. qubits
| No. Fermionic
---
Terms
| No. Pauli
---
Terms
| Hamiltonian’s
---
Circuit Depth
Full | Full | 48 | 730612 | 12207792 | 254329
5 HOMO - 5 LUMO | 5h5l | 20 | 20136 | 142940 | 7147
4 HOMO - 4 LUMO | 4h4l | 16 | 8292 | 46608 | 2913
3 HOMO - 3 LUMO | 3h3l | 12 | 2664 | 11076 | 923
2 HOMO - 2 LUMO | 2h2l | 8 | 564 | 1544 | 193
1 HOMO - 1 LUMO | 1h1l | 4 | 36 | 60 | 15

The following procedure was used to select the orbitals: Molecular Orbitals
(MOs) index was generated using the PySCF programming tool [12] and the
electronic Hamiltonian corresponding to the desired AS configuration was
constructed using the AS transformation module.

Energy Computation on Hybrid Quantum-Classical Framework. The PES computation
was performed using the IBM Qiskit Nature Module [13] with a PySCF interface.
We employed the VQE algorithm and designed it using the Qiskit library. The
VQE algorithm was executed on the statevector backend using perfect qubits but
also on quantum hardware using the Ibmq–nairobi backend via superconducting
qubits. Henceforth, we will refer to simulations based on perfect qubits as
"VQE - perfect q." and simulations based on physical qubits as "VQE - physical
q."

The singlet Hartree-Fock state was defined as the initial state for all
quantum computing-based computations. The variational form was constructed by
the Unitary Coupled Cluster Single and Double excitation (UCCSD) using fully-
entangled hardware-efficient ansatz. The ansatz implemented on the
ibmq–nairobi backend (VQE-physical q.) was parameterized by the rotational
parameters obtained from VQE-perfect q. computations. The fermionic terms of
selected AS configurations (1h1l, 2h2l, 3h3l, and 4h4l; see Table 1) were
mapped to a qubit Hamiltonian in the minimal basis set (STO-3G) using the
Bravy–Kitaev transformation. The two-qubit reduction technique was employed to
decrease the number of qubits by two in each generated qubit Hamiltonian
circuit. The SLSQP and SPSA optimizers were used to optimize the ansatz
parameters toward minimizing the electronic structure energy in statevector
and ibmq–nairobi simulators, respectively. To mitigate errors incurred during
qubit readout in physical q.-based implementation, we applied readout error
mitigation at resilience level one.

Energy Computation on Classical Framework. The PES of the full
$\mathrm{CO_{2}-Mg^{2+}}$ electronic structure was calculated in PySCF,
applying RHF and CCSD methods on the STO-3G basis set to obtain the full
electronic structure energy references. To obtain a relevant total energy
reference, with respect to the quantum computing-based method and AS
Transformation applied, the localized Full Configuration Interaction analysis
was performed. This was done by exact Diagonalization of different AS
configuration matrices embedded into the Hartree-Fock spin orbitals, which is
known as CASCI. The diagonalization was done using NumpyEigensolver.

##  4 Results and Discussion

The choice of an energy calculation method in quantum chemistry depends on the
balance between computational resources and the desired level of accuracy. The
system state is defined by the Schrödinger equation and several methods are
employed to solve it, varying in accuracy level. The mean-field approach, such
as HF, provides a simplified but computationally efficient representation of
electron-electron interactions, considering electrons as independent in an
average field generated by all other electrons. Since it neglects the
electrons correlation effect, it leads to a low accuracy calculation,
especially for strongly correlated systems [14]. CCSD is an extension or
correction with respect to the HF. It accounts for some electron correlation,
regarding various levels of excitation and often yields highly accurate
results for many molecular properties. The computational cost for CCSD is
expensive with respect to the other classical methods [15]. The most accurate
manner to simulate the molecular systems is by considering all possible
electronic configurations and interactions, this could be done by the FCI
approach. FCI becomes computationally intractable for large systems due to its
exponential scaling. Local application of FCI on a subset of important
configurations (CASCI), strikes a balance between accuracy and computational
cost [8]. The configuration interaction size selected is dictated by the
computational power of the machine used and even the most powerful classical
machines fail to perform FCI on large molecular systems.

VQE belongs to the category of variational methods, unlike FCI which is a
deterministic method. VQE aims to compute an upper bound for the lowest
possible expectation value of an observable with respect to a trial wave
function. The objective of the VQE is therefore to find a parametrization of
the trial wave function, such that the expectation value of the Hamiltonian
(energy equation) is minimized. VQE is assumed to achieve higher accuracy than
classical mean-field methods but at a lower cost compared to FCI. However, it
is worth mentioning that its accuracy depends on factors like the choice of
initial state, ansatz, optimization method, execution approach, etc. [16, 17].

The Total Energies at the Optimized Geometry. To compare the accuracy level of
energy computation for the optimal inter-nuclear distance, the total energies
of the electronic structure of $\mathrm{CO_{2}-Mg^{2+}}$ (ground state energy)
are shown as a function of AS configuration size (number of HOMOs and LUMOs
orbitals involved), see Fig. 2. Fig. 2 (A) schematically presents the
different AS configurations. The HOMO and LUMO orbitals chosen to create an AS
configuration are marked blue. The remaining inactive orbitals are marked red.
Different markings highlight that CASCI and VQE computation were applied on AS
configurations (blue orbitals), while the mean-field (HF) was applied on
inactive space (red). To obtain the total energy, CASCI/ VQE was embedded into
the mean field calculation.

Figure 2: A) The schematic presentation of AS configurations B) The
$\mathrm{CO_{2}-Mg^{2+}}$ relative energy with respect to the RHF($\Delta E$)
and absolute energy (Energy) at optimal bond length. For the methods where AS
transformation was applied, $\Delta E$ and Energy are presented as a function
of AS Configuration interaction size.

In Fig. 2 (B), RHF and CCSD results shown as dotted lines, illustrate
classical method thresholds of low and higher accuracy, respectively. One
should note that these two methods were applied to the entire
$\mathrm{CO_{2}-Mg^{2+}}$ electronic structure. Two different energy scales
were used to evaluate the data: total energy (right axes) and total energy
improvement with respect to the RHF accuracy level (left axes, $\Delta$E =
$\mathrm{E-E_{RHF}}$). The lowest computed energy, alongside its corresponding
methods and available computing resources, is used to represent the highest
accuracy level (solid line).

It can be observed that accuracy significantly increases when more than 1h1l
orbitals are involved. AS configuration size is expected to influence energy
value [18]. VQE-Perfect q. 2h2l result is comparable to CASCI 2h2l. This
demonstrates that the hybrid model performs comparably well as its classical
reference, and both are above the CCSD level. The deviation observed for VQE-
physical q. 2h2l implementation in comparison to the perfect qubit and
classical outcome is a result of erroneous hardware execution. With increasing
AS space size, considering the investigated range, the quantum computing-based
method shows no significant improvement in the accuracy level. Contrarily,
energy computed with CASCI decreases with increasing AS size. The lowest
energy was achieved for CASCI implementation with 4h4l.

Potential Energy Surface Scan. Fig. 3 (A-D) presents total energy evolution as
a function of inter-nuclear distance i.e. PES scan, obtained by various
methods (see section Method) for 1h1l, 2h2l, 3h3l, and 4h4l configurations.
RHF and CCSD are also shown in the figure representing full electronic system
classical computations.

Figure 3:  (A) (B), (C), and (D) are the $\mathrm{CO_{2}-Mg^{2+}}$ PES scans
for 1h1l, 2h2l, 3h3l, and 4h4l AS configuration as a function of internuclear
distance, respectively, including the RHF and CCSD results for the full
electronic structure. (E), (F), (G), and (H) represent the deviation VQE
results with 1h1l, 2h2l, 3h3l, and 4h4l configurations from their own
corresponding Localized CASCI energy reference, respectively. RHF and CCSD
deviations are related to the calculation of the full electronic structure
including the RHF and CCSD results.

For 1h1l computation, VQE and CASCI computations meet RHF accuracy level. This
can be expected, due to a low number of orbitals considered within 1h1l AS
configuration. In general, one can notice that as the AS size expands,
accuracy gradually increases. However, this trend does not hold true for the
VQE calculations, with both perfect q. and physical q. Nonetheless, the data
obtained by the VQE algorithm are well-converged and demonstrate a similar
trend to the energy evolution calculated on classical computers. 2h2l, 3h3l,
and 4h4l PES scan trends show a difference with respect to the full
configuration PES scans (RHF and CCSD) at about 2.6 Å. Namely, AS results
exhibit energy decay for distances larger than $\sim$ 2.6 Å.

The observed deviation between VQE and CASCI results is analyzed in Fig 3
(E-H). For a particular AS-Configuration, CASCI computation is taken as a
reference, and other computations are illustrated relative to their
corresponding Localized CASCI. For the 1h1l configuration, it can be observed
again that all the results overlapped with CASCI, except the CCSD, which
reaches higher accuracy. VQE-perfect q. 2h2l exhibits a reasonable overlap
with CASCI 2h2l, both being more accurate with respect to the RHF and CCSD.
VQE-physical q. deviation resulting from erroneous behavior of actual quantum
hardware can be observed, as expected. However, 3h3l and 4h4l configurations
suffer from discrepancies between VQE and CASCI data at low distances. It is
expected that these deviations become more obvious with the AS size
increasing. Still, it is interesting to notice that energies obtained by VQE,
for both perfect q. and physical q, seem to follow their corresponding CASCI
results. Moreover, for the STO-3G minimal basis set, these surpass the
accuracy of the CCSD, despite CCSD holding a reputation of a computationally
demanding approach.

The potential of quantum computing-based methods can be envisioned through the
fair agreement between the observed VQE-perfect q. and CASCI results. To
further improve energy computation, one should enlarge an AS size. However,
this is computationally expensive, and classical resources would be consumed
before reaching a significant number of orbitals, especially for large and
strongly correlated systems. With the assumption that quantum computing
resources will grow, it can be presumed that quantum computing-based methods
play a key role in reaching better accuracy levels for chemical system energy
calculations. It’s important to emphasize that the performance of quantum
hardware would need to improve concurrently, as indicated by the deviation in
the results between VQE-physical q. and VQE-perfect q. It was observed that
the discrepancy between VQE and CASCI increases with increasing the AS
configuration size, which suggests that VQE performance decreases with
enlarging an electronic structure matrix. This demonstrates that further
algorithm improvements and development are needed in order to reach an
industry-relevant quantum computing application.

3D Potential Energy Surface Scan. Different adsorption orientation angles for
the $\mathrm{CO_{2}-Mg^{2+}}$ interaction were incorporated with existing
calculations over the bond distance. This resulted in a 3D PES scan (Fig. 4),
which is obtained by VQE-Perfect q.

(a)

(b)

Figure 4: (a) 3D PES scan of the $\mathrm{CO_{2}-Mg^{2+}}$ obtained by
variation of two degrees of freedom and, (b) Schematic of different
orientations of $\mathrm{CO_{2}}$ on a hypothetical $\mathrm{CO_{2}-Mg^{2+}}$
embedded in a scaled-up MOF74 structure

Results presented in Fig 4. show that VQE can be used successfully to
manipulate two degrees of freedom. Still, one should note that the calculation
is not complete, with respect to the range of analysis needed for forming
chemistry-relevant conclusions about, for example MOF - CO2 reaction path.

The MOF74 Unit Cell Energy Computation. The idea behind the proposed
computation is to draft a direction toward the MOF74 full unit-cell energy
calculation. A local energy correction (LEC) method was applied as
schematically presented in Fig 5. The energy of the MOF74 unit cell, without
$\mathrm{CO_{2}}$ gases involved $(\mathrm{E(RHF_{S})})$ was calculated with
RHF, representing the low-cost method. The energy of the $\mathrm{CO_{2}}$
interaction with the most active molecule part, which is six open reaction
sites (6 $\mathrm{Mg^{2+}}$ ), was calculated with VQE-perfect q.
$(\mathrm{E(VQE_{S})})$. In addition, RHF was applied to the open reaction
sides $(\mathrm{E(RHF)_{C})})$, as it is required for final energy
determination. The final $\mathrm{E_{LEC}}$ was computed as presented in Eq.
2. Existence of preferential binding sites and entanglement between chemical
structures calculated by RHF and VQE were not considered within this study.

$\mathrm{E_{LEC}=E(RHF_{S})+E(VQE_{C})-E(RHF_{C})}$ (2)

Figure 5: Schematic of LEC method application for the PES computation of 6Mg-
MOF74 - $\mathrm{CO_{2}}$ complex Figure 6: (Left Axis:) The
$\mathrm{CO_{2}-6Mg-MOF74}$ relative total energies with respect to the RHF
($\Delta$E) and (Right Axis:) absolute total energies at optimal bond length.
For the methods where CASCI and VQE approaches were applied, both are
presented as a function of growing AS configurations.

Fig. 6 presents the total energy computation for CO2 - 6Mg MOF74 interaction
and relative energy with respect to the RHF reference value, $\Delta$E. In
addition, VQE-perfect q. 4h4l result is presented as a solid line to
illustrate thresholds in the same way as in Fig 2. It can be observed that the
RHF accuracy level is improved by applying quantum computing-based local
energy correction.

##  5 Conclusion and future work

This work aimed to evaluate if MOF design could be enhanced by quantum
computing in the future. To investigate that, a quantum computing-based energy
calculation model was constructed following the PISQ approach and a comparison
of perfect and physical q.-based implementation was done in parallel with ab-
initio (RHF, CCSD) and CASCI computations. The model was based on complete
active space configuration interaction. The potential of quantum has been
demonstrated through well agreement of VQE-perfect q. results with localized
FCI-based calculations (CASCI). However, an attempt to scale up the model to
involve more molecular orbitals and hardware implementation revealed drawbacks
of current quantum computing: low number of qubits available, erroneous
behavior of qubits, and algorithm design challenges. A pathway toward MOF 74
full unit cell energy computation was proposed. It is based on low and high-
cost model applications on different MOF unit cell building structures. Well-
converged results followed the trend as expected. However, the energy
computation needs to involve existing structure interactions in a more
representative manner. Perfect qubit-based simulation proved to be a valuable
approach to perform model evaluation. In addition, VQE-perfect q. results were
used to parameterize VQE-physical q. implementation, which shortens hardware
runtimes. Perfect qubit simulation was used to draft a pathway toward modeling
3D and full unit cell PES scans. These procedures do not reach chemistry-
relevant outcomes yet, and implementing them on hardware brings no advantage.
However, it is needed to validate the proof of concept and the PISQ approach
is suitable. These indicate some of the most important benefits of the PISQ
approach. The authors plan to expand the work by performing hardware
implementation beyond the 2h2l AS configuration, to enable further analysis of
the quantum potential. More investigation will be done to advance unit cell
energy computation.

## References

  * [1] John Preskill  “Quantum Computing in the NISQ era and beyond”  In _Quantum_ 2 Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften, 2018, pp. 79  DOI: 10.22331/q-2018-08-06-79
  * [2] Koen Bertels, Tamara Sarac, Aritra Sarkar and Imran Ashraf  “Quantum Computing – from NISQ to PISQ”, 2022  arXiv:2106.11840 [quant-ph]
  * [3] Richard P Feynman  “Simulating physics with computers”  In _Int. j. Theor. phys_ 21.6/7, 2018
  * [4] Allison L. Dzubak et al.  “Ab initio carbon capture in open-site metal–organic frameworks”  In _Nature Chemistry_ 4.10, 2012, pp. 810–816  DOI: 10.1038/nchem.1432
  * [5] Abhishek Mitra et al.  “Periodic Density Matrix Embedding for CO Adsorption on the MgO(001) Surface”  In _The Journal of Physical Chemistry Letters_ 13.32 American Chemical Society, 2022, pp. 7483–7489  DOI: 10.1021/acs.jpclett.2c01915
  * [6] Alberto Peruzzo et al.  “A variational eigenvalue solver on a photonic quantum processor”  In _Nature communications_ 5.1 Nature Publishing Group, 2014, pp. 1–7
  * [7] Rodney J. Bartlett and Monika Musiał  “Coupled-cluster theory in quantum chemistry”  In _Rev. Mod. Phys._ 79 American Physical Society, 2007, pp. 291–352  DOI: 10.1103/RevModPhys.79.291
  * [8] Björn O. Roos, Peter R. Taylor and Per E.M. Sigbahn  “A complete active space SCF method (CASSCF) using a density matrix formulated super-CI approach”  In _Chemical Physics_ 48.2, 1980, pp. 157–173  DOI: https://doi.org/10.1016/0301-0104(80)80045-0
  * [9] Dequan Wang et al.  “Accurate Potential Energy Surfaces for the Three Lowest Electronic States of N((2)D) + H(2)(X(1)$\sum$(g) (+)) Scattering Reaction”  In _ACS Omega_ 4.7, 2019, pp. 12167–12174
  * [10] Zijie Yan, Stephen K Gray and Norbert F Scherer  “Potential energy surfaces and reaction pathways for light-mediated self-organization of metal nanoparticle clusters”  In _Nature Communications_ 5.1, 2014, pp. 3751
  * [11] Loredana Valenzano, Bartolomeo Civalleri, Kaido Sillar and Joachim Sauer  “Heats of Adsorption of CO and CO2 in Metal–Organic Frameworks: Quantum Mechanical Study of CPO-27-M (M = Mg, Ni, Zn)”  In _The Journal of Physical Chemistry C_ 115.44 American Chemical Society, 2011, pp. 21777–21784  DOI: 10.1021/jp205869k
  * [12] Qiming Sun et al.  “PySCF: the Python-based simulations of chemistry framework”  In _WIREs Computational Molecular Science_ 8.1, 2018, pp. e1340  DOI: https://doi.org/10.1002/wcms.1340
  * [13] Qiskit contributors  “Qiskit: An Open-source Framework for Quantum Computing”, 2023  DOI: 10.5281/zenodo.2573505
  * [14] Attila Szabo and Neil S Ostlund  “Modern quantum chemistry: introduction to advanced electronic structure theory”  Courier Corporation, 2012
  * [15] Isaiah Shavitt and Rodney J Bartlett  “Many-body methods in chemistry and physics: MBPT and coupled-cluster theory”  Cambridge university press, 2009
  * [16] Sam McArdle et al.  “Quantum computational chemistry”  In _Reviews of Modern Physics_ 92.1 APS, 2020, pp. 015003
  * [17] Jules Tilly et al.  “The Variational Quantum Eigensolver: A review of methods and best practices” The Variational Quantum Eigensolver: a review of methods and best practices  In _Physics Reports_ 986, 2022, pp. 1–128  DOI: https://doi.org/10.1016/j.physrep.2022.08.003
  * [18] Valera Veryazov, Per Åke Malmqvist and Björn O. Roos  “How to select active space for multiconfigurational quantum chemistry?”  In _International Journal of Quantum Chemistry_ 111.13, 2011, pp. 3329–3338  DOI: https://doi.org/10.1002/qua.23068