# Deep Learning for Spatiotemporal Big Data: A Vision on Opportunities and
Challenges

Zhe Jiang  Department of Computer & Information Science & Engineering
University of Florida, Gainesville, FL, USA
<EMAIL_ADDRESS>
###### Abstract

With advancements in GPS, remote sensing, and computational simulation, an
enormous volume of spatiotemporal data is being collected at an increasing
speed from various application domains, spanning Earth sciences, agriculture,
smart cities, and public safety. Such emerging geospatial and spatiotemporal
big data, coupled with recent advances in deep learning technologies, foster
new opportunities to solve problems that have not been possible before. For
instance, remote sensing researchers can potentially train a foundation model
using Earth imagery big data for numerous land cover and land use modeling
tasks. Coastal modelers can train AI surrogates to speed up numerical
simulations. However, the distinctive characteristics of spatiotemporal big
data pose new challenges for deep learning technologies. This vision paper
introduces various types of spatiotemporal big data, discusses new research
opportunities in the realm of deep learning applied to spatiotemporal big
data, lists the unique challenges, and identifies several future research
needs.

###### Index Terms:

spatiotemporal big data, deep learning

##  I Introduction

With advancements in GPS, remote sensing, and computational simulation, a
rapidly expanding volume of spatiotemporal data is being collected from
diverse application domains, ranging from Earth sciences to agriculture, smart
cities to public safety [1]. For instance, NASA’s Earth Observing System Data
and Information System (EOSDIS) alone houses 90 petabytes of Earth science
data that continues to grow daily [2]. These emerging spatiotemporal data,
characterized by their volume, velocity, and variety, which surpass the
capabilities of traditional computational platforms, are referred to as
_spatiotemporal big data_ (STBD) [3].

Analyzing geospatial and spatiotemporal data has attracted researchers’
interest since early 2000, fostering the development of spatial and
spatiotemporal data mining [4, 5, 6, 7], a research field that aims to
discover non-trivial, previously unknown, but potentially useful patterns or
make accurate predictions on massive geospatial and spatiotemporal data. In
this context, the word “massive” refers to a data size that exceeds the
capability of traditional statistical methodologies (geostatistics, time
series analysis) [8], necessitating efficient algorithms such as computational
pruning. The rapid growth of this interdisciplinary field not only enhances
the awareness of available tools for domain scientists to make sense of their
spatiotemporal data but also inspires data mining researchers to develop novel
algorithms that address unique characteristics of spatiotemporal data, where
samples violate the common independent and identical distribution (i.i.d.)
assumption.

The emergence of STBD in the last decade has inspired the development of
spatial big data computational platforms (e.g., Hadoop-GIS [9], SpatialHadoop
[10], GeoSpark [11]). These platforms make it possible to analyze patterns in
STBD that are beyond the capability of traditional spatiotemporal data mining
algorithms. However, most of these systems focus on computational scalability
in simple query processing or visual pattern analysis, instead of complex
predictive learning tasks.

The recent advances in deep learning technologies have brought unique
opportunities. Thanks to the availability of big training data (e.g.,
ImageNet), efficient computational architectures (e.g., GPUs), and smarter
models and learning algorithms, deep learning has achieved tremendous success
in computer vision and language processing applications [12]. It has already
made promising progress in spatial and spatiotemporal domains such as Earth
system process modeling [13], drug discovery [14], transportation and urban
planning [15].

When STBD meets deep learning, new opportunities arise. First, deep neural
networks can provide end-to-end training without handcrafting spatial and
temporal contextual features [13]. Traditionally, in spatiotemporal data
mining, a major step is feature engineering based on domain knowledge. This
approach has limitations. For instance, researchers have to repeat this
tedious and time-consuming process for every new application, and the
handcrafted features may be insufficient to capture complex spatial and
temporal patterns. Second, some complex spatiotemporal patterns are only
visible in STBD and thus are not detectable by traditional spatial and
spatiotemporal data mining methods. In contrast, deep learning models are very
data-hungry with a large number of model parameters (e.g., millions to
billions). The huge capacity of deep learning models provides opportunities to
learn hidden patterns from STBD that are otherwise ignored by traditional
models.

This vision paper will briefly discuss the emerging opportunities of deep
learning for spatiotemporal big data, highlight the unique challenges, and
identify new research directions for this highly interdisciplinary area.

##  II Types of STBD

###  II-A Spatial data

_Point reference data_ : Spatial point reference data record non-spatial
attribute values at specific point locations in continuous space. The
assumption is that the attribute exists everywhere in the spatial domain but
is only observed at these sampled point locations. For example, this could
include the precipitation measurements obtained from various weather stations,
each representing a point location.

_Areal data_ : Areal data capture non-spatial attribute values across the
entire continuous space by dividing it into cells, either regular or
irregular, and assigning an aggregated value to each cell. An example of areal
data is a map displaying the number of COVID-19 cases in all U.S. counties.

_Spatial point process_ : Spatial point processes involve the recording of
spatial point locations where specific types of events occur, such as crimes
or traffic accidents. Unlike point reference data, the locations of events in
a spatial point process are treated as random variables rather than non-
spatial attributes. Examples include the spatial distribution of traffic crash
locations or crime events.

_Spatial network_ : Spatial networks are graphs in which nodes and edges
correspond to spatial objects. A common example is a road network, where nodes
represent road intersections, and edges represent road segments.

###  II-B Spatiotemporal data

When considering the additional time dimension on top of spatial data, we can
divide the time dimension into two cases: regular time intervals (discrete
time) or irregular time intervals (continuous time). Based on this
categorization, we can further identify the following common types of
spatiotemporal data.

_Spatiotemporal point reference data_ : Spatiotemporal point reference data
expand upon point reference data by including non-spatiotemporal attributes
sampled in both continuous space and continuous time. An example is water
quality samples collected by boats at various locations and timestamps.

_Spatial point time series_ : Spatial point time series extend point reference
data by incorporating regular time intervals for sampling non-spatial
attributes. For instance, hourly temperature measurements from a network of
weather stations, with spatial locations sampled in continuous space, fall
into this category.

_Areal temporal snapshots (or time series)_ : Areal temporal snapshots or time
series add a time dimension to areal models, with observations at irregular or
regular time intervals. For example, ecologists might collect Earth imagery
during different seasons (e.g., leaf-on and leaf-off) to monitor changes in
wetland locations.

_Spatiotemporal point process_ : Spatiotemporal point processes extend spatial
point processes by characterizing both the location and time of point events.
An example is traffic crash records, which include information about both the
location and time of each incident.

_Trajectories_ : Trajectories represent a specialized case of spatiotemporal
point processes where points are regularly sampled in time, and the locations
follow a sequential order. This is often seen in tracking systems for vehicles
or individuals.

_Spatial network time series_ : Spatial network time series expand upon
spatial network data by including time series data for non-spatiotemporal
attributes associated with network nodes or edges. For instance, sensors
monitoring traffic speed and flow volume at different road segments provide
insights into transportation dynamics over time.

##  III New Opportunities of Deep Learning for Spatiotemporal Big Data

###  III-A Why is deep learning for STBD different?

There are several aspects that make deep learning transformative compared with
traditional spatial and spatiotemporal data mining techniques for STBD.

First, deep neural networks enable end-to-end training, eliminating the need
for manually crafting spatial contextual features. In traditional spatial and
spatiotemporal machine learning, a significant step involves handcrafting
these features based on domain knowledge. Deep learning mitigates this by
learning relevant features from data directly. This approach offers advantages
such as reducing the tedium and time consumption associated with feature
engineering for each new application. For instance, in Earth science
applications like land cover mapping, deep convolutional neural networks can
learn rich spatial and spectral signatures directly from remote sensing image
pixels.

Second, deep learning models, characterized by their extensive number of
parameters (often ranging from millions to billions), are data-hungry and
possess a significantly larger capacity to uncover intricate spatial and
temporal patterns. This is particularly crucial for understanding complex
spatiotemporal processes involving multiple interacting physical variables
operating at various spatial and temporal scales. Deep neural networks are
well-suited to capture such complex patterns from STBD.

Third, the recent success of large language models [16], also known as
foundation models, in language processing tasks offers insights into
addressing a persistent challenge in spatiotemporal data mining: model
generalizability across heterogeneous space and time. Traditionally, due to
spatial heterogeneity, models trained in one region may perform poorly in
another, and the same holds for different application tasks [6]. Large
language models have demonstrated that with sufficient training data, it’s
possible to train a general model capable of performing well across diverse
geospatial regions and various downstream geoscience tasks. For example, in
remote sensing, where different applications may share common spectral
signatures in satellite imagery, it becomes feasible to train a single neural
network model to extract shared feature representations from the vast Earth
imagery big data.

###  III-B Illustrative Transformative Applications

_Foundation model for Earth imagery big data_ : A foundation model is a large
deep neural network model trained on a vast quantity of data at scale (often
by self-supervised or semi-supervised learning) such that it can be easily
adapted to a wide range of downstream tasks [17]. For example, a visual
foundation model such as DINOv2 can produce all-purpose features from an image
that work across image distributions and tasks without the need for finetuning
[18]. Given the petabytes of Earth imagery collected by NASA and other
commercial satellite companies, there is growing interest in developing such
foundation models in Earth sciences [19]. Such an Earth imagery foundation
model is transformative to geoscience research. For many years, a geoscience
researcher has to first handcraft spectral signatures for specific tasks. Even
with the popularity of deep convolutional neural networks, people still need
to collect a large number of training labels for different application [20].
With a foundation model, a user can potentially do any downstream task in any
region with only a small number of labels for fine-tuning. However, several
caveats exist. Thematic classes in geoscience are often not limited to
discrete objects (e.g., buildings, roads, cars). Techniques for learning
general features for object and image classes in computer vision may not be
effective for land cover and land use classes. Moreover, scientists often
prefer a model they can understand, interpret, and trust.

_Spatiotemporal AI surrogate for Earth system models:_ Another potential
application is the development of surrogate models for numerical simulation
models in Earth sciences (e.g., weather forecasting, hydrology, oceanography)
[21, 22]. For example, coastal circulation refers to the movement and
circulation patterns of water in the coastal regions of oceans and seas [23,
24]. Effective and efficient coastal circulation modeling plays a crucial role
in building early warning and forecasting systems for coastal hazards (e.g.,
floods and storm surges, marine heatwaves, red tides) [25, 26]. Such models
often contain hundreds of thousands of lines of code that run hundreds of CPUs
in a high-performance computing platform. For instance, the Regional Ocean
Model System (ROMS) is a free-surface, terrain-following, primitive equations
ocean model widely used by the scientific community to create high-resolution
simulations of ocean and coastal processes, including circulation,
temperature, salinity, and other variables [27]. Very often, scientists need
to run the simulation models multiple times for parameter calibration and
uncertainty quantification. With the popularity of deep learning and GPU
computational architectures, scientists have become increasingly interested in
training machine learning surrogates for numerical models based on a vast
number of simulation data [28]. The advantage is that once a model is well-
trained, it can run on a desktop or even a laptop with one or a few GPUs. This
is very important for researchers and stakeholders in developing countries
without access to significant computational resources.

##  IV Unique Challenges

Although STBD shares similar formats with image and video data (spatial raster
resembles imagery, and temporal snapshots resemble video frames), it poses
some unique challenges for deep learning technologies, as summarized below.

First, as the space is continuous, the input size of a spatial raster map can
be easily many thousand by many thousand pixels. This is much larger than the
common input shapes of existing deep learning models, which are a few hundred
by a few hundred pixels. Moreover, there exist complex long-range spatial
dependencies across multiple spatial scales. Thus, modeling the complex
spatial dependency within a large input map poses computational challenges
such as a large GPU memory footprint.

Second, spatial data contain multiple map layers co-registered in the same
spatial reference system. Each individual map layer can have a distinct
spatial resolution and format (vector versus raster) and diverse noise and
errors. For example, in water forecasting, the input feature layers can
include satellite images for land cover maps and climate model projections,
sensor observations at weather stations, river network polylines, and lake
polygons.

Third, training deep neural networks requires a large amount of ground truth
labels (e.g., ImageNet). Unfortunately, in spatiotemporal applications,
collecting ground truth is often slow, tedious, and expensive due to the need
for field surveys or visually interpreting high-resolution Earth imagery.

Fourth, spatiotemporal problems are often highly interdisciplinary and require
the integration of deep learning with physical knowledge [29]. This is because
machine learning problems are often under-constrained. Without incorporating
domain physics, pure data-driven deep learning can easily learn spurious
relationships that look deceptively good on the training data (even on
validation and test data) but do not generalize well to complex real-world
scenarios. Indeed, how to integrate universal physics into machine learning to
enhance model generalizability for a new region with limited labels has been
listed as an open research question by the CRA’s report on _20-Year Roadmap
for Artificial Intelligence_ [30].

Finally, compared with many physics-based models in scientific domains (e.g.,
hydrology, oceanography), deep neural network models are black-box in nature
with limited explainability, which is a concern as the end goal is to advance
knowledge and understanding of a subject. Moreover, as it is difficult to
explain why a deep neural network works, it is also hard to determine when it
will fail. The issue becomes more serious when it comes to model
generalizability in a complex real-world environment. Thus, uncertainty
quantification is of great importance [31].

##  V Future Research Directions

Addressing the challenges listed above requires additional research in the
design of novel deep-learning models and learning algorithms for STBD. Several
future research directions are listed here.

###  V-A Learning complex spatial dependency in continuous space

In traditional spatiotemporal data mining, the common practice of prediction
is to first establish a spatial neighborhood graph and then build a model to
capture the dependency structure within the neighborhood graph [32]. For deep
learning models, neighborhood dependency can be learned by convolutional
filters. However, for long-range dependency beyond the neighborhood,
convolutional operations are no longer sufficient. One potential direction is
to leverage the more recent attention mechanisms that are popular in vision
transformers [33]. However, due to the high computational costs, we cannot
compute all pair attention across all pixels (or even all pixel blocks). This
problem is even worse for a large spatial raster with many thousands by many
thousands of pixels. One potential direction is to capture such long-range
spatial dependency at different scales with some model components to learn
spatial dependency within smaller patches and other model components to learn
the dependency across patches. Another non-trivial issue exists for a large
number of spatial point samples in continuous space. This requires new
positional encoding and efficient design of self-attention layers across a
large number of point samples.

###  V-B Weakly supervised spatiotemporal deep learning

Although there is a paucity of high-quality ground truth in spatiotemporal
applications, there are often abundant weak labels that are sparse, coarse,
and noisy. For example, in flood mapping applications, the ground truth from
water sensors can be spatially sparse, and human-annotated flood extent
boundary polygons on high-resolution imagery are also imprecise (e.g., due to
obscurity from tree canopies). Similarly, in land cover mapping at a large
scale (e.g., continental or global), the weak labels can come from existing
classification products based on low-resolution satellite imagery. Addressing
this issue requires novel weak supervised learning [34] algorithms that can
train effective spatiotemporal deep learning models from sparse, coarse, and
noisy labels. For sparse labels, semi-supervised learning or domain-knowledge-
guided learning can potentially help fill the gap of incomplete labels. For
coarse training labels, the learning algorithms may require inferring the
hidden true label in a higher resolution iteratively while training an end
model. The process may involve multi-instance learning. For noisy labels, a
further noise model needs to be assumed to capture the probabilistic
relationship between observed noisy labels and the hidden true labels to
facilitate label inference during model training. The process becomes more
complex if the input labels are in the vector format (e.g., points, polylines,
and polygons) while the input features are in spatial raster imagery [35, 36].

###  V-C Knowledge-guided spatiotemporal deep learning

There are different sources of domain physical knowledge. One source is
symbolic knowledge in a geographical knowledge graph or knowledge base (with
spatial logic rules). To integrate such knowledge into spatiotemporal deep
learning requires neurosymbolic models. One potential challenge is that the
underlying phenomenon is in continuous space without explicit object
definitions. How to ground a geospatial knowledge base (logic rules) into
discrete locations involves a trade-off between computational efficiency and
spatial granularity [37]. Another source of geographical knowledge is
numerical physical models based on partial differential equations (PDEs).
Examples include weather forecasting models, ocean circulation models, and
hydrologic models. There are already extensive studies in using deep neural
networks to solve PDEs, such as physics-informed neural networks (PINNs) [38]
and neural operator learning [39]. PINN aims to train a neural network as a
function of spatial and temporal coordinates. Samples on the initial and
boundary conditions are considered as ground truth and other randomly sampled
locations in between are considered as regularization terms based on the PDE
constraint. Neural operator learning aims to train a neural network as an
operator for a family of PDEs, taking in samples on the initial and boundary
conditions as network inputs and producing a solution function as the neural
network outputs. The advantage of PINNs is that it does not require solving a
PDE first, as the neural network training process is solving the PDE, possibly
at a faster speed due to the efficient GPU parallelization and automatic
differentiation. However, the neural network itself is still a simple function
of coordinates, often multilayer perceptrons, without explicitly encoding the
spatial dependency between different sample locations. Training a PINN for
complex real-world PDEs can be challenging. The network needs to be retrained
for a new PDE instance. As for neural operator learning, it only requires
model training once for all but requires a significant amount of simulations
from numerical solvers. It is often challenging to collect a large number of
training instances for real-world scenarios due to the high time costs of
solving those PDEs.

###  V-D Model generalizability across heterogeneity

There are several issues to be addressed here. First, to address the challenge
of spatial and temporal heterogeneity, one potential idea is to design new
learning algorithms to be able to train an ensemble of models (each base model
is trained for one homogeneous sub-area). For example, meta-transfer learning
could be used to effectively train a shared module across different domain
tasks [40]. For spatial problems, the definition of different domains (e.g.,
homogeneous sub-regions) can be implicit. Several works have been done to
address this issue [41, 42]. Another more challenging issue is the out-of-
distribution (OOD) scenario. A well-trained model with high training and
validation accuracy can still show poor performance when being deployed in a
complex real-world environment. To address this issue, the model needs to be
able to identify when a test case is an OOD case and accurately quantify the
uncertainty of its prediction, e.g., based on sample distance or density in
feature space or geographical space [35].

## Acknowledgment

This material is based upon work supported by the National Science Foundation
(NSF) under Grant No. IIS-2147908, IIS-2207072, CNS-1951974, OAC-2152085.

## References

  * [1] E. Eftelioglu, Z. Jiang, X. Tang, and S. Shekhar, “The nexus of food, energy, and water resources: Visions and challenges in spatial computing,” in _Advances in Geocomputation: Geocomputation 2015–The 13th International Conference_. Springer, 2017, pp. 5–20.
  * [2] NASA, “NASA’s Earth Science Data Systems (ESDS) Program,” https://www.earthdata.nasa.gov/technology/open-science, Last Accessed on August 20, 2023.
  * [3] Z. Jiang and S. Shekhar, _Spatial Big Data Science - Classification Techniques for Earth Observation Imagery_. Springer, 2017.
  * [4] S. Shekhar, Z. Jiang, R. Y. Ali, E. Eftelioglu, X. Tang, V. Gunturi, and X. Zhou, “Spatiotemporal data mining: A computational perspective,” _ISPRS International Journal of Geo-Information_ , vol. 4, no. 4, pp. 2306–2338, 2015.
  * [5] G. Atluri, A. Karpatne, and V. Kumar, “Spatio-temporal data mining: A survey of problems and methods,” _ACM Computing Surveys (CSUR)_ , vol. 51, no. 4, pp. 1–41, 2018.
  * [6] Z. Jiang, “A survey on spatial prediction methods,” _IEEE Transactions on Knowledge & Data Engineering_, vol. 31, no. 09, pp. 1645–1664, 2019.
  * [7] A. Sharma, Z. Jiang, and S. Shekhar, “Spatiotemporal data mining,” _Handbook of Spatial Analysis in the Social Sciences_ , p. 352, 2022.
  * [8] N. Cressie and C. K. Wikle, _Statistics for spatio-temporal data_. John Wiley & Sons, 2015.
  * [9] A. Aji, F. Wang, H. Vo, R. Lee, Q. Liu, X. Zhang, and J. Saltz, “Hadoop-gis: A high performance spatial data warehousing system over mapreduce,” in _Proceedings of the VLDB endowment international conference on very large data bases_ , vol. 6, no. 11. NIH Public Access, 2013.
  * [10] A. Eldawy and M. F. Mokbel, “Spatialhadoop: A mapreduce framework for spatial data,” in _2015 IEEE 31st international conference on Data Engineering_. IEEE, 2015, pp. 1352–1363.
  * [11] J. Yu, J. Wu, and M. Sarwat, “Geospark: A cluster computing framework for processing large-scale spatial data,” in _Proceedings of the 23rd SIGSPATIAL international conference on advances in geographic information systems_ , 2015, pp. 1–4.
  * [12] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” _nature_ , vol. 521, no. 7553, pp. 436–444, 2015.
  * [13] M. Reichstein, G. Camps-Valls, B. Stevens, M. Jung, J. Denzler, N. Carvalhais, and f. Prabhat, “Deep learning and process understanding for data-driven earth system science,” _Nature_ , vol. 566, no. 7743, pp. 195–204, 2019\.
  * [14] J. Jumper, R. Evans, A. Pritzel, T. Green, M. Figurnov, O. Ronneberger, K. Tunyasuvunakool, R. Bates, A. Žídek, A. Potapenko _et al._ , “Highly accurate protein structure prediction with alphafold,” _Nature_ , vol. 596, no. 7873, pp. 583–589, 2021.
  * [15] X. Yin, G. Wu, J. Wei, Y. Shen, H. Qi, and B. Yin, “Deep learning on traffic prediction: Methods, analysis, and future directions,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 23, no. 6, pp. 4927–4943, 2021.
  * [16] T. Brown, B. Mann, N. Ryder, M. Subbiah, J. D. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell _et al._ , “Language models are few-shot learners,” _Advances in neural information processing systems_ , vol. 33, pp. 1877–1901, 2020.
  * [17] R. Bommasani, D. A. Hudson, E. Adeli, R. Altman, S. Arora, S. von Arx, M. S. Bernstein, J. Bohg, A. Bosselut, E. Brunskill _et al._ , “On the opportunities and risks of foundation models,” _arXiv e-prints_ , pp. arXiv–2108, 2021.
  * [18] M. Oquab, T. Darcet, T. Moutakanni, H. Vo, M. Szafraniec, V. Khalidov, P. Fernandez, D. Haziza, F. Massa, A. El-Nouby _et al._ , “Dinov2: Learning robust visual features without supervision,” _arXiv preprint arXiv:2304.07193_ , 2023.
  * [19] NASA, “NASA and IBM Openly Release Geospatial AI Foundation Model for NASA Earth Observation Data ,” https://www.earthdata.nasa.gov/news/impact-ibm-hls-foundation-model, Last Accessed on August 20, 2023.
  * [20] A. Karpatne, Z. Jiang, R. R. Vatsavai, S. Shekhar, and V. Kumar, “Monitoring land-cover changes: A machine-learning perspective,” _IEEE Geoscience and Remote Sensing Magazine_ , vol. 4, no. 2, pp. 8–21, 2016.
  * [21] T. Kurth, S. Subramanian, P. Harrington, J. Pathak, M. Mardani, D. Hall, A. Miele, K. Kashinath, and A. Anandkumar, “Fourcastnet: Accelerating global high-resolution weather forecasting using adaptive fourier neural operators,” in _Proceedings of the Platform for Advanced Scientific Computing Conference_ , 2023, pp. 1–11.
  * [22] N. Shi, J. Xu, S. W. Wurster, H. Guo, J. Woodring, L. P. Van Roekel, and H.-W. Shen, “Gnn-surrogate: A hierarchical and adaptive graph neural network for parameter space exploration of unstructured-mesh ocean simulations,” _IEEE Transactions on Visualization and Computer Graphics_ , vol. 28, no. 6, pp. 2301–2313, 2022.
  * [23] N. Kumar, G. Voulgaris, J. C. Warner, and M. Olabarrieta, “Implementation of the vortex force formalism in the coupled ocean-atmosphere-wave-sediment transport (coawst) modeling system for inner shelf and surf zone applications,” _Ocean Modelling_ , vol. 47, pp. 65–95, 2012.
  * [24] M. Olabarrieta, J. C. Warner, and N. Kumar, “Wave-current interaction in willapa bay,” _Journal of Geophysical Research: Oceans_ , vol. 116, no. C12, 2011.
  * [25] M. Olabarrieta, J. C. Warner, B. Armstrong, J. B. Zambon, and R. He, “Ocean–atmosphere dynamics during hurricane ida and nor’ida: An application of the coupled ocean–atmosphere–wave–sediment transport (coawst) modeling system,” _Ocean Modelling_ , vol. 43, pp. 112–137, 2012\.
  * [26] C.-E. Hsu, K. Serafin, X. Yu, C. Hegermiller, J. C. Warner, and M. Olabarrieta, “Total water levels along the south atlantic bight during three along-shelf propagating tropical cyclones: relative contributions of storm surge and wave runup,” _Natural Hazards and Earth System Sciences Discussions_ , vol. 2023, pp. 1–31, 2023.
  * [27] D. O. M. Group, “Regional Ocean Modeling System (ROMS),” https://www.myroms.org/, 2023.
  * [28] W. He, Z. Jiang, T. Xiao, Z. Xu, S. Chen, R. Fick, M. Medina, and C. Angelini, “A hierarchical spatial transformer for massive point samples in continuous space,” in _Proceedings of the 2023 Conference on Neural Information Processing Systems (NeurIPS)_ , 2023.
  * [29] A. Karpatne, R. Kannan, and V. Kumar, _Knowledge Guided Machine Learning: Accelerating Discovery Using Scientific Knowledge and Data_. CRC Press, 2022.
  * [30] Y. Gil and B. Selman, “A 20-year community roadmap for artificial intelligence research in the us,” _arXiv preprint arXiv:1908.02624_ , 2019.
  * [31] W. He and Z. Jiang, “A survey on uncertainty quantification methods for deep neural networks: An uncertainty source perspective,” _arXiv preprint arXiv:2302.13425_ , 2023.
  * [32] Z. Jiang, S. Shekhar, P. Mohan, J. Knight, and J. Corcoran, “Learning spatial decision tree for geographical classification: a summary of results,” in _Proceedings of the 20th International Conference on Advances in Geographic Information Systems_ , 2012, pp. 390–393.
  * [33] K. Han, Y. Wang, H. Chen, X. Chen, J. Guo, Z. Liu, Y. Tang, A. Xiao, C. Xu, Y. Xu _et al._ , “A survey on vision transformer,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 45, no. 1, pp. 87–110, 2022.
  * [34] Z.-H. Zhou, “A brief introduction to weakly supervised learning,” _National science review_ , vol. 5, no. 1, pp. 44–53, 2018.
  * [35] W. He, Z. Jiang, M. Kriby, Y. Xie, X. Jia, D. Yan, and Y. Zhou, “Quantifying and reducing registration uncertainty of spatial vector labels on earth imagery,” in _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_ , 2022, pp. 554–564.
  * [36] Z. Jiang, W. He, M. S. Kirby, A. M. Sainju, S. Wang, L. V. Stanislawski, E. J. Shavers, and E. L. Usery, “Weakly supervised spatial deep learning for earth image segmentation based on imperfect polyline labels,” _ACM Transactions on Intelligent Systems and Technology (TIST)_ , vol. 13, no. 2, pp. 1–20, 2022.
  * [37] Z. Xu, T. Xiao, W. He, Y. Wang, and Z. Jiang, “Spatial knowledge-infused hierarchical learning: An application in flood mapping on earth imagery,” in _The 31st ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems (GIS)_ , 2023.
  * [38] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” _Journal of Computational physics_ , vol. 378, pp. 686–707, 2019.
  * [39] N. B. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. M. Stuart, and A. Anandkumar, “Neural operator: Learning maps between function spaces with applications to pdes.” _J. Mach. Learn. Res._ , vol. 24, no. 89, pp. 1–97, 2023.
  * [40] Q. Sun, Y. Liu, T.-S. Chua, and B. Schiele, “Meta-transfer learning for few-shot learning,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2019, pp. 403–412.
  * [41] Z. Jiang, Y. Li, S. Shekhar, L. Rampi, and J. Knight, “Spatial ensemble learning for heterogeneous geographic data with class ambiguity: A summary of results,” in _Proceedings of the 25th ACM SIGSPATIAL international conference on advances in geographic information systems_ , 2017, pp. 1–10.
  * [42] Z. Jiang, A. M. Sainju, Y. Li, S. Shekhar, and J. Knight, “Spatial ensemble learning for heterogeneous geographic data with class ambiguity,” _ACM Transactions on Intelligent Systems and Technology (TIST)_ , vol. 10, no. 4, pp. 1–25, 2019.


# Computation of NMR shieldings at the CASSCF level using gauge-including
atomic orbitals and Cholesky decomposition

Tommaso Nottoli<EMAIL_ADDRESS>Dipartimento di Chimica e
Chimica Industriale, Università di Pisa, Via G. Moruzzi 13, I-56124 Pisa,
Italy  Sophia Burger<EMAIL_ADDRESS>Department Chemie, Johannes
Gutenberg-Universität Mainz, Duesbergweg 10-14, D-55128 Mainz, Germany  Stella
Stopkowicz<EMAIL_ADDRESS>Department Chemie, Johannes
Gutenberg-Universität Mainz, Duesbergweg 10-14, D-55128 Mainz, Germany
Fachrichtung Chemie, Universität des Saarlandes, Campus B2.2, D-66123
Saarbrücken, Germany  Jürgen Gauss<EMAIL_ADDRESS>Department Chemie,
Johannes Gutenberg-Universität Mainz, Duesbergweg 10-14, D-55128 Mainz,
Germany  Filippo Lipparini<EMAIL_ADDRESS>Dipartimento di Chimica
e Chimica Industriale, Università di Pisa, Via G. Moruzzi 13, I-56124 Pisa,
Italy

###### Abstract

We present an implementation of coupled-perturbed complete active space self-
consistent field (CP-CASSCF) theory for the computation of nuclear magnetic
resonance chemical shifts using gauge-including atomic orbitals and Cholesky
decomposed two-electron integrals. The CP-CASSCF equations are solved using a
direct algorithm where the magnetic Hessian matrix-vector product is expressed
in terms of one-index transformed quantities. Numerical tests on systems with
up to about 1300 basis functions provide information regarding both the
computational efficiency and limitations of our implementation.

††preprint: AIP/123-QED

##  I Introduction

Ab initio calculations of nuclear magnetic resonance (NMR) chemical shifts[1]
represent a powerful tool for the interpretation of experimental NMR spectra
and for elucidating chemical structures and conformations. It is well known
that the accurate prediction of NMR shieldings requires to take into account
dynamic correlation[2, 3, 4] most often through the use of the existing
hierarchy of coupled-cluster methods.[5] Nevertheless, even the use of such
accurate techniques may give erroneous results for systems that possess a
significant multi-reference character. An established starting point for the
treatment of multi-reference – also called strongly correlated – systems is
given by the complete active space–self-consistent field (CASSCF) method.[6,
7, 8, 9] With an adequate choice for the active space, this approach is able
to qualitatively describe transition-metal complexes, radicals, molecules with
stretched bonds and other systems where single-reference approaches fail.

During the eighties and nineties, multiconfigurational (MC) SCF linear-
response theory was developed[10, 11] and applied to the study of excitation
energies, transition moments, and other excited state properties.[12, 13] At
the same time, analytical second derivative formulations for a general MCSCF
wave function were developed[14, 15, 16] where perturbed molecular and
configurational parameters are obtained by solving a set of linear equations
called coupled-perturbed (CP). The solution of this linear system of equations
was later implemented using a direct formulation[17, 18, 19] where the
response Hessian matrix was never explicitly calculated. Concerning the
calculation of magnetic properties for MCSCF references through the use of
perturbation-dependent basis functions, namely gauge-including atomic orbitals
(GIAOs), also called London orbitals,[20, 21, 22, 23] implementations have
been reported for NMR chemical shifts, magnetizabilities and the related
rotational g-tensor.[24, 25, 26] In particular, the first description of the
GIAO-MCSCF theory for shieldings and its implementation exploiting direct
techniques was presented by Ruud et al.[24] It should be also noted that a
MCSCF scheme for the computation of shieldings with a different treatment of
the gauge-origin problem, that is by using individual gauges for localized
orbitals (IGLOs), was reported by van Wüllen and Kutzelnigg[27] in the same
years.

CASSCF calculations and the subsequent solution of the corresponding CP
equations are hampered by the combinatorial scaling exhibited by the full-CI
(FCI) problem in the active space. For this reason, standard direct FCI
implementations as the one reported in the present paper can be routinely
applied only to small/medium sized active spaces, that is, up to 14 electrons
in 14 orbitals. Under these circumstances, the main bottleneck in the
calculation is given by the storage and manipulation of the two-electron
repulsion integral (ERI) matrix, i.e., the operations needed for the orbital
optimization. The computational cost required to solve the CP-CASSCF equations
are similar to the one of a conventional second-order CASSCF optimization
algorithm. In particular, ERIs with two external indices are needed and the
computational scaling is $\mathcal{O}((N_{\rm int}+N_{\rm act})N^{4}_{\rm
b})$, $N_{\rm int}$ and $N_{\rm act}$ being the number of inactive and active
orbitals, respectively, and $N_{\rm b}$ the total number of basis functions.

In the last few years special emphasis has been placed on the development and
rewriting of efficient and fast quantum-chemistry codes that are able to
describe larger molecular systems. As for what concerns CASSCF second-order
properties, one of the first efforts in this direction was made by Dudley et
al.[28] who reported the implementation of a parallelized scheme for the
solution of the CP-MCSCF equations. More recently, Snyder et al.[29, 30, 31]
convincingly presented an implementation that exploits the sparsity in the
atomic orbital (AO) basis and parallelizes the calculation over graphical
processing units (GPU). In their paper, they report an observed reduction of
the computational scaling from $\mathcal{O}(N^{4}_{\rm b})$ to
$\mathcal{O}(N^{2}_{\rm b})$, thus allowing calculations on systems with up to
almost 7000 basis functions. Furthermore, Helmich-Paris[32] recently reported
results of linear-response CASSCF calculations on systems with more than 2000
basis functions. Despite being a linear-response formulation, such work shares
similarities with the present implementation since it implements a direct
algorithm based on the computation of one-index transformed Fock matrices, as
it will also be shown in our derivation. In order to reduce computational
cost, the work of Helmich-Paris exploits the resolution of the identity[33,
34, 35] and the overlap-fitted chain-of-spheres (COSX) approximation.[36, 37]

In this contribution, we exploit the Cholesky decomposition (CD) of the two-
electron integrals.[38, 39] This technique possesses the remarkable property
of giving approximate results that can be easily improved by lowering the
decomposition threshold. The formal scaling of the AO to MO integral
transformation is reduced and is no longer the bottleneck. Conversely, special
care has to be given to the contractions required to build Fock matrices and
intermediates. Most importantly, the memory requirement is highly reduced with
respect to a standard calculation. The CD has been widely used to accelerate
the computation of energies and gradients at various levels of theory,[40, 41,
42, 43, 44, 45, 46, 47] but so far little has been done for the computation of
second-order properties, especially magnetic ones. The first application of
the CD for the calculation of NMR shieldings has been only recently reported
by some of us within a GIAO based second-order Møller-Plesset perturbation
theory (GIAO-MP2) treatment.[48]

In this paper, we discuss the implementation of a CD-based GIAO-CASSCF theory
for the calculation of NMR chemical shifts within the CFour program
package.[49, 50] In section II, we briefly review the main equations of GIAO-
CASSCF theory and we describe the implementation of CD for the solution of the
CP-CASSCF equations focusing on the handling of perturbed Cholesky vectors. In
section III, we report our results showing both the accuracy of CD for the
computation of NMR shieldings and the efficiency of our implementation in
treating medium-sized systems. Section IV presents some concluding remarks.

##  II Theory and Implementation

In this section we briefly review GIAO-CASSCF theory outlining the main
equations. Then the Cholesky decomposition of perturbed integrals is
introduced and the implementation of CD-GIAO-CASSCF is discussed. The
following conventions for the indices are used in the following:

  * •

$\mu,\nu,\rho,\dots$ refer to AOs;

  * •

$p,q,r,\dots$ refer to generic MOs;

  * •

$i,j,k,\dots$ refer to inactive MOs;

  * •

$u,v,x,\dots$ refer to active MOs;

  * •

$a,b,c,\dots$ refer to external MOs;

  * •

$I,J,\dots$ refer to either state functions or Slater determinants;

  * •

$P,Q,\dots$ refer to the Cholesky vectors.

###  II.1 GIAO-CASSCF theory

The CASSCF wave function can be defined using the exponential unitary ansatz

$\ket{\Psi}=e^{i\hat{\kappa}}e^{i\hat{S}}\ket{0},$ (1)

where $\ket{0}$ is the multideterminantal reference state, $i$ the imaginary
unit, and $\hat{\kappa}$ and $\hat{S}$ are the Hermitian orbital and
configurational rotation operators defined as

$\displaystyle\hat{\kappa}=\sum_{p>q}\kappa_{pq}\hat{E}^{+}_{pq},$ (2)
$\displaystyle\hat{S}=\sum_{J\neq 0}S_{J}(\ket{J}\bra{0}+\ket{0}\bra{J}).$ (3)

Here, $\hat{E}^{+}_{pq}=\hat{E}_{pq}+\hat{E}_{qp}$ is the symmetric singlet
excitation operator, while $\ket{J}$ represents a generic state orthogonal to
the reference one. This functional form is more general than the one used in
Ref. 51 and 52, as we aim at describing the effect of a magnetic perturbation
that introduces imaginary terms into the Hamiltonian.[53] Furthermore, to
remove the gauge-origin dependence that affects the calculation of magnetic
properties with finite basis sets, we resort as usual to the use of GIAOs.[20,
21, 22, 23] This means that the basis set explicitly depends on the
perturbation, i.e., the external magnetic field ($\mathbf{B}$). To treat the
magnetic-field dependence of the AO basis set, we follow the formalism
proposed by Helgaker and Jørgensen[54] that transfers the AO dependence on the
perturbation to the MO overlap matrix. In other words, we require the MOs to
stay orthonormal for any value of the perturbing field. It can be shown that
this is achieved by introducing a modified gradient that enforces
orthonormality of the MOs:

$\nabla_{x}h=\partialderivative{h}{x}-\frac{1}{2}\left\\{\partialderivative{\ln{S}}{x},h\right\\}.$
(4)

In Eq. (4), curly braces are a compact representation for a one-index
transformation, e.g.

$\\{A,B\\}_{pq}=\sum_{r}\left(A_{pr}B_{rq}+A^{*}_{qr}B_{pr}\right).$ (5)

In the definition of the modified gradient in Eq. (4), we adopt the so-called
symmetric connection. An extended discussion on the topic of orbital
connections can be found in the literature.[55, 54, 56, 57]

All terms that stem from the second part of Eq. (4) can be viewed as
originating from the reorthogonalization of the MOs. As usual, differentiated
quantities are then evaluated in the limit of vanishing perturbation. In this
limit, the derivative of the logarithm of the overlap matrix can be shown to
be equal to the derivative of the overlap matrix itself. We also note that Eq.
(4) becomes the usual gradient when the overlap matrix does not depend on the
perturbation $x$. The NMR chemical shieldings can be defined via the
analytical second derivative of the variational CASSCF energy, $\mathcal{E}$,
with respect to the components of an external magnetic field and the nuclear
magnetic dipole moment:[24]

$\sigma_{ji}=\nabla_{B_{i}}\partialderivative{\mathcal{E}}{m_{N_{j}}}+\sum_{p>q}\partialderivative{\mathcal{E}}{m_{N_{j}}}{\kappa_{pq}}\derivative{\kappa_{pq}}{B_{i}}+\sum_{J\neq
0}\partialderivative{\mathcal{E}}{m_{N_{j}}}{S_{J}}\derivative{S_{J}}{B_{i}}.$
(6)

In Eq. (6), the first term is a static contribution, while the other two terms
stem from chain rule differentiation and include the wave function response to
the perturbation and will be called dynamic terms. Given the fact that the
GIAOs do not depend on the magnetic dipole moments, it is convenient to
perform first the differentiation with respect to $m_{N_{j}}$ and then with
respect to $B_{i}$ such that the static contribution takes the following form

$\nabla_{\mathbf{B}}\partialderivative{\mathcal{E}}{\mathbf{m}}=\sum_{pq}\left[h_{pq}^{(\mathbf{B},\mathbf{m})}-\frac{1}{2}\sum_{r}(S^{(\mathbf{B})}_{pr}h^{(\mathbf{m})}_{rq}+{S^{(\mathbf{B})}}^{*}_{qr}h^{(\mathbf{m})}_{pr})\right]\gamma_{pq},$
(7)

where the superscript between parentheses is a shorthand notation for the
explicit derivative of the integrals in the GIAO basis, e.g.

$h_{pq}^{(\mathbf{B},\mathbf{m})}=\sum_{\mu\nu}c^{*}_{\mu
p}\partialderivative{h_{\mu\nu}}{\mathbf{B}}{\mathbf{m}}c_{\nu q}$ (8)

with $c_{\mu p}$ as the MO coefficients. Explicit expressions for the
differentiated one-electron GIAO integrals are given by

$\displaystyle\partialderivative{h_{\mu\nu}}{m_{N_{j}}}=\frac{1}{c}\bra{\chi_{\mu}}\frac{L_{N_{j}}}{|\mathbf{r}-\mathbf{R}_{N}|^{3}}\ket{\chi_{\nu}},$
$\displaystyle\partialderivative{h_{\mu\nu}}{B_{i}}{m_{N_{j}}}=\frac{1}{2c^{2}}\bra{\chi_{\mu}}\frac{\mathbf{r}\cdot(\mathbf{r}-\mathbf{R}_{N})\delta_{ij}-\mathbf{r}_{i}(\mathbf{r}-\mathbf{R}_{N})}{|\mathbf{r}-\mathbf{R}_{N}|^{3}}\ket{\chi_{\nu}},$

$L_{N_{j}}$ being the $j$-th component of the angular-momentum operator around
nucleus $N$, $\mathbf{R}_{N}$ the position vector of nucleus $N$, and $c$ the
speed of light in atomic units. We use bold symbols to refer to the whole
vector or array and implicitly assume that all derivatives are evaluated at
zero field.

The calculation of the response terms of Eq. (6) requires the solution of the
coupled-perturbed CASSCF (CP-CASSCF) equations

$\mathbf{G}\derivative{\bm{\lambda}}{\mathbf{B}}=-\nabla_{\mathbf{B}}\partialderivative{\mathcal{E}}{\bm{\lambda}},$
(9)

where $\bm{\lambda}$ is a joint parameter including both $\kappa_{pq}$ and
$S_{J}$ and where
$\mathbf{G}=\partialderivative{\mathcal{E}}{\bm{\lambda}}{\bm{\lambda}}$ is
the magnetic CASSCF Hessian.

The CP-CASSCF equations are solved in a direct way without the explicit
computation and storage of the magnetic CASSCF Hessian. For NMR computations
we use a plain preconditioned conjugate gradient (PCG) algorithm being the
mathematically optimal choice for the solution of a symmetric, positive-
definite linear system of equations. The diagonal Hessian is used as
preconditioner.

The right hand side of Eq. (9) is the derivative of the CASSCF gradient
($\mathbf{g}=\partialderivative{\mathcal{E}}{\bm{\lambda}}$) with respect to
the external magnetic field, which can be written as follows:

$\displaystyle\nabla_{\mathbf{B}}g_{pq}=2i(\nabla_{\mathbf{B}}F_{pq}+\nabla_{\mathbf{B}}F_{qp}),$
(10)
$\displaystyle\nabla_{\mathbf{B}}g_{I}=2i\bra{\Phi_{I}}\sum_{uv}\hat{E}_{uv}\nabla_{\mathbf{B}}F^{I}_{uv}+\frac{1}{2}\sum_{uvxy}\hat{e}_{uvxy}\nabla_{\mathbf{B}}g_{uvxy}\ket{0}.$
(11)

Eq. (10) represents the molecular part of the gradient while Eq. (11)
represents the configurational one, with $\ket{\Phi_{I}}$ being a generic
Slater determinant. $g_{uvxy}$ is the active part of two-electron integral
matrix written using Mulliken’s convention. $F_{pq}$ is the generalized Fock
matrix that can be defined in terms of the inactive and active Fock matrices,
$F^{I}_{pq}$ and $F^{A}_{pq}$ respectively, and of the $Q$ matrix. The
derivatives of these three intermediates with respect to $\mathbf{B}$, using
the rule specified by Eq. (4), are given by

$\displaystyle\nabla_{\mathbf{B}}F^{I}_{pq}$ $\displaystyle=$ $\displaystyle
F^{I(\mathbf{B})}_{pq}-\frac{1}{2}\left\\{S^{(\mathbf{B})},F^{I}\right\\}_{pq}$
(12)
$\displaystyle-\frac{1}{2}\sum_{ir}\Big{(}g_{priq}-g_{pirq}\Big{)}S^{(\mathbf{B})}_{ir},$
$\displaystyle\nabla_{\mathbf{B}}F^{A}_{pq}$ $\displaystyle=$ $\displaystyle
F^{A(\mathbf{B})}_{pq}-\frac{1}{2}\left\\{S^{(\mathbf{B})},F^{A}\right\\}_{pq}$
(13)
$\displaystyle-\frac{1}{4}\sum_{uv}\sum_{r}\Big{(}g_{pruq}-g_{purq}\Big{)}S^{(\mathbf{B})}_{vr}\gamma_{uv},$
$\displaystyle\nabla_{\mathbf{B}}Q_{xp}$ $\displaystyle=$ $\displaystyle
Q^{(\mathbf{B})}_{xp}-\frac{1}{2}\sum_{r}S^{(\mathbf{B})}_{pr}Q_{xr}-\frac{1}{2}\sum_{uvy}\sum_{r}\Big{[}(\Gamma_{xuvy}-\Gamma_{xuyv})g_{pury}$
(14) $\displaystyle-\Gamma_{xvuy}g_{pruy}\Big{]}S^{(\mathbf{B})}_{vr}.$

We note that these equations consist of an explicit derivative term and
various reorthogonalization quantities involving the differentiated metric.
Eq. (11) does not include, unlike for real perturbations, the difference
between the perturbed reference and inactive energy since the trace of an
imaginary operator vanishes. This last fact holds also for the configurational
gradient perturbed with respect to nuclear magnetic dipole moments.

The implementation of the direct matrix-vector product between the magnetic
CASSCF Hessian and the imaginary perturbed parameters
($\derivative{\lambda}{B_{i}}$) exploits the same techniques as described in
Ref. 58 and 52, where all the relevant details can be found. We report in the
following the final expressions that need to be evaluated in order to assemble
the matrix-vector product between the Hessian and the perturbed parameters

$\displaystyle\sum_{r>s}G_{pq,rs}i\derivative{\kappa_{pq}}{B_{i}}=-2i(\tilde{F}_{pq}+\tilde{F}_{qp}),$
(15) $\displaystyle\sum_{I}G_{I,pq}i\derivative{c_{I}}{B_{i}}=2i({}^{\rm
T}F_{pq}+{}^{\rm T}F_{qp}),$ (16)
$\displaystyle\sum_{p>q}G_{pq,I}i\derivative{\kappa_{pq}}{B_{i}}=2i\bra{0}\tilde{\mathscr{H}}\ket{\Phi_{I}},$
(17)
$\displaystyle\sum_{J}G_{I,J}i\derivative{c_{I}}{B_{i}}=2i\left(\matrixelement{\Phi_{I}}{\mathscr{H}}{\derivative{\mathbf{c}}{B_{i}}}-E_{0}\derivative{c_{I}}{B_{i}}\right).$
(18)

In Eq. (15), we introduced the one-index transformed Fock matrices:

$\tilde{F}^{I}_{pq}=\left\\{\derivative{\kappa}{B_{i}},F^{I}\right\\}_{pq}+\sum_{ir}\Big{(}g_{priq}-g_{pirq}\Big{)}\derivative{\kappa_{ir}}{B_{i}},$
(19)
$\tilde{F}^{A}_{pq}=\left\\{\derivative{\kappa}{B_{i}},F^{A}\right\\}_{pq}+\frac{1}{2}\sum_{uvr}\Big{(}g_{pruq}-g_{purq}\Big{)}\gamma_{uv}\derivative{\kappa_{vr}}{B_{i}},$
(20) $\displaystyle\tilde{Q}_{xp}$ $\displaystyle=$
$\displaystyle\sum_{r}Q_{xr}\derivative{\kappa_{pr}}{B_{i}}$ (21)
$\displaystyle+\sum_{uvy}\sum_{r}\Big{[}(\Gamma_{xuvy}-\Gamma_{xuyv})g_{pury}-\Gamma_{xvuy}g_{pruy}\Big{]}\derivative{\kappa_{vr}}{B_{i}},$

that has a similar structure and symmetries as Eqs. (12) to (14). In Eq. (16),
we have introduced the generalized Fock matrix built with antisymmetrized
transition density matrices, that is

${}^{\rm T}F_{pq}=\sum_{r}{}^{\rm T}\gamma_{pr}h_{qr}+\sum_{rst}{}^{\rm
T}\Gamma_{prst}g_{qrst},$ (22)

with

$\displaystyle{}^{\rm
T}\gamma_{uv}=\matrixelement{0}{\hat{E}_{uv}}{\derivative{\mathbf{c}}{B_{i}}}-\matrixelement{\derivative{\mathbf{c}}{B_{i}}}{\hat{E}_{uv}}{0},$
(23) $\displaystyle{}^{\rm
T}\Gamma_{uvxy}=\matrixelement{0}{\hat{e}_{uvxy}}{\derivative{\mathbf{c}}{B_{i}}}-\matrixelement{\derivative{\mathbf{c}}{B_{i}}}{\hat{e}_{uvxy}}{0},$
(24)

and where
$\ket{\derivative{\mathbf{c}}{B_{i}}}=\sum_{I}\ket{\Phi_{I}}\derivative{c_{I}}{B_{i}}$.
Finally, in Eq. (17), $\tilde{\mathscr{H}}$ is the one-index transformed
Hamiltonian:

$\tilde{\mathscr{H}}=\sum_{uv}\hat{E}_{uv}\tilde{F}^{I}_{uv}+\frac{1}{2}\sum_{uvxy}\hat{e}_{uvxy}\tilde{g}_{uvxy}.$
(25)

The integrals referenced in Eqs. (11) and (17) show different permutational
symmetries compared to the real case. In particular, the one-electron
integrals are antisymmetric while the two-electron integrals only exhibit
fourfold permutational symmetry. The direct FCI routines used in our CASSCF
code have been generalized accordingly. We also note that Eqs. (15) to (18)
have been expressed in the determinantal basis – i.e., $I$ and $J$ are running
over all Slater determinants.

###  II.2 Implementation of GIAO-CASSCF with Cholesky-decomposed integrals

The two-step Cholesky decomposition algorithm proposed by Folkestad et
al.,[59] and further refined by Zhang et al.[60] has been implemented in the
Mainz INTegral (MINT) package[61] by some of us and extended for the
calculation of magnetic integral derivatives.[48, 62] The derivative of the
two-electron repulsion integral matrix can be written in a form similar to
corresponding density-fitting expressions

$\begin{split}\left(\partialderivative{(\mu\nu|\rho\sigma)}{\mathbf{B}}\right)_{\mathbf{B}=0}=&\sum_{PQ}^{N_{ch}}\left(\partialderivative{\mu\nu}{\mathbf{B}}\Big{|}P\right)_{\mathbf{B}=0}M_{PQ}^{-1}\left(Q|\rho\sigma\right)+\\\
&+\sum_{PQ}^{N_{ch}}\left(\mu\nu|P\right)M_{PQ}^{-1}\left(Q\Big{|}\partialderivative{\rho\sigma}{\mathbf{B}}\right)_{\mathbf{B}=0},\end{split}$
(26)

$M^{-1}_{PQ}$ being the inverse metric matrix in the non-orthogonal Cholesky
basis. From Eq. (26), the perturbed and unperturbed Cholesky vectors can be
defined respectively as

$\displaystyle\left(\partialderivative{L^{P}_{\mu\nu}}{\mathbf{B}}\right)_{\mathbf{B}=0}=\sum_{Q}\left(\partialderivative{\mu\nu}{\mathbf{B}}\Big{|}Q\right)_{\mathbf{B}=0}K_{QP}^{-T},$
(27) $\displaystyle L^{P}_{\mu\nu}=\sum_{Q}\left(\mu\nu|Q\right)K_{QP}^{-T},$
(28)

with $\mathbf{K}^{-T}$ – a shorthand notation for $(\mathbf{K}^{-1})^{T}$ – as
the inverse Cholesky factor of the metric, i.e.,
$\mathbf{M}=\mathbf{K}^{T}\mathbf{K}$. As shown in Ref. 48 and further
demonstrated in Ref. 62, the perturbed Cholesky vectors constitute
antisymmetric matrices and hence only the lower triangular part is stored.

All the equations presented in the previous section have been implemented
exploiting Cholesky-decomposed integrals and using the techniques previously
described in Ref. 52. In particular, the loop over the Cholesky vectors is
kept external in the code and is parallelized by means of shared-memory OpenMP
instructions. Furthermore, for each Cholesky vector, efficient level 2 and 3
BLAS routines are used to perform various contractions. As an example, let us
illustrate the implementation of the first terms of Eqs. (12), (13), and (14).
In particular, $F^{I(\mathbf{B})}_{pq}$ and $F^{A(\mathbf{B})}_{pq}$ are built
in the atomic-orbital basis.

$\displaystyle F^{I(\mathbf{B})}_{\mu\nu}$ $\displaystyle=$ $\displaystyle
h^{(\mathbf{B})}_{\mu\nu}+\sum_{P}\sum_{\rho\sigma}D^{I}_{\rho\sigma}\Bigg{[}\partialderivative{L^{P}_{\mu\nu}}{\mathbf{B}}L^{P}_{\rho\sigma}-L^{P}_{\mu\nu}\partialderivative{L^{P}_{\rho\sigma}}{\mathbf{B}}$
(29)
$\displaystyle-\frac{1}{2}\partialderivative{L^{P}_{\mu\rho}}{\mathbf{B}}L^{P}_{\sigma\nu}+\frac{1}{2}L^{P}_{\mu\rho}\partialderivative{L^{P}_{\sigma\nu}}{\mathbf{B}}\Bigg{]},$
$\displaystyle F^{A(\mathbf{B})}_{\mu\nu}$ $\displaystyle=$
$\displaystyle\sum_{P}\sum_{\rho\sigma}D^{A}_{\rho\sigma}\Bigg{[}\partialderivative{L^{P}_{\mu\nu}}{\mathbf{B}}L^{P}_{\rho\sigma}-L^{P}_{\mu\nu}\partialderivative{L^{P}_{\rho\sigma}}{\mathbf{B}}$
(30)
$\displaystyle-\frac{1}{4}\partialderivative{L^{P}_{\mu\rho}}{\mathbf{B}}L^{P}_{\sigma\nu}+\frac{1}{4}L^{P}_{\mu\rho}\partialderivative{L^{P}_{\sigma\nu}}{\mathbf{B}}\Bigg{]},$

with $D^{I}_{\mu\nu}=2\sum_{i}c^{*}_{\mu i}c_{\nu i}$ and
$D^{A}_{\mu\nu}=\sum_{uv}\gamma_{uv}c^{*}_{\mu u}c_{\nu v}$ being the AO
inactive and active density matrices, respectively. The Coulomb contribution
is easily computed by first contracting the density matrix with either a
perturbed or an unperturbed Cholesky vector. Concerning the exchange term, as
usual, it is assembled by first half transforming both the perturbed and
unperturbed vectors to the MO basis and then contracting them between each
other. Therefore, the formal scaling for the computation of the Coulomb
contributions is equal to $\mathcal{O}(N_{ch}N_{\rm b}^{2})$. The exchange
part of $F^{I(\mathbf{B})}_{\mu\nu}$ and $F^{A(\mathbf{B})}_{\mu\nu}$ requires
instead $\mathcal{O}(N_{\rm ch}N_{\rm b}^{2}N_{\rm act}+N_{\rm ch}N_{\rm
b}N^{2}_{\rm act})$ and $\mathcal{O}(N_{\rm ch}N_{\rm b}^{2}N_{int})$
floating-point operations, respectively.

The calculation of $Q^{(\mathbf{B})}_{xp}$ is easily done in the MO basis:

$Q^{(\mathbf{B})}_{xp}=\sum_{P}\sum_{uvy}\Gamma_{xuvy}\left(L^{P}_{pu}\partialderivative{L^{P}_{vy}}{\mathbf{B}}-\partialderivative{L^{P}_{pu}}{\mathbf{B}}L^{P}_{vy}\right),$
(31)

but we need here the perturbed Cholesky vectors in the MO basis with an active
index and a generic one. The cost for this transformation is
$\mathcal{O}(N_{\rm ch}N^{2}_{\rm b}N_{\rm act})$, which of course has to be
done for the three component of the magnetic field. Using perturbed Cholesky
vectors in the MO representation, the computational scaling of
$Q^{(\mathbf{B})}_{xp}$ is $\mathcal{O}(N_{\rm ch}N_{\rm act}^{4}+N_{\rm
ch}N_{\rm b}N_{\rm act}^{2})$. Therefore, the evaluation of Eqs. (29) to (31)
shows the same scaling as the construction of the unperturbed CD-based Fock
matrices with an additional prefactor of two (c.f. Eq. (26)). The active-
active block of the perturbed Cholesky vectors in the MO representation is
required also to build $\partialderivative{g_{uvxy}}{\mathbf{B}}$ for the
evaluation of Eq. (11).

The implementation has been carried out in the CFour program suite.[49, 50]
The workflow for the computation of CASSCF NMR shieldings is the following:
First, the MINT program is called where conventional one-electron integrals
and both the perturbed and unperturbed Cholesky vectors are calculated and
stored on disk. Second, the SCF equations are solved, and if requested, a
stability analysis and the computation of Unrestricted Natural Orbitals[63,
64] (UNOs) as a guess for the subsequent CASSCF calculation is performed.
Third, the CASSCF calculation is done. Fourth, MINT is called again for the
calculation of one-electron GIAO integrals, building differentiated Fock
matrices, and computing the static contribution to the chemical shieldings.
Eventually, the CP-CASSCF program is called for the solution of the coupled-
perturbed CASSCF equations and the evaluation of the dynamic terms. The
workflow is depicted in Fig. 1 for ease of reference.

Figure 1: Workflow of a CD-GIAO-CASSCF calculation in CFour.

##  III Numerical Results

###  III.1 Accuracy of GIAO-CASSCF NMR shieldings computed with Cholesky-
decomposed integrals

As a first analysis we tested the accuracy of isotropic nuclear magnetic
shieldings computed with Cholesky-decomposed integrals. The calculations have
been carried out on two structural isomers, namely formamide and formaldoxime
– whose molecular representation is shown in Fig. 2. In Table 1, we report the
maximum absolute error for the four different types of nuclei between a
reference calculation – i.e., where conventional integrals were used – and one
that exploits the CD with thresholds chosen to be $10^{-4}$ and $10^{-5}$,
respectively. To evaluate the effects of the CD approximation with different
basis sets, we perform these calculations with increasingly large basis sets
from Dunning’s hierarchy of correlation-consistent functions.[65] The
geometries have been optimized at the CASSCF(6,5)/cc-pVDZ level of theory
(using conventional integrals) and can be found in the supporting material
together with the resulting isotropic nuclear magnetic shieldings.

Figure 2:  The two molecules used to compare NMR shieldings between an exact calculation and one that uses Cholesky decomposed two-electron integrals. On the left, formamide, on the right formaldoxime. Table 1:  Maximum absolute error (in ppm) for the isotropic nuclear magnetic shieldings of [13]C, [17]O, [15]N, and [1]H in the two structural isomers formamide and formaldoxime. We report the result of various calculations performed with two different Cholesky-decomposition thresholds ($\delta$) and an increasing number of basis functions. |  $\delta$ = $10^{-4}$ |  $\delta$=$10^{-5}$
---|---|---
Nucleus | cc-pVDZ | cc-pVTZ | cc-pVQZ | cc-pVDZ | cc-pVTZ | cc-pVQZ
[13]C | 0.017 | 0.011 | 0.001 | 0.002 | 0.000 | 0.000
[17]O | 0.074 | 0.009 | 0.009 | 0.005 | 0.002 | 0.000
[15]N | 0.061 | 0.005 | 0.002 | 0.005 | 0.001 | 0.000
[1]H | 0.001 | 0.003 | 0.001 | 0.000 | 0.000 | 0.000

The order of magnitude of the errors is consistent with what has been
presented by Burger et al.;[48] i.e., using a Cholesky-decomposition threshold
equal to $10^{-4}$, the isotropic shieldings exhibit errors of no more than a
few hundredths of ppm. When using higher thresholds, as expected, the error
decreases. Specifically, with a threshold equal to $10^{-5}$, the errors
observed are of the order of $10^{-3}$ ppm. In all cases, the maximum errors
are of no relevance for the computational predictions of NMR shieldings.
Interestingly, it seems that the error is smaller when using larger basis
sets, which is consistent with an observed better agreement between the CASSCF
and CD-CASSCF energies in these cases.

###  III.2 Benchmark calculations on small- and medium-sized systems

In order to test the performance of the CD-GIAO-CASSCF algorithm we run
benchmark calculations exploiting the same set of molecules already used in
Ref. 52. We used a Cholesky-decomposition threshold of $10^{-5}$. The CASSCF
equations are considered converged when the root-mean-square (RMS) norm of the
configurational and orbital gradient is below $10^{-10}$, while for the CP-
CASSCF equations we used a threshold of $10^{-8}$ in the RMS of the residual.
These convergence criteria can be considered conservative and are more than
sufficient to obtain converged values for the isotropic NMR chemical shifts
within the error discussed in the previous subsection. As starting orbitals
for the CASSCF calculation we used the UNOs; i.e., we look for unstable
solution(s) in the RHF wave function, we optimize all the broken-symmetry UHF
solution(s), and we get the averaged density matrix from which we can compute
the natural orbitals. The results are shown in Table 2 and the isotropic
nuclear magnetic shieldings are reported in the supporting material.

Molecule | active space | $N_{\rm b}$ | 1st MINT | SCF111It includes also the timing required to search instabilities and solve all the UHF unstable solutions. The superscript numbers between parentheses represent the number of instabilities that were found and optimized. | CASSCF | 2nd MINT | CP-CASSCF | Total
---|---|---|---|---|---|---|---|---
adrenaline | 6,6 | 572 | 52.34 | 1.74(1) | 2.06 | 1.96 | 4.25 | 62.36
anthracene | 14,14 | 560 | 52.58 | 2.80(2) | 14.67 | 1.32 | 27.65 | 99.02
azulene | 10,10 | 412 | 21.97 | 0.45(1) | 0.72 | 0.70 | 1.38 | 25.21
biphenyl | 12,12 | 500 | 37.46 | 0.81(1) | 2.18 | 1.22 | 3.77 | 45.43
catechol | 6,6 | 324 | 11.49 | 0.21(1) | 0.28 | 0.38 | 0.40 | 12.75
dopamine | 6,6 | 484 | 33.14 | 0.79(1) | 1.40 | 1.10 | 1.78 | 38.2
fluorene | 12,12 | 530 | 45.72 | 1.32(1) | 2.56 | 2.20 | 3.96 | 55.78
indole | 8,8 | 368 | 16.47 | 0.35(1) | 0.48 | 0.42 | 0.77 | 18.49
l-dopamine | 6,6 | 574 | 44.9 | 2.08(1) | 3.28 | 2.64 | 4.40 | 57.31
naphthalene | 10,10 | 412 | 23.92 | 0.51(1) | 0.84 | 0.73 | 1.24 | 27.25
niacin | 6,6 | 340 | 12.99 | 0.29(1) | 0.35 | 0.41 | 0.52 | 14.55
niacinamide | 6,6 | 354 | 14.47 | 0.33(1) | 0.36 | 0.45 | 0.62 | 16.23
nicotine | 6,6 | 556 | 50.79 | 1.39(1) | 1.79 | 1.95 | 3.28 | 59.20
nor-adrenaline | 6,6 | 514 | 39.35 | 1.27(1) | 1.55 | 1.45 | 2.55 | 46.18
picolinic acid | 6,6 | 340 | 12.84 | 0.28(1) | 0.34 | 0.41 | 0.49 | 14.36
pyridine | 6,6 | 250 | 5.22 | 0.10(1) | 0.12 | 0.14 | 0.17 | 5.75
pyridoxal | 8,8 | 486 | 34.09 | 0.88(1) | 1.29 | 1.22 | 2.49 | 39.97
pyridoxamine | 6,6 | 528 | 42.72 | 1.21(1) | 1.82 | 1.47 | 2.9 | 50.12
pyridoxin | 6,6 | 514 | 39.97 | 1.09(1) | 1.48 | 1.45 | 2.84 | 46.83
resveratrol | 14,14 | 678 | 74.71 | 5.96(2) | 19.89 | 2.89 | 34.79 | 138.24
serotonin | 8,8 | 558 | 48.28 | 1.61(1) | 2.41 | 1.99 | 4.36 | 58.65
tryptophan | 8,8 | 618 | 64.02 | 2.40(1) | 3.58 | 2.38 | 6.87 | 79.24
2Me2HSdiox | 4,4 | 474 | 34.40 | 0.79(1) | 1.20 | 1.15 | 1.80 | 39.35
2Me4HSdiox | 6,6 | 446 | 30.34 | 0.64(1) | 0.91 | 0.96 | 1.39 | 34.24
coumarin dye | 12,12 | 872 | 149.86 | 13.14(2) | 12.35 | 7.77 | 21.14 | 204.27
BODIPY dye | 8,8 | 1280 | 386.94 | 77.53(3) | 96.05 | 35.57 | 71.96 | 668.06
Table 2: CD-GIAO-CASSCF computational timings for all the systems in the
benchmark set. For each molecule we report the active space, the total number
of basis functions ($N_{\rm b}$), and the CPU wall time (in minutes) for the
various programs required to compute NMR shieldings at the CD-GIAO-CASSCF
level. In particular, the time needed to compute one-electron integrals,
perturbed and unperturbed Cholesky vectors (1st MINT); the solution of the SCF
equations and the optimization of UHF unstable solution(s) (SCF); the
optimization of the CASSCF wave function (CASSCF); the computation of one-
electron GIAO integrals, $F^{I(\mathbf{B})}$ and $F^{A(\mathbf{B})}$, and the
evaluation of the static term (2nd MINT); the solution of the CP-CASSCF
equations and evaluation of the dynamic terms (CP-CASSCF). In the last column,
the overall time is also reported. All calculations were carried out on a
single cluster node equipped with 4 Intel Xeon Gold 6140 M CPUs running at
2.30 GHz and sharing the work over 28 OpenMP threads.

At the moment, we have not implemented the out-of-core handling of the
(un)perturbed Cholesky vectors. A NMR calculation on chlorophyll, one of the
benchmark systems presented in Ref. 52 is, however, only possible with such an
out-of-core implementation. We therefore replaced chlorophyll with the smaller
BODIPY dye, the structure of which is depicted in Fig. 3.

Figure 3:  New BODIPY dye molecule added to the benchmark set. Oxygen in red,
carbon in cyan, nitrogen in blue, hydrogen in white, boron in pink, and
fluorine in yellow.

As we did for the other systems in the set, we optimized the new structure at
the B3LYP/6-31G(d)[66, 67] level of theory using the Gaussian 16 software;[68]
the optimized geometry is given in the supporting material. We report the
timings for each step required for the CD-CASSCF calculation of nuclear
magnetic shieldings in Table 2. Comparing the timings between CASSCF and CP-
CASSCF, we notice that the two operations require about the same amount of
time. This is to be expected since the two methods exhibit the same scaling
and overall computational cost. More specifically, a PCG iteration in CP-
CASSCF is about as expensive as a micro-iteration in a norm-extended
optimization[69] (NEO) second-order optimization step. The ratio between the
cost of the two operations depends on the system, and is mainly influenced by
the convergence rate of the NEO macro- and micro-iterations on the CASSCF
side, which in turn depends on the quality of the starting guess, and on the
convergence of the PCG ones in the CP-CASSCF part of the calculation.
Moreover, we note that the CP-CASSCF problem is linear, but it has to be
solved for the three components of the magnetic field. It is important to
remark here the very good performance of the UNO guess, that affords robust
and smooth convergence of the nonlinear CASSCF optimization problem. This
comes, of course, at a price, as the calculation of the symmetry-broken UHF
solutions can be as expensive as the actual CASSCF calculation. Nevertheless,
besides the good quality of the guess, using the UNO approach provides
additional advantages, namely, a black-box procedure to select the active
space, and insights on the multireference character of the system, and in
particular on whether a multireference treatment is indeed warranted.
Nevertheless, when the UNO procedure becomes too expensive, one could resort
to one of the many other alternative approaches to select active spaces that
have been proposed in the literature.[70, 71, 72]

A more comprehensive analysis of the computational cost of the overall
calculation shows, however, that the real bottleneck in our implementation is
neither the CASSCF nor the CP-CASSCF step. In fact, the most expensive step,
at the current stage, is the calculation of the unperturbed and perturbed
Cholesky vectors, as performed by the first MINT call. With respect to our
previous work,[48] we are using a more efficient implementation of the
Cholesky decomposition, which is based on the two-step algorithm proposed by
Folkestad et al. [59] and which was extended to the calculation of the
perturbed vectors. Looking more in detail into the cost of this step, we
observed that it is strongly dominated by the calculation of the two-electron
integrals and integral derivatives, the former of which are required in both
steps of the CD calculation – i.e., both to determine the Cholesky basis and
to compute the actual Cholesky vectors. As discussed in Ref. 62, there is room
for further optimization of the integrals evaluation. Most importantly, this
part of the calculation is, at the moment, not parallelized: a parallel
evaluation of the shell-quartets required in the procedure would dramatically
reduce the overall cost of this step. An optimized, parallel implementation of
the evaluation of the Cholesky vectors and their derivatives is being actively
pursued. The only systems for which the integral evaluation is not completely
dominating the calculation are resveratrol, anthracene, and the BODIPY dye,
where we observe that the first MINT call takes about 60% of the overall
computer time, to be compared with 80% or more for the other calculations.
Concerning the first two systems, that are described using a larger (14,14)
active space, this can be explained by the cost stemming from the
configurational part of the CASSCF and CP-CASSCF solvers, which starts to
become non-negligible. For the BODIPY dye, on the other hand, the molecular
part is dominating the calculation, as it can be seen by the non-negligible
impact of the UNO procedure (as well represented by the SCF timing reported in
Table 2) as well as of the cost of the CASSCF and CP-CASSCF steps.

##  IV Conclusions

In this contribution we presented the implementation of a CASSCF procedure to
compute NMR chemical shifts that uses gauge-including atomic orbitals (GIAO)
and the Cholesky decomposition (CD) of the two-electron integral matrix. The
implementation closely resembles the one reported in Ref. 52 and follows the
direct formulation proposed by Jensen et al.[58] The computational cost for
the solution of the coupled-perturbed CASSCF (CP-CASSCF) equations is
asymptotically equivalent to the cost of a conventional second-order CASSCF
optimization algorithm. In particular, when using small active spaces, the
expensive operations are the contractions required to build the Fock matrices.
We showed that in the overall workflow for the computation of NMR shieldings
at the CD-GIAO-CASSCF level, the bottleneck consists in the calculation of the
perturbed and unperturbed Cholesky vectors and in particular the evaluation of
the two-electron integrals and the differentiated two-electron integrals. For
this reason, special effort will be invested in the future in the optimization
and parallelization of this part of the code.

At the moment, our implementation relies on the in-core storage of the
Cholesky vectors and, for the calculation of the property gradients, both the
unperturbed Cholesky vectors and one (out of three) set of perturbed vectors.
This is reasonable on standard computer nodes for calculations with up to
about 1500 basis functions. Larger calculations would require an out-of-core
handling of the Cholesky vectors. We are currently working on an efficient
implementation that streamlines and minimizes as much as possible slow disk
I/O.

As a future development, we are actively working on the implementation of
geometrical gradients and magnetizabilities, using doubly differentiated
Cholesky vectors, both at the CASSCF level of theory and at other correlated,
single-reference levels.

## Supplementary Material

See supplementary material for geometries and NMR isotropic shieldings of the
systems used to test the accuracy of CD-GIAO-CASSCF and the performance of the
overall algorithm.

###### Acknowledgements.

S.S. and J.G. acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG)
within project B5 of the TRR 146 (project no. 233 630 050). S.S. also
acknowledges support from the DFG via grant STO 1239/1-1.

## Data Availability Statement

The data that support the findings of this study are available within the
article and its supplementary material.

## References

  * Kaupp, Bühl, and Malkin [2004] M. Kaupp, M. Bühl, and V. G. Malkin,  _Calculation of NMR and EPR Parameters: Theory and Applications_ (Wiley-VCH, Weinheim, 2004).
  * Gauss [1993] J. Gauss, “Effects of electron correlation in the calculation of nuclear magnetic resonance chemical shifts,” J. Chem. Phys. 99, 3629–3643 (1993).
  * Gauss [1995] J. Gauss, “Accurate calculation of NMR chemical shifts,” Ber. Bunsenges. Phys. Chem. 99, 1001–1008 (1995).
  * Gauss and Stanton [2003] J. Gauss and J. F. Stanton, “Electron-correlated approaches for the calculation of NMR chemical shifts,” Adv. Chem. Phys 123, 355–422 (2003).
  * Shavitt and Bartlett [2009] I. Shavitt and R. J. Bartlett,  _Many-body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory_ (Cambridge University Press, 2009).
  * Roos, Taylor, and Siegbahn [1980] B. O. Roos, P. R. Taylor, and P. E. M. Siegbahn, “A complete active space SCF method (CASSCF) using a density matrix formulated super-CI approach,” Chem. Phys. 48, 157–173 (1980).
  * Werner [1987] H.-J. Werner, “Matrix-formulated direct multiconfiguration self-consistent field and multiconfiguration reference configuration-interaction methods,” Adv. Chem. Phys. 69, 1–62 (1987).
  * Shepard [1987] R. Shepard, “The multiconfiguration self-consistent field method,” Adv. Chem. Phys. 69, 63–200 (1987).
  * Roos [1987] B. O. Roos, “The complete active space self-consistent field method and its applications in electronic structure calculations,” Adv. Chem. Phys. 69, 399–445 (1987).
  * Dalgaard [1980] E. Dalgaard, “Time-dependent multiconfigurational Hartree–Fock theory,” J. Chem. Phys. 72, 816–823 (1980).
  * Yeager and Jørgensen [1979] D. L. Yeager and P. Jørgensen, “A multiconfigurational time-dependent Hartree-Fock approach,” Chem. Phys. Lett. 65, 77–80 (1979).
  * Olsen and Jørgensen [1985] J. Olsen and P. Jørgensen, “Linear and nonlinear response functions for an exact state and for an MCSCF state,” J. Chem. Phys. 82, 3235–3264 (1985).
  * Jørgensen, Jensen, and Olsen [1988] P. Jørgensen, H. J. Aa. Jensen, and J. Olsen, “Linear response calculations for large scale multiconfiguration self-consistent field wave functions,” J. Chem. Phys. 89, 3654–3661 (1988).
  * Jørgensen and Simons [1983] P. Jørgensen and J. Simons, “Ab initio analytical molecular gradients and Hessians,” J. Chem. Phys. 79, 334–357 (1983).
  * Camp  _et al._ [1983] R. N. Camp, H. F. King, J. W. McIver Jr., and D. Mullally, “Analytical force constants for MCSCF wave functions,” J. Chem. Phys. 79, 1088–1089 (1983).
  * Hoffmann  _et al._ [1984] M. R. Hoffmann, D. J. Fox, J. F. Gaw, Y. Osamura, Y. Yamaguchi, R. S. Grev, G. Fitzgerald, H. F. Schaefer III, P. J. Knowles, and N. C. Handy, “Analytic energy second derivatives for general MCSCF wave functions,” J. Chem. Phys. 80, 2660–2668 (1984).
  * Almlöf and Taylor [1985] J. Almlöf and P. R. Taylor, “Molecular properties from perturbation theory: A unified treatment of energy derivatives,” Int. J. Quantum Chem. 27, 743–768 (1985).
  * Yamamoto  _et al._ [1996] N. Yamamoto, T. Vreven, M. A. Robb, M. J. Frisch, and H. B. Schlegel, “A direct derivative MC-SCF procedure,” Chem. Phys. Lett. 250, 373–378 (1996).
  * Bernhardsson  _et al._ [1999] A. Bernhardsson, R. Lindh, J. Olsen, and M. Fülscher, “A direct implementation of the second-order derivatives of multiconfigurational SCF energies and an analysis of the preconditioning in the associated response equation,” Mol. Phys. 96, 617–628 (1999).
  * London [1937] F. London, “Théorie quantique des courants interatomiques dans les combinaisons aromatiques,” J. Phys. Radium 8, 397–409 (1937).
  * Hameka [1958] H. F. Hameka, “On the nuclear magnetic shielding in the hydrogen molecule,” Mol. Phys. 1, 203–215 (1958).
  * Ditchfield [1972] R. Ditchfield, “Molecular orbital theory of magnetic shielding and magnetic susceptibility,” J. Chem. Phys. 56, 5688–5691 (1972).
  * Wolinski, Hinton, and Pulay [1990] K. Wolinski, J. F. Hinton, and P. Pulay, “Efficient implementation of the gauge-independent atomic orbital method for NMR chemical shift calculations,” J. Am. Chem. Soc. 112, 8251–8260 (1990).
  * Ruud  _et al._ [1994] K. Ruud, T. Helgaker, R. Kobayashi, P. Jørgensen, K. L. Bak, and H. J. Aa. Jensen, “Multiconfigurational self-consistent field calculations of nuclear shieldings using London atomic orbitals,” J. Chem. Phys. 100, 8178–8185 (1994).
  * Ruud  _et al._ [1995a] K. Ruud, T. Helgaker, K. L. Bak, P. Jørgensen, and J. Olsen, “Accurate magnetizabilities of the isoelectronic series BeH-, BH, and CH+. The MCSCF-GIAO approach,” Chem. Phys. 195, 157–169 (1995a).
  * Ruud, Helgaker, and Jørgensen [1997] K. Ruud, T. Helgaker, and P. Jørgensen, “The effect of correlation on molecular magnetizabilities and rotational g tensors,” J. Chem. Phys. 107, 10599–10606 (1997).
  * van Wüllen and Kutzelnigg [1993] C. van Wüllen and W. Kutzelnigg, “The MC-IGLO method,” Chem. Phys. Lett. 205, 563–571 (1993).
  * Dudley  _et al._ [2006] T. J. Dudley, R. M. Olson, M. W. Schmidt, and M. S. Gordon, “Parallel coupled perturbed CASSCF equations and analytic CASSCF second derivatives,” J. Comp. Chem. 27, 352–362 (2006).
  * Snyder Jr.  _et al._ [2015a] J. W. Snyder Jr., E. G. Hohenstein, N. Luehr, and T. J. Martínez, “An atomic orbital-based formulation of analytical gradients and nonadiabatic coupling vector elements for the state-averaged complete active space self-consistent field method on graphical processing units,” J. Chem. Phys. 143, 154107 (2015a).
  * Snyder Jr.  _et al._ [2015b] J. W. Snyder Jr., E. G. Hohenstein, N. Luehr, and T. J. Martínez, “An atomic orbital-based formulation of analytical gradients and nonadiabatic coupling vector elements for the state-averaged complete active space self-consistent field method on graphical processing units,” J. Chem. Phys. 143, 154107 (2015b).
  * Snyder Jr.  _et al._ [2017] J. W. Snyder Jr., B. S. Fales, E. G. Hohenstein, B. G. Levine, and T. J. Martínez, “A direct-compatible formulation of the coupled perturbed complete active space self-consistent field equations on graphical processing units,” J. Chem. Phys. 146, 174113 (2017).
  * Helmich-Paris [2019] B. Helmich-Paris, “CASSCF linear response calculations for large open-shell molecules,” J. Chem. Phys. 150, 174121 (2019).
  * Whitten [1973] J. L. Whitten, “Coulombic potential energy integrals and approximations,” J. Chem. Phys. 58, 4496–4501 (1973).
  * Vahtras, Almlöf, and Feyereisen [1993] O. Vahtras, J. Almlöf, and M. W. Feyereisen, “Integral approximations for LCAO-SCF calculations,” Chem. Phys. Lett. 213, 514–518 (1993).
  * Neese [2003] F. Neese, “An improvement of the resolution of the identity approximation for the formation of the Coulomb matrix,” J. Comp. Chem. 24, 1740–1747 (2003).
  * Neese  _et al._ [2009] F. Neese, F. Wennmohs, A. Hansen, and U. Becker, “Efficient, approximate and parallel Hartree–Fock and hybrid DFT calculations. a ‘chain-of-spheres’ algorithm for the Hartree–Fock exchange,” Chem. Phys. 356, 98–109 (2009).
  * Izsák and Neese [2011] R. Izsák and F. Neese, “An overlap fitted chain of spheres exchange method,” J. Chem. Phys. 135, 144105 (2011).
  * Beebe and Linderberg [1977] N. H. F. Beebe and J. Linderberg, “Simplifications in the generation and transformation of two-electron integrals in molecular calculations,” Int. J. Quantum Chem. 12, 683–705 (1977).
  * Koch, Sánchez de Merás, and Pedersen [2003] H. Koch, A. Sánchez de Merás, and T. B. Pedersen, “Reduced scaling in electronic structure calculations using Cholesky decompositions,” J. Chem. Phys. 118, 9481–9484 (2003).
  * Bozkaya [2014] U. Bozkaya, “Derivation of general analytic gradient expressions for density-fitted post-Hartree-Fock methods: An efficient implementation for the density-fitted second-order Møller-Plesset perturbation theory,” J. Chem. Phys. 141, 124108 (2014).
  * Aquilante  _et al._ [2008] F. Aquilante, T. B. Pedersen, R. Lindh, B. O. Roos, A. Sánchez de Merás, and H. Koch, “Accurate ab initio density fitting for multiconfigurational self-consistent field methods,” J. Chem. Phys. 129, 024113 (2008).
  * Bozkaya [2016] U. Bozkaya, “Orbital-optimized MP3 and MP2.5 with density-fitting and Cholesky decomposition approximations,” J. Chem. Theory Comput. 12, 1179–1188 (2016).
  * Blaschke and Stopkowicz [2021] S. Blaschke and S. Stopkowicz, “Cholesky decomposition of complex two-electron integrals over GIAOs: Efficient MP2 computations for large molecules in strong magnetic fields,” J. Chem. Phys. 156, 044115 (2021).
  * Boström  _et al._ [2014] J. Boström, V. Veryazov, F. Aquilante, T. B. Pedersen, and R. Lindh, “Analytical gradients of the second-order Møller–Plesset energy using Cholesky decompositions,” Int. J. Quantum Chem. 114, 321–327 (2014).
  * Delcey  _et al._ [2014] M. G. Delcey, L. Freitag, T. B. Pedersen, F. Aquilante, R. Lindh, and L. González, “Analytical gradients of complete active space self-consistent field energies using Cholesky decomposition: Geometry optimization and spin-state energetics of a ruthenium nitrosyl complex,” J. Chem. Phys. 140, 174103 (2014).
  * Feng  _et al._ [2019] X. Feng, E. Epifanovsky, J. Gauss, and A. I. Krylov, “Implementation of analytic gradients for CCSD and EOM-CCSD using Cholesky decomposition of the electron-repulsion integrals and their derivatives: Theory and benchmarks,” J. Chem. Phys. 151, 014110 (2019).
  * Schnack-Petersen  _et al._ [2022] A. K. Schnack-Petersen, H. Koch, S. Coriani, and E. F. Kjønstad, “Efficient implementation of molecular CCSD gradients with Cholesky-decomposed electron repulsion integrals,” J. Chem. Phys. , in press (2022).
  * Burger  _et al._ [2021] S. Burger, F. Lipparini, J. Gauss, and S. Stopkowicz, “NMR chemical shift computations at second-order Møller–Plesset perturbation theory using gauge-including atomic orbitals and Cholesky-decomposed two-electron integrals,” J. Chem. Phys. 155, 074105 (2021).
  * Matthews  _et al._ [2020] D. A. Matthews, L. Cheng, M. E. Harding, F. Lipparini, S. Stopkowicz, T.-C. Jagau, P. G. Szalay, J. Gauss, and J. F. Stanton, “Coupled-cluster techniques for computational chemistry: The CFOUR program package,” J. Chem. Phys. 152, 214108 (2020).
  * [50] J. F. Stanton, J. Gauss, L. Cheng, M. E. Harding, D. A. Matthews, and P. G. Szalay, “CFOUR, Coupled-Cluster techniques for Computational Chemistry, a quantum-chemical program package,” With contributions from A. Asthana, A.A. Auer, R.J. Bartlett, U. Benedikt, C. Berger, D.E. Bernholdt, S. Blaschke, Y. J. Bomble, S. Burger, O. Christiansen, D. Datta, F. Engel, R. Faber, J. Greiner, M. Heckert, O. Heun, M. Hilgenberg, C. Huber, T.-C. Jagau, D. Jonsson, J. Jusélius, T. Kirsch, M.-P. Kitsaras, K. Klein, G.M. Kopper, W.J. Lauderdale, F. Lipparini, T. Metzroth, L.A. Mück, D.P. O’Neill, T. Nottoli, J. Oswald, D.R. Price, E. Prochnow, C. Puzzarini, K. Ruud, F. Schiffmann, W. Schwalbach, C. Simmons, S. Stopkowicz, A. Tajti, J. Vázquez, F. Wang, J.D. Watts, C. Zhang, X. Zheng, and the integral packages MOLECULE (J. Almlöf and P.R. Taylor), PROPS (P.R. Taylor), ABACUS (T. Helgaker, H.J. Aa. Jensen, P. Jørgensen, and J. Olsen), and ECP routines by A. V. Mitin and C. van Wüllen. For the current version, see http://www.cfour.de.
  * Lipparini and Gauss [2016] F. Lipparini and J. Gauss, “Cost-effective treatment of scalar relativistic effects for multireference systems: A CASSCF implementation based on the spin-free Dirac-Coulomb Hamiltonian,” J. Chem. Theory Comput. 12, 4284–4295 (2016).
  * Nottoli, Gauss, and Lipparini [2021] T. Nottoli, J. Gauss, and F. Lipparini, “Second-order CASSCF algorithm with the Cholesky decomposition of the two-electron integrals,” J. Chem. Theory Comput. 17, 6819–6831 (2021).
  * Helgaker and Jørgensen [1991] T. Helgaker and P. Jørgensen, “An electronic Hamiltonian for origin independent calculations of magnetic properties,” J. Chem. Phys. 95, 2595–2601 (1991).
  * Helgaker and Almlöf [1984] T. U. Helgaker and J. Almlöf, “A second-quantization approach to the analytical evaluation of response properties for perturbation-dependent basis sets,” Int. J. Quantum Chem. 26, 275–291 (1984).
  * Simons, Jørgensen, and Helgaker [1984] J. Simons, P. Jørgensen, and T. U. Helgaker, “Higher molecular-deformation derivatives of the configuration-interaction energy,” Chem. Phys. 86, 413–432 (1984).
  * Olsen  _et al._ [1995] J. Olsen, K. L. Bak, K. Ruud, T. Helgaker, and P. Jørgensen, “Orbital connections for perturbation-dependent basis sets,” Theor. Chimica Acta 90, 421–439 (1995).
  * Ruud  _et al._ [1995b] K. Ruud, T. Helgaker, J. Olsen, P. Jørgensen, and K. L. Bak, “A numerically stable orbital connection for the calculation of analytical Hessians using perturbation-dependent basis sets,” Chem. Phys. Lett. 235, 47–52 (1995b).
  * Jensen and Ågren [1986] H. J. Aa. Jensen and H. Ågren, “A direct, restricted-step, second-order MC SCF program for large scale ab initio calculations,” Chem. Phys. 104, 229–250 (1986).
  * Folkestad, Kjønstad, and Koch [2019] S. D. Folkestad, E. F. Kjønstad, and H. Koch, “An efficient algorithm for Cholesky decomposition of electron repulsion integrals,” J. Chem. Phys. 150, 194112 (2019).
  * Zhang  _et al._ [2021] T. Zhang, X. Liu, E. F. Valeev, and X. Li, “Toward the minimal floating operation count Cholesky decomposition of electron repulsion integrals,” J. Phys. Chem. A 125, 4258–4265 (2021).
  * Gauss  _et al._ [2021] J. Gauss, F. Lipparini, S. Burger, S. Blaschke, M.-P. Kitsaras, and S. Stopkowicz, (2021), Johannes Gutenberg-Universität Mainz, unpublished.
  * Gauss  _et al._ [2022] J. Gauss, S. Blaschke, S. Burger, T. Nottoli, F. Lipparini, and S. Stopkowicz, “Cholesky decomposition of two-electron integrals in quantum-chemical calculations with perturbative or finite magnetic fields using gauge-including atomic orbitals,” (2022), arXiv:2205.07035.
  * Pulay and Hamilton [1988] P. Pulay and T. P. Hamilton, “UHF natural orbitals for defining and starting MC-SCF calculations,” J. Chem. Phys. 88, 4926–4933 (1988).
  * Tóth and Pulay [2020] Z. Tóth and P. Pulay, “Comparison of methods for active orbital selection in multiconfigurational calculations,” J. Chem. Theory Comput. 16, 7328–7341 (2020).
  * Dunning, Jr. [1989] T. H. Dunning, Jr., “Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen,” J. Chem. Phys. 1007, 4572–4585 (1989).
  * Hehre, Ditchfield, and Pople [1972] W. J. Hehre, R. Ditchfield, and J. A. Pople, “Self-consistent molecular orbital methods. XII . Further extensions of Gaussian-type basis sets for use in molecular orbital studies of organic molecules,” J. Chem. Phys. 56, 2257–2261 (1972).
  * Becke [1993] A. D. Becke, “A new mixing of Hartree-Fock and local density-functional theories,” J. Chem. Phys. 98, 1372–1377 (1993).
  * Frisch  _et al._ [2016] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, “Gaussian 16 Revision A.03,” (2016), Gaussian Inc., Wallingford CT.
  * Jensen and Jørgensen [1984] H. J. Aa. Jensen and P. Jørgensen, “A direct approach to second-order MCSCF calculations using a norm extended optimization scheme,” J. Chem. Phys. 80, 1204–1214 (1984).
  * Sayfutyarova  _et al._ [2017] E. R. Sayfutyarova, Q. Sun, G. K. L. Chan, and G. Knizia, “Automated construction of molecular active spaces from atomic valence orbitals,” J. Chem. Theory Comput. 13, 4063–4078 (2017).
  * Stein and Reiher [2016] C. J. Stein and M. Reiher, “Automated selection of active orbital spaces,” J. Chem. Theory Comput. 12, 1760–1771 (2016).
  * Sayfutyarova and Hammes-Schiffer [2019] E. R. Sayfutyarova and S. Hammes-Schiffer, “Constructing molecular $\pi$-orbital active spaces for multireference calculations of conjugated systems,” J. Chem. Theory Comput. 15, 1679–1689 (2019).