Modern quantum chemistry with [Open]Molcas

Francesco Aquilante, Jochen Autschbach, Alberto Baiardi, Stefano Battaglia, Veniamin Borin, Liviu Chibotaru, Irene Conti, Luca de Vico, Mickaël Delcey, Ignacio Fdez. Galván, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Francesco Aquilante, Jochen Autschbach, Alberto Baiardi, Stefano Battaglia, Veniamin Borin, et al. Modern quantum chemistry with [Open]Molcas. The Journal of Chemical Physics, 2020, 152 (21), pp.214117. 10.1063/5.0004835. hal-02814121

HAL Id: hal-02814121

https://amu.hal.science/hal-02814121

Submitted on 6 Jun 2020

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Modern Quantum Chemistry with [Open]Molcas

Francesco Aquilante,¹, a) Jochen Autschbach,², b) Alberto Baiardi,³, c) Stefano Battaglia,⁴, d) Veniamin A. Borin,⁵, e) Liviu F. Chibotaru,⁶, f) Irene Conti,⁷, g) Luca De Vico,⁸, h) Mickaël Delcey,⁹, i) Ignacio Fdez. Galván,⁴, j) Nicolas Ferré,¹⁰, k) Leon Freitag,³, l) Marco Garavelli,⁷, m) Xuejun Gong,¹¹, n) Stefan Knecht,³, o) Ernst D. Larsson,¹², p) Roland Lindh,⁴, q) Marcus Lundberg,⁹, r) Per Åke Malmqvist,¹², s) Artur Nenov,¹³, t) Jesper Norell,¹⁴, u) Michael Odelius,¹⁴, v) Massimo Olivucci,⁸, 15, w) Thomas B. Pedersen,¹⁶, x) Laura Pedraza-González,⁸, y) Quan M. Phung,¹⁷, z) Kristine Pierloot,⁶, aa) Markus Reiher,³, bb) Igor Schapiro,⁵, cc) Javier Segarra-Martí,¹⁸, dd) Francesco Segatta,⁷, ee) Luis Seijo,¹⁹, ff) Saumik Sen,⁵, gg) Dumitru-Claudiu Sergentu,², hh) Christopher J. Stein,³, ii) Liviu Ungur,¹¹, jj) Morgane Vacher,²⁰, kk) Alessio Valentini,²¹, ll) and Valera Veryazov¹², mm)

¹) Theory and Simulation of Materials (THEOS) and National Centre for Computational Design and Discovery of Novel Materials (MARVEL), École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

²) Department of Chemistry, University at Buffalo, Buffalo, NY 14260-3000, USA

³) ETH Zurich, Lab. f. Phys. Chem., Vladimir-Prelog-Weg 2, 8093 Zurich, Switzerland

⁴) Department of Chemistry - BMC, Uppsala University, P.O. Box 576, SE-751 23 Uppsala, Sweden

⁵) Fritz Haber Center for Molecular Dynamics Research, Institute of Chemistry, The Hebrew University of Jerusalem, Jerusalem 9190401, Israel

⁶) Department of Chemistry, KU Leuven, Celestijnenlaan 200F, 3001 Leuven, Belgium

⁷) Dipartimento di Chimica Industriale “Toso Montanari”, Università di Bologna, Viale del Risorgimento 4, Bologna, I-40136 Italy

⁸) Dipartimento di Biotecnologie, Chimica e Farmacia, Università degli Studi di Siena, via Aldo Moro 2, 53100, Siena, Italy

⁹) Department of Chemistry - Ångström Laboratory, Uppsala University, SE-751 21 Uppsala, Sweden

¹⁰) Aix-Marseille Univ, CNRS, Institut Chimie Radiculaire, Marseille, France

¹¹) Department of Chemistry, University of Singapore, 3 Science Drive 3, 117543 Singapore

¹²) Division of Theoretical Chemistry, Lund University, POB 124, Lund 22100, Sweden

¹³) Dipartimento di Chimica Industriale Toso Montanari, Università di Bologna, Viale del Risorgimento 4, 40136 Bologna, Italy

¹⁴) Department of Physics, AlbaNova University Center, Stockholm University, SE-106 91 Stockholm, Sweden

¹⁵) Department of Chemistry, Bowling Green State University, Bowling Green, Ohio 43403, United States of America

¹⁶) Hylleraas Centre for Quantum Molecular Sciences, Department of Chemistry, University of Oslo, P.O. Box 1033 Blindern, N-0315 Oslo, Norway

¹⁷) Institute of Transformative Bio-Molecules (WPI-ITbM), Nagoya University, Chikusa, Nagoya 464-8602, Japan

¹⁸) Department of Chemistry, Molecular Sciences Research Hub, Imperial College London, White City Campus, 80 Wood Lane, London W12 0BZ, UK

¹⁹) Departamento de Química, Instituto Universitario de Ciencia de Materiales Nicolás Cabrera, and Condensed Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid, 28049 Madrid, Spain

²⁰) Laboratoire CEISAM - UMR CNRS 6230, Université de Nantes, 44300 Nantes, France

²¹) Theoretical Physical Chemistry, Research Unit MolSys, Université de Liège, Allée du 6 Août, 11, 4000 Liège, Belgium

(Dated: 11 May 2020)

[Open]Molcas is an ab initio electronic structure program providing large set of computational methods from Hartree-Fock and Density Functional Theory to various implementations of Multiconfigurational Theory. This article provides a comprehensive overview of the main features of the code, specifically reviewing the use of the code in previously reported chemical applications as well as more recent applications including the calculation of magnetic properties from optimized density matrix renormalization group wave functions.

<footnote>a) Electronic mail: francesco.aquilante@epfl.ch</footnote>

<footnote>b) Electronic mail: jochena@buffalo.edu</footnote>
**Contributions:** All authors contributed equally to this work.

**Author to whom correspondence should be addressed:** valera.veryazov@teokem.lu.se

# I. INTRODUCTION

Modern quantum chemistry is impossible without a versatile computational software, which includes calculation of integrals, optimization of wave functions and properties computation. Not surprisingly, computational codes in the field are grouping into large packages, which simplifies the development process and their usage.

The list of quantum chemistry computer programs, maintained at Wikipedia¹, contains almost one hundred different computational codes with a large overlap in the functionality. Even in a narrow field of the codes, computing the electronic structure of the ground and excited states of molecules, a researcher can choose between a

large set of either well-established or newly developed codes. In addition to functionality from the scientific point of view, the codes are very different with respect to performance, the need of hardware resources, documentation and user friendliness.

As a general purpose package for quantum chemical calculations, MOLCAS, has made a long journey from a collection of homemade codes to a professional software with distributed development, automatic verification, user support, etc. Recently the vast majority of the codes in MOLCAS have been released as open source - the OpenMolcas project. This package is user-friendly and ready-to-use, and is also a developers' platform. While OpenMolcas is a free-of-charge package, with web based community support, the MOLCAS package is a licensed distribution of a refurbished version of the OpenMolcas package, with some additional external utilities. The MOLCAS package is user-oriented and offers explicit support in installation, performance enhancement and basic calculations. For all practical purposes both packages offer the same or similar computational tools. In what follows we will make the distinction between MOLCAS and OpenMolcas when so is called for, else-wise we will refer to both distributions as [Open]Molcas.

There are several reviews describing new and recent features of MOLCAS²⁻⁴ and OpenMolcas⁵, as well as the development infrastructure of the code⁶, external software interfaced to MOLCAS and Graphical User Interface programs⁷,⁸. For a more comprehensive list of features, computed properties, implementation details and theoretical description, we would suggest the reader to consult these papers, together with the manual⁹.

The purpose of the current paper is to highlight not necessarily the newest, but the most commonly used features of [Open]Molcas with a focus to the general outcome for computational chemistry. In order to reduce the size of the paper, the computational details (including input examples) are presented in the Supplementary Information. Thus, we encourage the reader interested in a deeper understanding of the possibilities of the [Open]Molcas package, to consult the SI section of this paper.

# II. PROGRAM OVERVIEW

From the programming point of view, [Open]Molcas is a package that is composed of a large number of individual computational units, which run together driven by an input parser. The modularity simplifies the dis-

<footnote>c) Electronic mail: alberto.baiardi@phys.chem.ethz.ch</footnote>

<footnote>d) Electronic mail: stefano.battaglia@kemi.uu.se</footnote>

<footnote>e) Electronic mail: veniamin.borin@mail.huji.ac.il</footnote>

<footnote>f) Electronic mail: Liviu.Chibotaru@kuleuven.be</footnote>

<footnote>g) Electronic mail: irene.conti@unibo.it</footnote>

<footnote>h) Electronic mail: Luca.DeVico@unisi.it</footnote>

<footnote>i) Electronic mail: mickael.delcey@kemi.uu.se</footnote>

<footnote>j) Electronic mail: Ignacio.Fernandez@kemi.uu.se</footnote>

<footnote>k) Electronic mail: nicolas.ferre@univ-amu.fr</footnote>

<footnote>l) Electronic mail: leon.freitag@phys.chem.ethz.ch</footnote>

<footnote>m) Electronic mail: marco.garavelli@unibo.it</footnote>

<footnote>n) Electronic mail: xuejun.gong@nus.edu.sg</footnote>

<footnote>o) Electronic mail: stefan.knecht@phys.chem.ethz.ch</footnote>

<footnote>p) Electronic mail: ernst_dennis.larsson@teokem.lu.se</footnote>

<footnote>q) Electronic mail: roland.lindh@kemi.uu.se</footnote>

<footnote>r) Electronic mail: marcus.lundberg@kemi.uu.se</footnote>

<footnote>s) Electronic mail: Per-Ake.Malmqvist@teokem.lu.se</footnote>

<footnote>t) Electronic mail: artur.nenov@unibo.it</footnote>

<footnote>u) Electronic mail: jesper.norell@fysik.su.se</footnote>

<footnote>v) Electronic mail: odelius@fysik.su.se</footnote>

<footnote>w) Electronic mail: olivucci@unisi.it</footnote>

<footnote>x) Electronic mail: t.b.pedersen@kjemi.uio.no</footnote>

<footnote>y) Electronic mail: la.pedrazagonzalez@student.unisi.it</footnote>

<footnote>z) Electronic mail: quan.phung@itbm.nagoya-u.ac.jp</footnote>

<footnote>aa) Electronic mail: kristin.pierloot@kuleuven.be</footnote>

<footnote>bb) Electronic mail: markus.reiher@phys.chem.ethz.ch</footnote>

<footnote>cc) Electronic mail: igor.schapiro@mail.huji.ac.il</footnote>

<footnote>dd) Electronic mail: j.segarra-marti@imperial.ac.uk</footnote>

<footnote>ee) Electronic mail: francesco.segatta2@unibo.it</footnote>

<footnote>ff) Electronic mail: luis.seijo@uam.es</footnote>

<footnote>gg) Electronic mail: saumik.sen@mail.huji.ac.il</footnote>

<footnote>hh) Electronic mail: dumitruc@buffalo.edu</footnote>

<footnote>ii) Electronic mail: christopher.stein@phys.chem.ethz.ch</footnote>

<footnote>jj) Electronic mail: chmlu@nus.edu.sg</footnote>

<footnote>kk) Electronic mail: morgane.vacher@univ-nantes.fr</footnote>

<footnote>ll) Electronic mail: alessio.valentini@uliege.be</footnote>

<footnote>mm) Electronic mail: valera.veryazov@teokem.lu.se</footnote>
tributed development and allows to have a better con-
trol over computational resources, especially memory. As
a drawback of the module structure of [Open]Molcas
there is a need to organize communication between exe-
cutable codes, which includes handling of return codes,
transferring the data, and keeping the parallel infrastruc-
ture. Currently there are two implementations of input
parser, written in Perl/C (molcas.exe in MOLCAS) and
in Python (pymolcas in OpenMolcas). Both parsers use
Enhanced MOLCAS Input Language (EMIL), described
in details in the SI.

The main computer language used in [Open]Molcas
is Fortran. The majority of the code is compatible
with the Fortran 77 standard, but some newer contri-
butions are written in Fortran 90 or later. All system-
dependent routines, like handling of memory, I/O and
parallelization, are written in C language. Not surpris-
ingly, the main computational tasks are related to lin-
ear algebra problems. The code has an interface to stan-
dard BLAS/LAPACK libraries, providing the best per-
formance for production calculations, but at the same
time there is a possibility to control the usage of mathe-
matical libraries for development purposes.

Despite the general trend of reducing the variety of
hardware architecture and software, we continue to sup-
port [Open]Molcas not only on the mainstream setup
(x86.64 CPU with Linux OS and GNU or Intel compiler),
but also on marginal hardware and compilers. The com-
plete list of supported platforms, compilers and libraries
is provided in the SI. This policy restricts the usage of
some, usually nonstandard or recently standardized, fea-
tures of compilers, but at the same time it helps to reduce
the number of bugs in the code.

Since MOLCAS 6.0, a robust verification system is
used to maintain the stability of the code. An application
developer includes a call to a library function, specifying
as parameters: a label, a value and an allowed thresh-
old for any computed data. The verification procedure
can use this call either to store the computed data, or
to verify it against the reference. The verification can
be performed either locally, or at remote installations of
[Open]Molcas (with the use of a special selection of com-
pilers, compilation flags of libraries).

A breakthrough in performance originated from
the implementation of what is known as “Cholesky
infrastructure”.⁸,¹⁰,¹¹ In practice, every quantum chem-
istry method implemented in OpenMolcas is formulated
directly in terms of the so-called Cholesky vectors instead
of the traditional formulation based on two-electron in-
tegrals. Next to the massive savings in storage, reduced
scaling is achieved in the evaluation of most tensor quan-
tities used within the algorithms for mean-field and cor-
relation methods energy and properties methods.

[Open]Molcas computes the optimized wavefunction
and the corresponding energies for a large set of approxi-
mate Hamiltonians, including Hartree-Fock (HF), Kohn-
Sham Density Functional Theory (DFT), Multiconfigu-
rational methods (Complete Active Space (CASSCF),

Restricted Active Space (RASSCF), Generalized Ac-
tive Space (GASSCF), Density Matrix Renormaliza-
tion Group (DMRG)), Perturbation theory of the sec-
ond order (MP2, CASPT2, RASPT2), Valence Bond
theory (CASVB), Coupled Cluster theory (CCSD,
CCSD(T)), Multireference Configuration Interaction (M
RCI). [Open]Molcas is capable to optimize the geome-
try of a molecule (using either analytical or numerical
gradients). Most of the methods allow to compute the
electronic structure of a molecule in a medium (solution,
solids). Recent development of [Open]Molcas is also fo-
cused on computing various properties.

It is difficult to present strict and accurate limitations
of the codes in [Open]Molcas, but on current hardware,
it is possible to study a molecule with more than a hun-
dred atoms, with about two thousand basis functions.
The limits in terms of the active space selection depend
on the method and the hardware configuration. A user
should be aware that multiconfigurational methods re-
quire a very large amount of memory: 2 GB for a small
calculation and 16–32 GB for a large, but not extreme
calculation. Most of the modules are parallelized, how-
ever, a production run will require even more memory if
parallelization is used¹².

Recent development is also focused on creating data in-
terface with other computational and visualization codes.
In some cases the data exchange is customized (interface
to other computational codes, e.g. QCMAQUIS, Colum-
bus, etc.), but the work in progress also includes a gen-
eral HDF5-formatted interface and XML/json formatted
output.

The following sections present the main computational
codes in [Open]Molcas. The flowchart of modules, as well
as a short programming description of codes is presented
in the SI. In order to distinguish, say, the SCF method
and the SCF program, the latter is printed in small capital
letters.

A. Integrals, gradients, and second derivatives

The original MOLCAS package has since its in-
ception utilized Gaussian basis sets to expand the
one-particle space. While the original MOLCAS 1.0
used MOLECULE¹³, subsequent versions have used the
SEWARD,¹⁴ ALASKA,¹⁵ and MCKINLEY¹⁶ programs to
compute integrals and up to second order derivatives
of those. The hallmark of these modules is the efficient
computation of integrals in the Gaussian basis of gen-
erally contracted basis functions. Moreover, the integral
codes are implemented for the efficient use of real spher-
ical harmonics of arbitrary order. The two-electron inte-
grals, and associated derivatives are computed using the
Rys–Gauss quadrature, while one-electron integrals are
computed using the Hermite–Gauss quadrature. As any
standard package OpenMolcas supports a large array of
one-electron operators, but also some which are unique.
Here it is worth to mention integrals over the exact
semi-relativistic operator of the interaction between light and matter.¹⁷ The two-electron integrals are also implemented in an integral-direct fashion in the scf module.¹⁸ Recently the conventional integral generation to disc has been less used as this has been replaced with an efficient implementation of Cholesky decomposition and so-called density-fitting techniques, a method which reduces computational times one order of magnitude or more. This matter will to some extent be discussed in the following subsection.

## 1. Cholesky Infrastructure

Most quantum chemistry software need to cope with the computational bottlenecks arising from the evaluation and storage of the two-electron integrals. A strategy that has become popular over the years with the names of “Density Fitting” (DF) or “Resolution of the Identity” (RI) uses the following approximate tensor decomposition in terms of 3-center and 2-center integrals:

$$ (\mu\nu|\lambda\sigma) \approx \sum_{PQ} (\mu\nu|P)[(P|Q)]^{-1}(Q|\lambda\sigma), \quad (1) $$

where sets of so-called auxiliary basis functions P, Q are designed for each atomic basis set through data-fitting of specific energy contributions—e.g., second-order Møller-Plesset (MP2) correlation energy—and then used to calculate the two-electron integrals $(\mu\nu|\lambda\sigma)$ over atomic orbitals. Despite the need to compute the matrix inversion over the auxiliary functions $[(P|Q)]^{-1}$, significant computational advantages come from employing Eq. (1) due to the fact that the number of auxiliary functions is only about 3 to 5 times the total size of the atomic basis set, thus potentially orders of magnitude smaller than the leading dimension of the full integral matrix. Therefore, the number of 2- and 3-center integrals is much smaller than the number of otherwise needed two-electron (and up to 4-center) integrals.

On the other hand, it is now a well-established fact that such generation of the auxiliary basis set can be made without data-fitting or bias towards a specific quantum chemistry method, if one exploits the onset of numerical linear dependence in the product basis of the atomic orbitals.¹⁹ This type of ab initio DF originated in the earlier idea of employing the scalar product of vectors from the Cholesky decomposition (CD) of the integral matrix in order to approximate the integrals:

$$ (\mu\nu|\lambda\sigma) \approx \vec{\mathcal{L}}_{\mu\nu} \cdot \vec{\mathcal{L}}_{\lambda\sigma}, \quad (2) $$

with computational advantages arising from the fact that the length of each CD vector is only about the same as the size of a standard auxiliary basis set used in Eq. (1). The main appeal of standard CD over RI approximations is that one can, at the price of increasing the CD vector length, effectively control the accuracy of the integral representation by means of the threshold used for the

incomplete matrix decomposition. This is a very important property especially in the context of highly accurate methods, as it helps preserving the systematic improvement of such quantum chemical models. The use of ab initio DF combines the accuracy of the CD representation with the computational ease of Eq. (1) by requiring only CD of each atomic subblock (acCDs) of the integral matrix to define the auxiliary basis.²⁰ Once generated on the fly, such an auxiliary basis set is used in Eq. (1) and in this way one can also compute analytical derivatives for the two-electron integrals, as needed for potential energy surface (PES) exploration at any level of theory.²¹,²² The option to perform acCD based calculations is available in [Open]Molcas by means of the keyword RICD in GATEWAY and SEWARD. Additional features of the Cholesky infrastructure deal with exploiting a fully-localized reformulation of Eq. (1) for use in connection with linear-scaling correlation methods, especially non-variational ones, as any local DF two-electron integral approximation may cause variational collapse.²³,²⁴

## B. The RASSCF program

In a “Multi-Reference Self-Consistent Field” (MCSCF) calculation, also called “Multi-Configuration Hartree-Fock” (MCHF), the electronic wave function is composed of a number of different configurations, with different occupation numbers. Orbitals that are doubly occupied in all the configurations are variously called “core” or “inactive”, while those with varying occupation are called “active”. Moreover, the configuration functions, as well as the orbital functions, are optimized so as to make the energy of the MCSCF state optimal, i.e. minimized or stationary.

General MCSCF optimization is usually very hard to converge to a solution, and is also marred by multiple local minima. In contrast, the CASSCF, i.e. “Complete” active space, procedure²⁵ is usually well-behaved numerically. The term “complete” implies that all configurations that are formed by distributing electrons in all possible ways among the active orbitals are used in the Configuration Interaction (CI). Use of a fast and robust CI solver and selecting the energy eigenfunction then makes the energy a function of only the “active space”, rather than of the individual orbitals. This makes the optimization procedure simpler and, in most cases, rather robust.

However, the number of configurations, i.e. the size of the CI eigenvalue problem, will grow in a very inconvenient way beyond a dozen or so active orbitals. One way out is to restrict the active space a bit, by disallowing some of the configurations. By requiring certain orbital subspaces to have a restricted range of occupation numbers, we get a drastic reduction in the number of configurations, but also a more difficult optimization problem. This kind of restriction is called “Generalized” active space, i.e. the acronym is GASSCF²⁶,²⁷. A less problematic, but also less powerful, method is the “Re-
stricted” active space, RASSCF<sup>28</sup>. Here we can concep- tually start with a big CAS, but then restrict it by re- quiring the orbitals in one subspace, the RAS1 space, to be fully occupied with just a few exceptions, i.e. the number of “holes” in RAS1 is limited. Similarly, there is a RAS3 space, which is empty except for a maximum al- lowed number of electrons. The remaining active orbitals belong to the RAS2 space.

The RASSCF program was originally computing only CASSCF wave functions, but was later expanded by al- lowing RASSCF, and later also GASSCF. However, it has also become possible to do similar calculations with a huge number of active orbitals, when the optimization is performed by relaxation procedures called “Density Matrix Renormalization Group” methods, which grew out of methods used in lattice quantum dynamics. Some DMRG codes have now been interfaced as external mod- ules, which can be optionally linked to OpenMolcas.

Multireference calculations can in principle be ap- plied to arbitrary kinds of electronic systems, regard- less of charge, spin and point group symmetry and with molecules with any conformation and level of excitation. However, they can require large resources in terms of ba- sis sets, many active orbitals, and large CI expansions. The calculations outlined above will usually require also a calculation of residual dynamic correlation.

CASSCF cannot be used to optimize more than 12–16 active orbitals. Many more are possible using RASSCF or GASSCF. Even then, the optimization procedure yields a wave function that is optimized with the re- striction that the virtual orbitals have zero occupation in the wave function. Such a wave function includes non- dynamic correlation, but is lacking much dynamic corre- lation, which can also differ considerably between differ- ent states/geometries.

Thus, the CASSCF (as well as RASSCF and GASSCF) calculation will usually be followed by a perturbative cal- culation of the missing dynamic correlation. In fact, for CASSCF the accuracy is usually comparable to that of a Hartree–Fock (HF), when the state in case is well suited for such calculations. The difference is that CASSCF is able to give results of uniform quality, also for radicals, broken bonds, excited states, etc. RASSCF and GASSCF are much better in allowing at least some non-dynamic correlation, but often not enough. One procedure to go further will be briefly described in what follows.

C. The CASPT2 program

The CASPT2 program solves the Raleigh–Schrödinger perturbation equation

$$ (\hat{H}_0 - E_0)\Psi_1 = -(\hat{H}_1 - E_1)\Psi_0 \quad (3) $$

where $Ψ_0$, the zeroth-order wave function, is a CASSCF, or more generally, a RASSCF wave function. $Ψ_1$, the first- order perturbation to the wave function, is expressed as a

large number—typically 100 000–1 000 000—of excitation amplitudes. The zeroth-order Hamiltonian is

$$ \hat{H}_0 = \hat{P}_0 \hat{F} \hat{P}_0 + \hat{P}_1 \hat{F} \hat{P}_1 \quad (4) $$

Here, $\hat{P}_0$ and $\hat{P}_1$ are projectors onto the reference func- tion and the interacting space, respectively, and $\hat{F}$ is an effective one-electron spin-summed operator. Its action on a multiconfigurational wave function is not diago- nal, and the equation is solved by an iterative method Preconditioned Conjugate Gradient (PCG). The inter- acting space is usually huge. It is subdivided into eight symmetry-blocked parts, further partitioned into symme- try block by index permutations and point group sym- metry, and each such block represents the action of some linear combinations of two-electron operators acting on the reference CASSCF or RASSCF state.

The linear combinations, and the blocking, are used to transform the perturbation equation into a form where it can be solved by efficient PCG iterations. Thus, the diag- onal blocks become themselves diagonal matrices, and so a very efficient solver is obtained for the pre-conditioning, which also gives a starting diagonal solution. This is then refined by typically 6–20 PCG iterations.

The result is a first-order perturbation wave function $Ψ_1$. This is then used to express the final energy in the form of a Hylleraas functional, i.e. the energy has an error that scales with the square of errors in $Ψ_1$.

Developments include e.g. multi-state CASPT2 (MS- CASPT2), extended MS-CASPT2 (XMS-CASPT2) — which uses a state-average $\hat{F}$ and ensures invariance under unitary rotations of the reference states — and extended dynamically weighted CASPT2 (XDW- CASPT2)<sup>29</sup>, a method that interpolates between the for- mer two and retains advantages from both.

The CASPT2 program can accept CASSCF or RASSCF reference wave functions, DMRG wave functions, if an op- tional DMRG package has been linked with OpenMolcas.

D. The RASSI program

For multi-configurational methods, a common problem is to compute overlaps and matrix elements of various op- erators over a set of wave functions for different states. The orbitals are usually optimized for the various states. If these are not too dissimilar and have the same spin and the same point group representation they may be treated together, giving an ensemble optimization (State Aver- aged CASSCF). If any interstate properties are sought for, we have the following problem: How to compute ma- trix elements over CI wave functions that use different orbital bases?

This was the original reason for developing the State Interaction programs. The original one, CASSI, took pair- wise states from a set of states, and could compute their overlaps and Hamiltonian matrix element, as well as ma- trix elements of property integrals such as e.g. multi-
pole moments. Using the wave functions as basis func-
tions, it could then compute orthonormal states and non-
interacting linear combinations of these, and the energy
eigenvalues and proper dipole transition probabilities, for
example³⁰.

This was shown to be efficient and practical, pro-
vided that the state functions had been computed using
a common set of AO basis functions, and for CASSCF
wave functions. It was primarily used for the purposes
indicated above, but also e.g. for transforming to dia-
batic wave function components for crossing or nearly
crossing states. Afterwards, of course a large number of
other disparate uses were found. RASSCF wave func-
tions can be used in a similar way³¹. States that dif-
fer in spin and/or number of electrons may be treated,
leading to so-called Spin-Orbit RASSI³² (SO-RASSI)
where the spin-free RASSI wave functions are combined
with spin functions, and using so-called “Atomic Mean-
Field Integrals” (AMFI) spin-orbital Hamiltonian inte-
grals from SEWARD, properties end energies involving in-
dividual spin components can be computed. Also, using
different number of electrons, so-called Dyson amplitudes
can be obtained, yielding e.g. probabilities for core-hole
dynamics. RASSI can also compute so-called “binatural
orbitals”³³, which describe pairs of states. The imple-
mentation is new, and has not yet been much used.

In conclusion, the RASSI program is a fairly versatile
tool for handling multiconfigurational wave functions in
many different ways. In the future, it should also have
the capability to combine states having different basis
functions, or different number of inactive/active orbitals.

E. External programs for DMRG-based calculations

For a general description of DMRG and the inter-
faces available for performing these calculations with
[Open]Molcas, the reader is advised to refer to refer-
ence 5.

1. DMRG-based Electronic Structure with QCMAQUIS

QCMAQUIS³⁴⁻³⁶ is a stand-alone software module that provides algorithms utilizing the DMRG approach³⁷,³⁸ in quantum chemistry³⁹. Targeting large active orbital spaces beyond restrictions of traditional CAS approaches in [Open]Molcas calls for either specialized approaches such as the ones discussed in Section II F or new wave function parametrizations which allow one to circumvent the exponential growth of basis states with system size in a CASCI formulation. The optimization of matrix-product state (MPS) ground- and excited-state wavefunctions by the module QCMAQUIS, which stands for Quantum Chemical Modern Algorithms for Quantum Interacting Systems, exploits a genuine matrix-product operator formalism of the (relativistic) quantum chemical Hamiltonian³⁴,³⁵,⁴⁰.

By means of a recently established Fortran-to-C++
interface, QCMAQUIS is modularly integrated into the
current CASSCF solver of OpenMolcas (under the
alias DMRGSCF). Hence, DMRG-SCF calculations with
QCMAQUIS<sup>41</sup> allow for—in analogy to the traditional
CASSCF approach in OpenMolcas—(i) the inclusion of
equilibrium and non-equilibrium solvent models<sup>42</sup>, (ii)
the computation of state-specific<sup>43</sup> as well as state-
averaged<sup>44</sup> ground- and excited gradients, (iii) the ap-
proximation of dynamical electron correlation in post-
DMRG-SCF calculations second-order multi-reference
perturbation theory<sup>45</sup>, (iv) also multi-configurational
pair-density functional theory<sup>46</sup> can be employed for this
purpose), (v) moreover, the capabilities of the RASSI
approach (see Section IID) to handle non-orthogonal
CI-type wavefunctions as input states are matched with
the MPS state-interaction approach (MPSSI)<sup>47</sup> of QC-
MAQUIS. The latter is accessed in OpenMolcas with the
MPSSI program. This capability is decisive for the calcula-
tion of transition dipole moments as well as SO coupling
matrix elements between states optimized separately as
MPS-type wave functions. In Ref. 47 we showed with
the example of actinide complexes that DMRG-SCF in
combination with the MPSSI approach yields magnetic
properties such as ground- and excited-state *g*-tensors
in full agreement with data obtained from correspond-
ing CASSCF/RASSI calculations. For an example that
combines for the first time the DMRG-SCF/PT2/MPSSI
functionality with the SINGLE_ANISO approach, the
reader is referred to Section III B 6. In addition, by tak-
ing advantage of the same implementation framework as
the RASSI module, the MPSSI module facilitates the com-
putation of Spin Orbit (SO) -coupled wavefunctions that
are suitable for a subsequent ground- and excited-state
bonding analysis as described in Section IID 4 and/or
that can serve as input functions to the SINGLE_ANISO
module<sup>4</sup>.

The QCMAQUIS–OpenMolcas interface can target ex-
cited states within either a state-specific (SS)³⁴ or a state-
average (SA) formulation of DMRG. The SA formulation
is realized based on the multi-state matrix product state
idea (MS-MPS)⁴⁴,⁴⁸ that generalizes the MPS concept to
multiple electronic states. MS-MPS relies on a common
tensor to encode all target electronic states that there-
fore can be optimized simultaneously with a lower com-
putational cost and a higher stability in the presence of
nearly-degenerate electronic states, e.g. in the vicinity of
a conical intersection. We presented an extension of the
Lagrangian-based SA-CASSCF linear response theory⁴⁹
to MS-MPS wave functions which allows us to calculate
SA-DMRG-SCF energy gradients and non-adiabatic cou-
plings with QCMAQUIS. In future work, we plan to fur-
ther improve the efficiency of the SA-DMRG-SCF gra-
dients calculation by exploiting the structure of an MS-
MPS to efficiently evaluate the contractions required for
the calculations of one- and two-body reduced density
matrices for a given MS-MPS.

As indicated above, the QCMAQUIS–OpenMolcas in-
terface contains an implementation of the second-
order *n*-electron valence state perturbation theory
(NEVPT2)⁵⁰⁻⁵², which employs an MPS reference wave-
function and exploits the Cholesky decomposition of
two-electron integrals (CD-DMRG-NEVPT2)⁴⁵. Quasi-
degenerate NEVPT2 (QD-NEVPT2)⁵³, i.e., a genuine
multi-state formulation of NEVPT2, is also supported.
Prominent examples of a successful application of
DMRG-NEVPT2 include the calculation of spin-state
energetics⁵,⁴⁵ as well as dissociation energies⁵⁴ of sev-
eral large transition metal complexes. The main compu-
tational bottleneck present in NEVPT2 calculations, es-
pecially with large active orbital spaces, is the evalua-
tion of the four-particle reduced density matrix (4-RDM),
which scales as the eighth power with the number of ac-
tive orbitals *N*. Although DMRG supports active spaces
of up to about *N* = 100 orbitals, the prohibitive scal-
ing of the 4-RDM limits the current implementation in
OpenMolcas to *N* ≤ 22 orbitals. An interface for a mas-
sively parallel computation of 4-RDM matrix elements
on a cluster is provided, and in a future release efficiency
improvements in 4-RDM element computation will be im-
plemented, raising this limit to *N* ≤ 30 active orbitals.

Although the present implementation of QCMAQUIS
exploits an efficient, OpenMP-based shared-memory par-
allelization scheme, MPS optimizations for very large
active orbital spaces (N > 60) can suffer from severe,
memory-bound limitations for a given number of renor-
malized block states, the bond dimension *m*. To ad-
dress this issue adequately, future releases of QCMAQUIS
will benefit from a state-of-the-art hybrid OpenMP-MPI
parallelization framework that takes advantage of the
shared-memory functionalities across multiple compute
nodes provided by the MPI-3 standard. Developments
along these lines are currently ongoing.

2. DMRG with [Open]Molcas-CheMPS2 interface

[Open]Molcas-CheMPS2 is a simple, yet robust DMRG interface between [Open]Molcas and the CheMPS2 library,⁵⁵ allowing one to perform very efficient DMRG-CASSCF and (especially) DMRG-CASPT2 calculations. To the best of our knowledge, DMRG-CASPT2 is only available in three codes: [Open]Molcas, CheMPS2,⁵⁶ and orz.⁵⁷ As compared to the latter two, the uniqueness and strength of the implementation in OpenMolcas is that the Cholesky decomposition technique is fully supported. Thus, calculations of large molecules with several thousand basis functions are feasible.⁵⁸ DMRG-CASPT2 in OpenMolcas has been employed to tackle difficult problems involving mono-⁵⁹⁻⁶¹ and di-nuclear transition metal complexes.⁵⁸,⁶²⁻⁶⁴ In [Open]Molcas, the DMRG wave function and the RDMs required for CASPT2 calculations (2-, 3-RDMs, and the generalized Fock matrix contracted with the 4-RDM (F.4-RDM)) are calculated using the CheMPS2 library. In the current implementation, the computa-

tional cost of F.4-RDM is significantly reduced by em-
ploying the pseudocanonical orbital basis. This strategy,
however, is more suitable for highly symmetric or com-
pact molecules and limits the active space to ~30 or-
bitals. For general molecules and calculations requiring
a larger number of active orbitals, a cumulant approxi-
mation of F.4-RDM⁶⁵ (DMRG-cu(4)-CASPT2) has been
implemented⁶¹ and will be available in a future release.

F. Selection of the active orbital space

Multiconfigurational approaches require a selection of
the active space.⁶⁶ Traditionally, this procedure is de-
scribed as a set of recipes⁶⁷, but this set appears ever
expanding and following these rules requires some skills
from the user side. In the following section we described
the recent updates of the tools, which can simplify this
uneasy step in multiconfigurational calculations.

1. From small-basis active orbital selection to large-basis CAS/RAS with EXPBAS

Manual identification of active spaces is usually sim-
pler in a small basis set, largely due to the fact that the
virtual orbitals become more localized so that “chemi-
cal” interpretation is easier. In addition, since the com-
putational time for any quantum chemical method scales
with the size of the basis set, using a smaller basis allows
for a faster diagnosis of the stability of the chosen ac-
tive space. [Open]Molcas contains a module, called EXP-
BAS, which expands a wavefunction from a smaller basis
set into a larger one, giving the user the ability to first
find an active space in a small basis—either manually or
through the AUTOCAS algorithm described next, before
choosing a larger, more expensive basis, for production
calculations.

An example of how to use EXPBAS on a model
[Co(H₂O)(OH)₃]-system is given in the Supporting In-
formation. This system serves as a minimal model for
the interaction between a water molecule adsorbed on
a Co₂O₃-surface. In the Supporting Information, we de-
scribe how to first select a small as possible ANO-type
basis set for generating an initial guess for an active
space consisting of the water molecule's valence orbitals
as well as the Co 3d- and 3d'/4d-orbitals. Further, we
provide input examples for how to expand the wavefunc-
tion from the smaller basis into a larger one using EXP-
BAS and finish of with a small comparison between the
EXPBAS-strategy for identifying an active space compared
to starting with guess orbitals generated by GUESSORB.

2. AUTOCAS facilitates automated orbital selection

Since choosing an appropriate active orbital space for
a problem at hand is a non-trivial task, it ultimately
calls for a high degree of automation—already to overcome expertism and potential anthropogenic bias. To this end, QCMAQUIS and OpenMolcas can be steered by the graphical user interface AUTOCAS which provides a fully automated active orbital space selection protocol for multiconfigurational wavefunctions⁶⁸,⁶⁹. We developed this protocol⁷⁰,⁷¹ to cope with valence properties, for which it inspects the whole valence orbital space of a molecule in an approximate, but comparatively fast DMRG calculation to identify the strongly statically correlated orbitals. We note that the protocol will also work for other properties (e.g., Rydberg excited states) if relevant orbitals are also considered in this fast exploratory DMRG calculation.

The AUTOCAS selection algorithm based on orbital entanglement measures replaces empirical selection guidelines by physical quantities that are calculated from an approximate wave function—a procedure that eliminates arbitrariness and enhances reproducibility. The selection threshold⁷⁰ of the protocol is defined such that only the most strongly statically correlated orbitals are selected for the active space which yields compact active spaces. This guarantees a separation of the static and dynamic electron correlation, where the latter can be accounted for by subsequent CASPT2 or NEVPT2 calculations within OpenMolcas. For instance, the AUTOCAS selection algorithm reliably identifies relevant double-shell orbitals of 3d-transition metals⁷², confirming that suitable zeroth-order DMRG-SCF or CASSCF wave functions can be constructed from the automatically selected active spaces.

In cases where the active space is constructed only from valence orbitals, we recommend an automated active space selection in a minimal basis with AUTOCAS and subsequent expansion of these pre-optimized active orbitals to the final large basis set with EXPBAS.

### III. TYPICAL PROBLEMS SOLVED BY [OPEN]MOLCAS

In the following section several cases of the use of [Open]Molcas package in typical quantum chemical applications are described. In Supplementary information user can find corresponding inputs.

#### A. Precise wave function calculations

##### 1. RASSCF/RASPT2 electronic structure calculation

An incomplete time line illustrates the evolution from the early times when only CASPT2 was available for dynamic correlation of CASSCF wave functions. This combination was then used extensively for molecules with lighter elements up to transition metal compounds. Spin-orbit interaction could already be treated⁷³, e.g. in the

<figure><img src="image_3.png" /></figure>

<table><caption>FIG. 1. DMRG-CASPT2 standard binding enthalpy $\Delta H^o_{DMRG-CASPT2}$ (in kcal/mol) of O$_2$ to M$^{II}$P (M = Fe, Co) and $[Co^{II}(corrin)(Im)]^+$ . Experimental values were obtained for four-coordinate ferrous and cobaltous porphyrin sites in metal-organic frameworks.⁸²,⁸³</caption><thead><tr><th></th><th>Fe<sup>II</sup>P-O<sub>2</sub></th><th>Co<sup>II</sup>P-O<sub>2</sub></th><th>[Co<sup>II</sup>(corrin)(Im)-O<sub>2</sub>]<sup>+</sup></th></tr></thead><tbody><tr><td>ΔH<sup>o</sup><sub>DMRG-CASPT2</sub><br>(kcal/mol)</td><td>8.2</td><td>4.1</td><td>9.3</td></tr><tr><td>ΔH<sup>o</sup><sub>exp</sub><br>(kcal/mol)</td><td>8.1</td><td>3.6</td><td>-</td></tr></tbody></table>

study of EPR g-tensors⁷⁴, and by the Douglas-Kroll-Hess scalar relativistic integrals by⁷⁵ also calculations involving heavier elements including up to actinides became possible⁷⁶⁻⁷⁸. Using RASSCF, larger active spaces were possible, but technical complications did not allow a proper RASPT2 program to be developed. Ultimately, an approximate form of RASPT2 started being used experimentally⁷⁹, and it was accepted into Molcas. A benchmark study by Sauri et al.⁸⁰ showed the general usability of the approach.

Among the large number of more recent applications, we can mention the chromium dimer study⁸¹ as a recent example of the capability of the plain RASSCF/RASPT2. A number of other examples can be found in the papers³⁻⁵ and in other contributions in this paper.

##### 2. Electronic Structure and Energetic Properties Calculated with DMRG-CASPT2

To demonstrate typical problems that can be solved by DMRG-CASPT2 in [Open]Molcas, we extend a previous work⁶⁰ by studying the electronic structure and energetic properties of O$_2$ adducts of three transition metal macrocycles, i.e. $^3\text{Fe}^{II}\text{P}$, $^2\text{Co}^{II}\text{P}$ (P = porphin), and $^2[\text{Co}^{II}(\text{corrin})(\text{Im})]^+$ (Im = imidazole) (Figure 1). The latter complex, serving as a model for vitamin B$_{12}$r,⁸⁴ is of particular interest considering its potential application in oxygen reduction reaction (ORR).⁸⁵ The active space consisting of 19 active orbitals is computationally prohibited with conventional CASSCF-CASPT2 but can be easily treated with DMRG-CASPT2. The standard binding enthalpies predicted by DMRG-CASPT2 (including DFT thermal correction) are in excellent agreement with experimental data.⁸²,⁸³ The binding enthalpy of O$_2$ to $^2[\text{Co}^{II}(\text{corrin})(\text{Im})]^+$ (9.3 kcal/mol) is similar to experi-
mental values measured for five-coordinate CoN4 com-
plexes (7.8 kcal/mol to 13.1 kcal/mol).<sup>86</sup> This high bind-
ing enthalpy explains the capability of vitamin B<sub>12</sub> cat-
alyzing the four-electron ORR.<sup>85</sup>

3. Strategies for generating excited states

Excited states (ES) are often calculated using the
state-average formalism, which avoids separate optimiza-
tions of each state while ensuring a balanced description
of all states. In case the excited state belongs to a dif-
ferent irreducible representation, these states can be op-
timized separately by taking advantage of point-group
symmetry. States with different spin multiplicity are al-
ways optimized separately. [Open]Molcas makes it pos-
sible to further control the position of holes using re-
stricted active space wave functions.<sup>26,28,79</sup> Together with
energy-penalty terms or core-valence separation (CVS),
this makes it possible to target single and double hole
states, even in deep core orbitals.<sup>87,88</sup>

The exploitation of the RASSCF/RASPT2 protocol
benefiting from parallel execution, point-group symmetry
and integral calculation speed up via Cholesky decompo-
sition, has facilitated the calculation of the manifold of
higher lying ES with an accuracy which allows a direct
comparison with experiments.<sup>5</sup>

4. State interactions with RASSI

A main point with most multiconfigurational calcula-
tions is that the orbitals are optimized for the state of in-
terest. However, in calculations involving several or many
excited states, it is frequently difficult, sometimes almost
impossible, to optimize orbitals specifically for each state.
In most cases, orbitals are optimized for a set of states,
which is a stable procedure if this set contains all the
lowest states of any particular symmetry. If the molecule
has a higher symmetry (like e. g. C₃ᵥ) which causes de-
generacy of any states (like states with open e orbitals),
it is also necessary to make sure that the degenerate com-
ponents do not fall into different irreps of the actual sym-
metry used for the calculation.

In any case, the result is then subsets of states with
their own specific set of optimized orbitals, and it is usu-
ally necessary to compute various matrix elements which
couple wave functions with more or less different orbitals.
For transition amplitudes, it is necessary that the two
states are orthonormal and non-interacting, which is not
automatic if the states of a multiconfigurational calcula-
tion are obtained with different orbitals.

In e.g. MRCI, the states are obtained from a large
set of configurations, all described using one common
set of orbitals, and the problem disappears. In order to
get the advantage in MCSCF of having orbitals more or
less adapted to each state, the problem can be addressed
efficiently by the RASSI procedure. Its principle is the

ability for any two states to compute transition density
one- and two-electron matrices in the common set of AO
basis functions, and these are then combined with AO
integrals to form any matrix elements needed.

This ability has also been extended for special calcula-
tions, e.g. the Dyson amplitudes to model ionization pro-
cesses, and the capability to include spin-orbit coupling
in spin-free RASSCF calculations. Also, for higher accu-
racy of interactions with UV and X-ray photons, Open-
Molcas can compute transition probabilities using not
only the dipole (length or velocity) approximation, but
also using higher order or exponential forms¹⁷,⁸⁹ of the
electron-photon interaction. RASSI will simply combine
the transition densities with different molecular integrals
provided by SEWARD.

5. Computing excitation energies with large active spaces

As an example of what is described in Sections III A 3
and III A 4, we present here the strategy applied to com-
pute the vertical excitation energies of bacteriochlorophill
*a* (BChl) molecules. BChls, together with carotenoids,
represent the main chromophores in many light harvest-
ing complexes, such as the much studied LH1 and LH2
systems.⁹⁰ BChls are porphyrin-like molecules, which
present an extended π conjugated system surrounding
a central magnesium atom.

The complete π–π* system of a BChl unit extends over
24 molecular orbitals, containing 26 electrons. One empty
3p magnesium atomic orbital is also, often, considered as
part of the orbitals available to accept electronic density
upon photo-excitation, bringing the total number of, po-
tentially, active orbitals to 25. A complete active space³¹
treatment of such system represents a daunting task,
since it would require the handling of 3 863 302 870 000
configuration state functions (CSFs). Instead, the pur-
sued strategy was that of employing a restricted active
space²⁸ methodology, followed by a second order pertur-
bation theory correction.⁷⁹ However, the choice of which
orbitals to include in each RAS subspace is not straight-
forward.

The first step was that of employing a simple “singles
and doubles” calculation, i.e., to include 13 fully occupied
orbitals in RAS1, 0 orbitals in RAS2, 12 empty orbitals in
RAS3, and allowing 2 holes in RAS1 and 2 excitations in
RAS3, together with a state averaging over 2 electronic
states (ground and first electronically excited). Subse-
quently, the orbitals of the obtained trial wavefunction
(product of a linear combination of only 12 403 CSFs)
were ordered by their average occupation. After orbital
re-ordering, the final wavefunction was computed includ-
ing 11 orbitals in RAS1, 4 orbitals in RAS2, 10 orbitals
in RAS3, 3 holes and 3 excitations in RAS1 and RAS3,
respectively (for a total of 10 203 265 CSFs), and state
averaging over 2 states.

The energies of the ground and excited states were re-
evaluated at the multi-state RASPT2 (MS-RASPT2)⁸⁰
level of theory, employing a technique allowing the splitting of the calculation of each root correction and the corresponding Hamiltonian matrix elements into as many separate calculation as the number of roots.⁴ The so-obtained excitation energies, transition dipole moments and oscillator strengths were in line with what was expected from experimental data, while the obtained configuration interaction coefficients demonstrated the need for a multi-reference method for the accurate treatment of such large systems.⁹¹ Example inputs are provided in Section S3 as supporting material. The same methodology was successfully applied also to BChls of a different light harvesting system, namely LH3.⁹² Furthermore, the effects on the excitation energies due to neighboring amino acidic residues were evaluated by including also one amino acid moiety in the calculation.⁹³ While such calculation did not increase the number of CSFs, it increased the number of basis functions (from the starting 852 for BChl using ANO-RCC double zeta)⁹⁴ up to 1075 after addition of a tyrosine moiety. Finally, a similar procedure was successfully employed to compute the excitation energy to the second excited state, which however required the calculation of 5 states so as to resolve excited states mixing.⁹⁵

## 6. Lanthanide-activated luminescent materials

Inorganic phosphors made of solids optically activated with transition metal and, especially, lanthanide ions, are the subject of intense scientific and technological research because of their key role in a broad range of solid-state devices with high societal demand, like energy-efficient solid-state lighting devices, displays, ultraviolet to infrared solid-state lasers, scintillating detectors for medical imaging, security and high-energy physics calorimetry, remote pressure and temperature measurement systems, solar cell enhancers, energy storage phosphors, persistent luminescence, quantum information processing, and more.⁹⁶

All these applications are based on the rich photon engineering that can be built on the complex manifolds of excited states of the lanthanide ions in the solid hosts. Dense sets of tens to hundreds of excited states of different natures separated with gaps, and their different radiative and non-radiative transition probabilities, offer endless practical single- and multi-photon absorption and emission possibilities. Knowledge of the manifolds from experiments alone is limited and information from *ab initio* calculations is increasingly demanded.

In these complex manifolds of heavy elements in solids, states with partial filling of shells of molecular orbitals of several kinds co-exist, like $4f^n$, $4f^{n-1}5d$, and $4f^{n-1}6s$ configurations, together with impurity-trapped-exciton configurations (ITE) $4f^{n-1}\phi_{ITE}$, and all sort of charge transfer states, like ligand-to-metal (LMCT), inter-valence (IVCT), metal-to-metal (MMCT), or compensator-to-metal (CMCT) states. Hence, the mul-

tireference multiconfigurational wave function techniques of [Open]Molcas and its scalar and spin-orbit coupling relativistic Hamiltonians, are the right instruments for the theoretical tackling of these manifolds, provided they are used together with the *ab initio* model potential (AIMP) embedded tools that take properly into account the quantum mechanical interactions between the active center and the solid host.⁹⁷

Recent applications involving halide, oxide and sulfide hosts, and the complex manifolds of lanthanides like Pr, Eu, Tm, or Yb, together with the archetipely $R_1$ line of $Cr^{3+}$, are found in Refs. 98–102.

## B. Magnetic properties of mono- and poly- nuclear compounds

In recent years, the [Open]Molcas package has been successfully employed for describing magnetic properties of mono- and polynuclear compounds containing transition metals and/or lanthanides and actinides. In this section, the main functions allowing these studies are reviewed alongside with some results. Practical examples how to run such calculations are given in Supplementary Information. In all these studies— with the exception of the last example, which rested on DMRG-SCF/MPSSI/SINGLE_ANISO calculations⁴⁷ —the computational methodology employed was based on the RASSCF/RASPT2/RASSI/SINGLE_ANISO + POLY_ANISO calculations. A more detailed description of the capabilities of [Open]Molcas for magnetic properties can be found in references 4 and 5.

## 1. Ab initio computation of the parameters of Spin Hamiltonians for any dimension of the pseudospin

Experimental data obtained in electron paramagnetic resonance (EPR) measurements is commonly interpreted using a Spin Hamiltonian formalism¹⁰³. In the limit of weak effects of the spin-orbit coupling effects, the Spin Hamiltonian may be written as:

$$ \hat{H}_{EPR} = \mathbf{B} \cdot \mathbf{g} \cdot \tilde{\mathbf{S}} + \tilde{\mathbf{S}} \cdot \mathbf{D} \cdot \tilde{\mathbf{S}} + A \tilde{\mathbf{I}} \cdot \tilde{\mathbf{S}}, \quad (5) $$

where $\tilde{\mathbf{S}}$ is the effective spin (pseudospin) defined in the basis of the low-lying states. The first term in Eq. 5 is the Zeeman splitting, where $\mathbf{B}$ represents the applied field, while $\mathbf{g}$ is a $3 \times 3$ tensor describing the evolution of the energy states in applied field. The second term is the Zero-Field Splitting (ZFS) (for $\tilde{\mathbf{S}} \ge 1$), where $\mathbf{D}$ is a $3 \times 3$ tensor describing the splitting of the levels in the absence of the applied magnetic field. The last term in Eq. 5 is the hyperfine interaction, where $\tilde{\mathbf{I}}$ is the spin of the nucleus, $\mathbf{A}$ is the parameter describing the interaction between nuclear spin and the ground electronic spin. It is common to extract parameters for Eq. 5 by fitting experimental EPR spectra. It is noted here that the general form of Spin Hamiltonian equations may include tensorial terms
up to rank 2S̃ for Zeeman splitting and for zero-field splitting operators.

Comparison between the experimentally fitted parameters of the Spin Hamiltonians and the ab initio extracted parameters represents an active research area. In MOLCAS, two implementations exist for the computation of the g-tensor. One implementation is done inside the RASSI module (EPRG) allowing the computation of the g-tensor for a Kramers Doublet<sup>74</sup>. Another implementation, inside the SINGLE_ANISO module (MLTP), is more general, allowing the computation of the g-tensor for any size of the pseudospin (e.g. for triplet, quadruplet, etc.). The details of this implementation are given in<sup>104</sup>. Also, the SINGLE_ANISO implementation allows extraction of the parameters of the Zero-Field Splitting (i.e. the D-tensor) for any size of the pseudospin, as well as the parameters of higher-rank tensorial operators of Zeeman splitting and zero-field splitting<sup>104</sup>. In addition, the sign of the product of the main values g<sub>X</sub>g<sub>Y</sub>g<sub>Z</sub> of the g-tensor is computed<sup>105</sup>. It was shown<sup>106</sup> that this sign determines the direction of precession of the magnetic moment around the applied field.

## 2. Ab initio crystal field for lanthanides and transition metal compounds

The crystal field of lanthanide ions in complexes or in crystals is one of the main characteristics defining the structure of luminescence bands and magnetic properties<sup>107</sup>. Since [Open]Molcas allows direct computation of the electronic structure of lanthanide compounds by a RASSCF / RASPT2 / RASSI computational scheme with satisfying accuracy, a methodology for the first-principle derivation of crystal-field parameters for lanthanides was implemented inside SINGLE_ANISO. Two schemes were implemented: i) J-multiplet specific crystal field and ii) L-term specific crystal field. This implementation is based on the following steps: i) performing a rigorous ab initio calculation on a desired lanthanide / transition metal compound, ii) diagonalization of the (Z-component of the) magnetic moment in the basis of the eigenstates of the free-ion multiplet J or free ion term L, iii) put in correspondence ab initio states and the eigenstates of the magnetic moment, iv) transformation of the (diagonal) RASSI energy matrix to a new basis, where the Z-component of the magnetic moment is diagonal and each of the states is characterized with a definite projection |M<sub>J</sub>⟩ / |M<sub>L</sub>⟩ on the quantization axis, v) apply the irreducible tensor operator method to find all crystal field parameters.

The crystal field Hamiltonian is represented by:

$$ \hat{H}_{\text{CF}} = \sum_{k=0}^{2S \text{ or } 2L} \sum_{q=-k}^{+k} B_k^q \hat{O}_k^q, \quad (6) $$

where the $B_k^q$ are the parameters extracted from the performed ab initio calculation, $\hat{O}_k^q$ represent the irreducible

tensor operator of rank k and projection q which is written in the basis of $2J + 1$ or $2L + 1$ states. The extracted parameters $B_k^q$ allow to recompute the complete energy matrix $\hat{H}_{\text{CF}}$ via Eq. 6, which gives after diagonalization, the original RASSI or CASSCF energies, eigenstates and their properties as they were initially obtained in the performed ab initio calculation.

For both methods described above, the extracted parameters are given in several operator representations: Extended Stevens Operators<sup>108</sup>, Chibotaru-ITO<sup>104</sup>, Iwahara,<sup>109</sup> etc. In each case, a statistical analysis of the role played by each parameter is given, such that the importance of each parameter for the total crystal field splitting is quantified.

## 3. Blocking barriers of mono- and poly-nuclear single-molecule magnets (SMM)

Single molecule magnets are paramagnetic molecules which are able to preserve their intrinsic magnetization for a long time, after the applied magnetic field is switched off. The performance of the single-molecule magnets is indicated by the temperature at which the magnetic hysteresis is closing (blocking temperature), and the height of the magnetization blocking barrier. Magnetization blocking phenomenon is directly related to the axiality of the crystal field (i.e. the dominance of the axial crystal field parameters over non-axial ones, at least for the ground multiplet or term) and magnetic axiality of the ground and excited doublet states (i.e. the dominance of the main value g<sub>Z</sub> over the other two components, g<sub>X</sub> and g<sub>Y</sub>). The requirements and routes how to obtain performant molecular magnets have been previously discussed, e.g. Refs. 110,111, while it is an active research area<sup>112-119</sup>.

In the case of mono-nuclear single-molecule magnets, the direct application of the ab initio computational scheme based on RASSCF / CASPT2 / RASSI / SINGLE_ANISO allows for a direct evaluation of the crystal field for lanthanides and transition metal compounds<sup>107</sup> and of the g-tensor in the ground and excited states. The computational accuracy of this method proved sufficient to understand the origin of the slow magnetic relaxation, the subtle differences between similar compounds, etc.

In particular, the performed ab initio calculations allow for direct evaluation of the blocking barrier for single-molecule magnets. As example, Figure 2 shows the structure and blocking barrier of today's record holder for the best molecular magnet, computed with MOLCAS<sup>120</sup>, while this research area benefited a lot from this computational strategy.
<figure><img src="image_4.png" /><figcaption>FIG. 2. Magnetization blocking barrier of the best single-molecule magnet obtained to this date as computed with the CASSCF / CASPT2 / RASSI / SINGLE_ANISO computational tools¹²⁰. The short black lines represent spin-orbit eigenstates placed in the plot at the *ab initio* energy (y axis) and the intrinsic magnetic moment (x axis). The intensity of the blue and red arrows indicate value of the average transition dipole moment between the corresponding states. Similar plots are generated automatically by the PLOT key of SINGLE_ANISO.⁵.</figcaption></figure>

## 4. Thermodynamic magnetic functions: susceptibility, magnetization, torque

Thermodynamic functions like field- and temperature-dependent molar magnetization vector $\vec{M}_{\alpha}(\vec{B}, T)$, molar magnetic susceptibility tensor $\chi_{\alpha\beta}$, powder averaged susceptibility and powder magnetization, field-, temperature- and direction-dependent magnetization torque function $\tau_{\alpha}(\vec{B}, T)$, field- and temperature-dependent magnetic susceptibility $\chi_{\alpha\beta}(\vec{B}, T)$ are available inside the SINGLE_ANISO module.

From the computation point of view, the entire SINGLE_ANISO module is much faster than other calculations discussed here, usually done within a few minutes in most of the cases. Only the computation of the powder magnetization is the most time consuming, as the Zeeman Hamiltonian has to be built and diagonalized independently for a large number of field directions and field strength points.

## 5. Semi-ab initio description of magnetic exchange interaction and the simulation of magnetic properties of polynuclear compounds

Direct *ab initio* determination of the energy spectra and properties of polynuclear compounds containing lanthanides and/or transition metals with unpaired spins is still a challenge for the current computational methods. The difficulties arise due to several factors. First, the active space of the CASSCF method must comprise all frontier orbital shells (e.g., all sets of d and f shells) containing unpaired electrons. Second issue is related to the number of states required to be optimized for a decent description of the spin-orbit interaction. Last, the numerical accuracy of the dynamical correlation methods, like MRCI or CASPT2, would need to provide very accurate energy dif-

ferences, within 1 cm⁻¹ or lower, if possible. Obviously, it is not yet possible to perform such calculations routinely for any compound of interest. To this end, in order to tackle the electronic structure and magnetic behavior of polynuclear compounds, a semi-*ab initio* approach was proposed and implemented in POLY_ANISO¹²¹,¹²². In brief, the method consists of the following steps:

* For a compound consisting of *N* magnetic metals sites, build *N* mononuclear fragments, each fragment containing unpaired electrons only on one metal site. This is achieved by computational substitution of all other (magnetic) metal sites by their diamagnetic equivalents (e.g. $Dy^{3+} \rightarrow Lu^{3+}$). Ideally, the ligand framework is left unaltered in all fragments.

* Perform *ab initio* calculation for each fragment at the desired level of theory (e.g., CASSCF / CASPT2 / RASSI / SINGLE_ANISO). For each fragment, the local spin-orbit spectra, spin and magnetic moments are obtained.

* The total magnetic exchange and magnetic dipole-dipole interaction Hamiltonians are built and diagonalized in the basis of the products of the local *ab initio* wave functions obtained for the individual metal sites at previous step. The magnetic exchange interactions between metal fragments are included in an effective (phenomenological) way using the Lines approximation¹²³. Magnetic dipole-dipole interaction is added exactly, given that all local magnetic moments are available.

* Compute the magnetic properties of the entire polynuclear compound using the obtained exchange (coupled) eigenstates

* Extract (fit) the parameters of the inter-site exchange interaction from the direct comparisons between computed and measured magnetic data.

The algorithm was implemented in the POLY_ANISO software and is available as a module in the MOLCAS package and also as a stand-alone program. Figure 3 shows a particular example of application of the above described computational approach for the investigation of the magnetic properties of $Dy_3$ triangles¹²².

Recently, the projection of exchange interaction on products of irreducible tensor operators was implemented in the POLY_ANISO. This function allows a full derivation of all parameters of the multipolar magnetic exchange, in a way similar to the extraction of the crystal field. The feature was first reported in Ref. 5.

## 6. Magnetic properties of an actinide complex calculated with MPSSI

In this paragraph we will discuss the first application of the MPSSI approach⁴⁷ described in Section II E 1 that
<figure><img src="image_5.png" /><figcaption>FIG. 3. a) Energy structure on individual Dy sites in Dy$_{3}$ molecule; b) red dashed lines represent the local main anisotropy axes of the ground doublet states of each Dy site with respect to molecular frame. The blue arrows represent one of the two components of the ground state of the en- tire Dy$_{3}$ molecule; c) comparison between measured and cal- culated molar magnetization of the entire Dy$_{3}$ triangle; d) comparison between measured and calculated molar magnetic susceptibility of the entire Dy$_{3}$ triangle$^{122}$.</figcaption></figure>

makes a combined use of both the AUTOCAS function-
ality described in Section IIF 2 in order to rationalize
the selected active orbital space and the SINGLE_ANISO
approach described above to calculate magnetic prop-
erties. In particular, we consider the magnetic proper-
ties of a hexakisamidouranium(V) complex, [U(NH$_{2}$)$_{6}$]$^{-1$
which serves as a model complex for a unique, pseudo-
octahedral hexakisamidouranium(V) complex synthe-
sized and experimentally characterized by K. Meyer et
al.$^{124}$.

The active orbital space consisting of seven electrons
distributed in 10 active orbitals, CAS(7,10), is sufficiently
small to facilitate a direct comparison of the DMRG-
SCF/MPSSI data with the corresponding numbers ob-
tained from CASSCF/CASSI calculations. Dynamical
correlation effects in the MPSSI (CASSI) calculations
were taken into account by standard multi-state CASPT2
calculations for the seven lowest-lying electronic states
of [U(NH2)6]1 which also define the ensemble of states
for the preceding state-averaged DMRG-SCF (CASSCF)
calculations. Accordingly, the number of renormalized
block states m which determines the accuracy of the
DMRG-SCF calculations was chosen such as to exactly
reproduce the corresponding CASSCF calculations for
each of the active orbital spaces considered in this ex-
ample.

Figure 4 illustrates the state-averaged orbital entan-
glement diagram obtained by means of the AUTOCAS

program based on DMRG-SCF calculations with a
CAS(7,10). Here, orbitals #4–10 correspond to predom-
inantly 5f orbitals of the U central ion while orbitals
#1–3 exhibit mainly U 6p character with little contri-
butions from each of the N center of the (NH2)⁻ amido
ligands. Interestingly, orbital #8 (of U 5f<sub>s</sub> character)
which displays a seemingly small orbital entanglement,
will be automatically included by the automatic selec-
tion protocol because it is singly occupied in the elec-
tronic ground state. The latter finding is in perfect agree-
ment with conclusions drawn from scalar relativistic DFT
calculations for the same model complex¹²⁴. The situa-
tion is, however, less clear for orbitals #1–3. Inclusion
of the latter in the zeroth-order wave function (corre-
sponding to a CAS(7,10)) leads to a lowering of the MS-
CASPT2 excitation energies for the two sets of lowest-
lying three-fold degenerate spin-free states of approx. 10
and 5 %, respectively. More importantly, the isotropic
ground-state *g*-factor of this highly-symmetric molecular
complex changes from *g*<sub>xx</sub> = *g*<sub>yy</sub> = *g*<sub>zz</sub> = −1.06 to −1.10
in going from a CAS(1,7) to a CAS(7,10) zeroth-order
DMRG-SCF (or equivalently CASSCF) wavefunction.
The latter *g*-value of −1.10 obtained with our DMRG-
SCF/MS-PT2/MPSSI approach is in very good agree-
ment with the measured isotropic |*g*|-factor of 1.12 for
the hexakisamidouranium(V) complex of Ref. 124. More-
over, the calculated DMRG-SCF/MS-PT2/MPSSI effec-
tive magnetic moment *μ*<sub>eff</sub> = 0.95 *μ*<sub>B</sub> (within a tem-
perature range of 5 K to 35 K) of the present hexak-
isamidouranium(V) model complex [U(NH<sub>2</sub>)<sub>6</sub>]<sup>−1</sup> agrees
reasonably well with the corresponding effective mag-
netic moment *μ*<sub>eff</sub> = 1.16 *μ*<sub>B</sub> experimentally determined
in SQUID measurements within the same temperature
range for the larger hexakisamidouranium(V) complex
of Ref. 124.

C. Beyond single point calculations: Geometry optimization, Potential Energy Surfaces, Molecular Dynamics

Quantum chemistry calculations are rarely reduced to
computing the energy or other properties of a system at
some predefined structure. More often than not, the user
would like to optimize the structure, find other significant
points on the potential energy surface (saddle points,
crossing points...), obtain plausible reaction mechanisms,
simulate the change of structure with time, or otherwise
explore the relation between structure and properties.
MOLCAS, like many other quantum chemistry packages,
includes some tools to allow and facilitate these tasks.

1. Energy derivatives

The key property for exploring potential energy sur-
faces is the energy gradient, or the derivative of the
energy with respect to the nuclear coordinates. In
<figure><img src="image_6.png" /><figcaption>FIG. 4. Orbital entanglement diagram for the final active space, CAS(7,10), averaged over the seven lowest doublet electronic states of $[U(NH_2)_6]^{-1}$ generated with AUTOCAS. Orbital #8 is singly occupied in the electronic ground state and hence automatically included by the automatic selection protocol despite its small orbital entanglement. Exclusion of orbitals #1-3 with similarly small entanglement has a minor effect on the energetics of the seven lowest-lying electronic states as well as on the calculated magnetic properties.</figcaption></figure>

[Open] Molcas, it is possible to compute analytical gradients for many wave function methods, including: SCF (HF or DFT), MP2,¹²⁵ CASSCF (state-specific or state-averaged,⁴⁹) MC-PDFT,¹²⁶ DMRG-SCF; also when using Cholesky-decomposed two-electron integrals in the acCD flavor.²²,¹²⁵,¹²⁷⁻¹²⁹ For state-averaged CASSCF, analytical nonadiabatic couplings are available too.²¹ The advantage of analytical gradients is that their computational cost is usually of the order of magnitude of an energy calculation. Even in cases where [Open] Molcas does not currently support analytical gradients (notably CASPT2), numerical differentiation is supported and works transparently for the user.

Similarly, the second derivatives of the energy are available analytically for some variational methods (HF, state-specific CASSCF), and can be computed numerically for the rest. Second derivatives provide a means to obtain vibrational frequencies and thermochemical properties, or characterize critical points on a potential energy surface.

## 2. Geometry optimization

Geometry optimizations are performed by using nonredundant internal coordinates to represent the molecular structure, thus eliminating translational and rotational degrees of freedom whenever the energy is invariant to them. By default, the internal coordinates are generated automatically based on empirical force constants,¹³⁰ such that the most relevant coordinates are prioritized and the user is not required to choose them. In addition, arbitrary constraints can be specified, and they will be satisfied at the converged geometry.¹³¹ These constraints can be any combination of “primitive” constraints, including interatomic distances, angles, dihedrals constraints to hypersphere or hyperplane surface, or treating fragments as rigid bodies. Moreover, numerical energy differentiation can be made aware of constraints, reducing the number of energy calculations required.¹³²

Judicious use of constraints permits the exploration of specific regions of a potential energy surface.¹³¹ One can, for instance, perform a “relaxed scan” by setting the values of some coordinates and optimizing all the other degrees of freedom. Saddle points (transition states) can be located starting from any structure by specifying a guiding constraint that forces the optimization towards a region where a saddle point is presumed to exist; once the right curvature of the surface is detected, the constraint is automatically released and a standard saddle-point optimization proceeds. Minimum energy path optimization and intrinsic reaction coordinate analysis are performed automatically as a series of constrained optimizations. An additional energy difference constraint allows the location of minimum energy crossing points or conical intersections, and in combination with other constraints, other points or paths in the intersection seam can be optimized.²¹

## 3. Molecular dynamics

An alternative way to explore the potential energy surface (PES) is to run molecular dynamics (MD) simulations. In this type of simulation the molecular system of interest is evolving in time and its kinetic energy is explicitly taken into account. This is in contrast to a geometry optimization or other static calculations like the minimum energy path (MEP) where the kinetic energy is zero. As a result, MD allows to overcome small barriers and reaches regions of the PES which are otherwise not accessible. This is accomplished by the DYNAMIX module in [Open] Molcas that solves Newton's equations of motion, to propagate the system in time. In order to start an MD simulation the user needs to provide the initial conditions, which are the nuclear coordinates and velocities. Based on the selected electronic structure method and basis set, forces are computed which are required to extract the acceleration to solve Newton's equations of motion. The DYNAMIX module works with any elec-
<figure><img src="image_7.png" /><figcaption>FIG. 5. Molecular Dynamics simulation of proteorhodopsin</figcaption></figure>

tronic structure method in [Open]Molcas. However, it has
been more widely used to study ultrafast photochemi-
cal or chemiluminescent processes that are non-radiative
and are completed within a few hundred femtoseconds
(10<sup>-15</sup> s).<sup>133-141</sup> Photochemical processes start in the
electronic excited state and then pass non-radiatively
through a conical intersection (Fig. 5) to the ground
state, while chemiluminescent processes behave the other
way around. The non-radiative transition is modelled by
a trajectory surface hopping algorithm<sup>142</sup> that works for
multiconfigurational methods in [Open]Molcas. The sim-
ulation of these events provides a molecular level insight
into the rearrangements responsible for the excited state
relaxation.

Furthermore, the effect of external forces can also be
included. ¹³⁶ The MD can be combined with other mod-
ules in [Open]Molcas such as espf¹⁴³ to run QM/MM
molecular dynamics in complex environments. One ex-
ample is a study of the photoisomerization in the prote-
orhodopsin protein (Fig. 5) by Borin et al. ¹⁴⁴ Using the
input file (see SI) a QM/MM trajectory was calculated at
the CASSCF level of theory. Starting from the Franck-
Condon point, the trajectory evolved towards the coni-
cal intersection within 227 fs. While this conformational
change did not require much space in the tight protein
binding pocket, the completion of the *trans-to-cis* iso-
merization was hindered specifically by the interaction of
the C14 of the retinal and the amino acid Tyr200. It was
found that this is due to a twist around the C13 = C14
and C15 = N in opposite directions. Hence, it was de-
duced that tyrosine 200 is a critical amino acid for a
successful isomerization.

In addition to the native DYNAMIX module,
[Open]Molcas is interfaced to several standalone
molecular dynamics programs. One option is SHARC¹⁴⁵,
that provides the possibility to include the effect of
an external electric field and transitions between state
of different spin multiplicity¹⁴⁶. Another option is
COBRAMM,¹⁴⁷ that, with its combined quantum
mechanical/molecular mechanical (QM/MM) framework
and surface hopping scheme, allows to simulate ground
and excited state reactions in vacuum, explicit solvent

environments, or complex macromolecular systems (e.g.,
biopolymers), including non-adiabatic molecular dynam-
ics and the underlying transient spectral signatures. A
more sophisticated option is explicitly considering the
nuclear wave function and simulating quantum molecular
dynamics. This approach automatically includes pure
quantum effects such as tunneling. In this case the inter-
face of [Open]Molcas with Quantics¹⁴⁸ is used to apply
the direct dynamics variational multiconfigurational
Gaussian (DD-vMCG) method.¹⁴⁹

4. Automated QM/MM Modeling

A specialized protocol for the automated construc-
tion of hybrid quantum mechanical/molecular mechan-
ical (QM/MM) models, which uses OpenMolcas as the
electronic structure calculation engine, has been intro-
duced. This is the Automatic Rhodopsin Modeling (a-
ARM) protocol, a command-line oriented computational
tool, designed for the congruous and reproducible gener-
ation of monomeric, gas-phase and globally uncharged
QM/MM models of rhodopsins based on electrostatic
embedding and the hydrogen-link-atom frontier between
the QM and MM subsystems. ¹⁵⁰,¹⁵¹ Although a-ARM
currently only constructs rhodopsin-like models, it pro-
vides a template for future development and generation
of an automatic QM/MM building strategy for other,
more general systems. Furthermore, the tool is already
much in use for research, as rhodopsins represent one
case of a widespread family of biological photoreceptors
capable of carrying out different biological functions. In
fact, members of the rhodopsin family are found in many
diverse organisms and, thus, constitute an exception-
ally widespread class of light-responsive proteins, driv-
ing fundamental functions in vertebrates, invertebrates
and microorganisms. ¹⁵²⁻¹⁵⁴ a-ARM has been shown to
be able to generate models suitable for the prediction of
trends in photochemical properties (i.e., maximum ab-
sorption (λ<sub>max</sub><sup>a</sup>) and emission (λ<sub>max</sub><sup>f</sup>) wavelengths) of wild
type rhodopsin-like photoreceptors and their variants,
with an error bar of 3.0 kcal/mol (0.13 eV). ¹⁵⁰,¹⁵¹,¹⁵⁵⁻¹⁵⁸

The current version of the protocol¹⁵¹ comprehends
two phases, called *input file generator* and *QM/MM
model generator*, respectively (see Figure S15 in the
SI). The first phase (Figure S15A) allows the automatic
preparation of a 3D structure in PDB format, which
contains: information on the protein structure, includ-
ing the retinal (RET) chromophore and excluding mem-
brane lipids and ions; assigned protonation states of all
the ionizable residues; optimized positions of Cl⁻/Na⁺
external counterions needed to neutralize both inner and
outer protein surfaces; and a file containing the list of
amino acid residues forming the cavity hosting the RET
protonated Schiff base (rPSB). Such a structure is a suit-
able input to be used in the subsequent QM/MM model
generator phase. Given the options selected in the first
phase, *a-ARM* is sub-divided in *a-ARM<sub>default</sub>* and *a*-
ARM<sub>customized</sub>. The former refers to a fully automatic input generation, which uses default parameters as suggested by the code (i.e., chain, pH, protonation states, cavity), whereas the latter allows the customization of some of such parameters. The novelty of these two approaches is that, regardless of the user or the computational facilities, reproducible inputs, and consequently reproducible QM/MM models, are guaranteed when the same parameters are imposed.

a-ARM generates QM/MM models of rhodopsin proteins (Figure S15B in the SI) using some of the [Open]Molcas features described above, e.g., the speed-up due to the Cholesky infrastructure. Actually, a-ARM employs a particular QM/MM approach only available on the [Open]Molcas/TINKER interface (ESPF module). In the framework of the ESPF method,<sup>143</sup> the QM part of the chromophore directly interacts with the MM electrostatic potential (usually generated by permanent and sometimes induced atomic multipoles) through one-electron operators whose expectation values represent the QM charge distribution of the chromophore. It is noteworthy that the ESPF method scales with the size of the QM subsystem only, making it particularly efficient for very large molecular systems like rhodopsin proteins. As described in Ref. 150, the chromophore nonbonding and bonding interactions are modeled based on AMBER94 rules, where a customized set of AMBER94 parameters is used for the van der Waals parameters for RET. Partial charges are computed as AMBER-like RESP charges, by using a different set of parameterized charges for each conformation of the RET.

Briefly, in order to obtain 10 independent QM/MM models, the *a*-ARM input is treated starting with *N* = 10 different random seeds that provide 10 independent sets of initial velocities, by the following actions (see Figure S15 in the SI): optimization of crystallographic water molecules and addition of hydrogen atoms to polar residues, using DOWSER;<sup>159</sup> molecular mechanics energy minimization and simulated annealing/molecular dynamics relaxation, employing GROMACS;<sup>160</sup> and geometry optimization and single point energy calculation at the CASSCF(12,12)/AMBER and CASPT2(12,12)/6-31G(d) [also MS-] levels respectively,<sup>161</sup> using the [Open]Molcas<sup>4</sup>/TINKER<sup>162</sup> interface. For the output 10 replicas of the final equilibrated gas-phase and globally uncharged monomer QM/MM *a*-ARM model, the first vertical excitation energies (Δ*E*<sub>S1–S0</sub>) between ground state (S<sub>0</sub>) and the first singlet excited state (S<sub>1</sub>) are calculated. The final *a*-ARM result is the average of the 10 Δ*E*<sub>S1–S0</sub> values. Detailed explanation of the *a*-ARM protocol workflow is provided in Refs 150 and 151.

The validation of the *a*-ARM protocol for the prediction of trends on Δ*E*<sub>S1–S0</sub> (λ<sup>a</sup><sub>max</sub>), was performed by using a benchmark set of 44 rhodopsins (25 wild type and 19 mutants).<sup>151</sup> The *a*-ARM<sub>default</sub> approach proved to be capable of reproducing the Δ*E*<sub>S1–S0</sub> values for 79% (35/44), with an error lower than 3.0 kcal/mol (0.13 eV), whereas the other 21% were successfully obtained with

the *a*-ARM<sub>customized</sub> approach (i.e., changing the protocation states) (see Figure 5 in Ref. 151). Figure 6 reports the current state of the art for *a*-ARM-obtained QM/MM models of rhodopsins. As shown, the *a*-ARM protocol demonstrates a very good agreement for the prediction of trends in excitation energies for rhodopsins which structure was obtained from either X-ray crystallography or comparative modeling, as well as for rhodopsin variants. Some of these *a*-ARM QM/MM models have been used as preliminary inputs for more sophisticated calculations such as pH-constant dynamics,<sup>163</sup> the construction of fully relaxed QM/MM models embedded in a biological membrane, population analysis and one/two photon absorption spectrum simulation.<sup>158</sup> Finally, *a*-ARM is available as a user-friendly web interface, called Web-ARM (web-arm.org),<sup>164</sup> which provides the same features as the *a*-ARM protocol to the general audience.

D. From states to spectra

Electronically excited states play important roles in many chemical reactions and spectroscopic techniques. Theoretical simulations of spectra make it possible to

<figure><img src="image_8.png" /><figcaption>FIG. 6. (A) Computed excitation energies Δ*E*<sub>S1–S0</sub> in both kcal/mol (left axis) and eV (right axis) for various rhodopsins (x axis; M and A indicates microbial and animal rhodopsins, respectively). The employed protein structures where obtained from X-ray crystallography (green panel, left) or through comparative modeling (white panel, center). Two sets of variants for bovine rhodopsins (Rh) and bacteriorhodopsin (bR) are also reported (red panel, right). The computed data was obtained using the ARM method (yellow dots),<sup>150</sup> *a*-ARM<sub>default</sub><sup>151,164</sup> (gray up-turned triangles) and *a*-ARM<sub>customized</sub><sup>151</sup> (red squares). Experimental data, as energy difference corresponding to the wavelength of the absorption maxima, is also reported (blue down-turned triangles). (B) Differences between computed and experimental excitation energies ΔΔ*E*<sub>S1–S0</sub> in both kcal/mol (left axis) and eV (right axis). All computed data is within a 3 kcal/mol (0.13 eV) error, apart from a number of outliers that were corrected using *a*-ARM<sub>customized</sub>. Further details can be found in Ref. 151.</figcaption></figure>
directly compare with experiments and to predict their outcomes. Spectra can either be generated through response theory on the ground state or by explicit calculations of excited states. The latter approach is particularly useful for high-level wave function methods like CASPT2 and/or multiphoton processes where response equations are not yet available. In this section the unique opportunities to study excited state properties with the RASSI module are outlined,³¹,³² followed by recent examples from transient non-linear spectroscopy and X-ray scattering.

## 1. Multichromophoric systems and Frenkel Excitonic Model

Multichromophoric complexes are conveniently studied by reconstructing the electronic Hamiltonian of the complete system in terms of independent but interacting units. A general theoretical framework to describe such systems is the so called Frenkel excitonic model, for which more details are reported in the SI to the present paper. In brief, one first performs quantum chemistry computations on single sites, i.e. subregions of the entire system, and then computes site-site coupling terms to reconstruct the multichromophoric manifold of electronic states. An OpenMolcas application of this computational model to light harvesting complexes, with site-site coupling terms computed at the transition density level (beyond the dipole-dipole approximation) has been demonstrated in recent years.¹⁶⁵

## 2. Transient non-linear spectral signatures

Transient spectroscopy is a versatile tool for resolving molecular dynamics on both the ground and the excited state. It groups an ample body of techniques that combine laser pulses ranging from the IR to the X-ray that allow to follow dynamical events from the femtosecond to the hour time scale. Its simulation requires to describe the electronic structure of the system, its coupling to the nuclear degrees of freedom and non-adiabatic effects. Regarding the electronic structure calculations, a bottleneck is associated with the computation of the manifold of excited states and the transitions between them at different molecular geometries. These are the physical observables that determine energy position and intensity of the recorded signals and reveal the system dynamics in the experimental measurements. Novel unique functionalities in OpenMolcas aim at facilitating the computation of these observables. Recently, we exploited some of these features in simulations of transient spectroscopy on a number of conjugated and aromatic compounds¹⁶⁵⁻¹⁷¹.

On the example of the ultrafast photoinduced isomerization of trans-azobenzene we present two types of transient absorption spectroscopies: one utilizing UV-Vis pulses (a)¹⁶⁶ and one utilizing a combination of UV-Vis

and X-ray (b) pulses to probe the dynamics¹⁶⁷.

(a) UV-Vis transient spectroscopy: The simulation of UV-Vis signals in transient absorption spectroscopy requires to compute the manifold of higher lying bound valence excited states.¹⁷² Despite not being directly involved in the photoinduced dynamics, they serve as state specific probes (i.e. “spectator” states) that allow to follow the system evolution in real time. The electronic structure was computed at few selected geometries (excited state minima and conical intersections) within the RASSCF/RASPT2 framework. This strategy allows to expand the active space beyond the full π-valence space, a requirement of paramount importance for obtaining quantitative results when dealing with higher lying excited states.¹⁷²,¹⁷³ Point-group symmetry ($C_{2h}$) and the parallelized Cholesky decomposition facilitated the simultaneous treatment of tens of excited states in the state-averaging formalism. Transition dipole moments, RASPT2 energies and (numerical) gradients were then used to simulate the system time-resolved non-linear spectroscopy via the *Spectron code*.¹⁷⁴,¹⁷⁵ The resulting spectra (Figure 7(a))) show excellent agreement with state-of-the-art experiments, revealing photoinduced absorption (PA) features of the “spectator” states whose intensity beatings are an indicator of the vibrational dynamics¹⁶⁶. On the basis of these simulations we could show that near-UV irradiation of azobenzene involves a radiative-to-kinetic energy transfer into CNN in-plane bending modes, triggering an ultrafast nonproductive decay channel which explains the experimentally observed reduced quantum yield of the *trans* → *cis* isomerization.

(b) X-ray based transient spectroscopy: The photoinduced dynamics of a variety of systems can also be followed employing X-ray pulses. In NEXAFS (near edge X-ray absorption fine structure) spectroscopy these are sent on a already excited sample with a controlled time delay, inducing a transition from the core orbitals to empty valence orbitals. These core-to-valence transitions may be used to reveal the system dynamics in a state-specific way.¹⁷⁶ To this aim, we explored the azobenzene $n\pi^*$ state dynamics by computing the N1s → valence transitions at selected geometries along the minimum energy path connecting the Franck-Condon point and the lowest energy conical intersection with the ground state. We made use of the recently implemented projection technique called highly excited states (HEXS) which allows the selective target of core-hole configuratons and the straightforward calculation of the transition dipole moments between valence and core-excited states (via the RASSI program).¹⁷¹

The simulated spectra, shown in Figure 7(b)), demonstrate the sensitivity of NEXAFS to the key molecular modes involved in the azobenzene isomerization on the $n\pi^*$ surface: eventually we showed that some specific pre-K-edge features are particularly suitable for tracking coherent ultrafast CI-mediated non-adiabatic dynamics, as they show molecular mode specific changes, and appear
<figure><img src="image_9.png" /></figure>

FIG. 7. Ultrafast photoinduced dynamics of the trans-azobenzene revealed by (a) transient spectroscopy in the UV-Vis spectral region¹⁶⁶, and (b) time-resolved NEXAFS spectroscopy, employing a combination of UV and X-ray pulses¹⁶⁷. In (a), the $\pi\pi^* \rightarrow n\pi^*$ internal conversion dynamics along the symmetric bending coordinate is revealed (in both theory and experiment) by the disappearance of the $\pi\pi^*$ PA signals (PA1 and PA2), accompanied by the appearance of the $n\pi^*$ PA signal. In (b) the key molecular motions (CNN symmetric and asymmetric bending and CNNC torsion) that bring the molecule from the $n\pi^*$ Franck-Condon point to the minimum energy conical intersection (CI) with the ground state evoke changes in two intense NEXAFS pre-edge features. These correspond to $1s \rightarrow n$ and $1s \rightarrow \pi^*$ transitions. Note that the shortest distance between the two peaks is reached at the CI.

in a background-free spectral window.

### 3. X-ray scattering in transition-metal complexes

Coordination complexes of first-row transition metals participate in a wide range of processes, like electron transfer and spin crossover, and can act as efficient catalysts. X-ray spectroscopy is widely used to study these complexes due to its elemental specificity and sensitivity to the electronic structure. An excellent example of the strength of the RASSI approach is modelling of the two-photon resonant inelastic X-ray scattering (RIXS) process. When performed at the metal L-edge a 2p→3d excitation is followed by decay of a 3d electron refilling the 2p hole; the energy difference between incident and emitted photons, i.e. the energy loss, gives the valence-excitation energies of the final state. The two-dimensionality of RIXS grants access to intense metal-centered transitions and spectrally separable metal-ligand excitations, which makes it a powerful probe of metal-ligand interactions.

RIXS spectra are affected by both 2p and 3d spin-orbit coupling, as well as correlation between 2p and 3d electrons and strong relaxation due to the core-hole, which makes the RASSI approach highly suitable.¹⁷⁸⁻¹⁸⁰ Further, in response theory the RIXS process is a relativistic two-photon process, but with RASSI the spectrum can be

FIG. 8. Schematic exemplification of the excited-state RIXS three-photon process: (i) Optical pumping of ground state $|0\rangle$ yields valence-excitation $|v\rangle$. (ii) X-ray absorption into a manifold of core-excited states $|c\rangle$. (iii) X-ray emission from decay of $|c\rangle$ into $|v\rangle$ or $|0\rangle$. The measured energy loss $\omega_{loss} = \omega_{in} - \omega_{out}$ in RIXS stems from the energy difference between initial state of (ii) and final state of (iii). A sample difference spectrum, with respect to RIXS from $|0\rangle$, is shown for an excited quartet state of iron(III)hexacyanide at the Fe $L_3$ edge, from RASPT2 simulations.¹⁷⁷

generated by combining individual transition dipoles moments between initial, intermediate and final states. The same procedure also works for valence excited species, effectively simulating three-photon optical pump-X-ray

<figure><img src="image_10.png" /></figure>
probe processes, as illustrated in Fig. 8. RASPT2 sim-
ulations have been used for orbital-specific mapping of
the femtosecond ligand exchange dynamics of Fe(CO)₅
in solution.¹⁸¹ Another example is the prediction of the
fingerprints of electronic, spin and structural dynamics
in photoexcited ferricyanide.¹⁷⁷ By including terms be-
yond the electric dipole approximation, metal K pre-edge
(1s→3d) RIXS can also be described.¹⁸²

4. X-ray absorption spectra of actinide complexes and excited state bonding analyses

The Natural Bond Orbital (NBO) toolkit<sup>183,184</sup> is widely popular in all branches of Chemistry and used, for instance, for chemical bonding analyses within electronic structure calculations. NBO is most often used in conjunction with DFT calculations, but it can accommodate the one-particle (spin-) density matrices from post-HF calculations as well. This includes the complex multi-reference spin-orbit (SO) states generated by the RASSI module of [Open]Molcas. Elsewhere<sup>4,185,186</sup> we described a modification of the RASSI property code to generate, from the SO states, various types of one-particle (transition) density matrices in the AO basis that are associated with different types of one-electron property operators (real spin free (type 1), imaginary or/and spin dependent (types 2,3,4,)). The type 1 density matrices, for example, can be diagonalized in order to obtain spin-independent natural orbitals (NOs) and their occupations for the SO states, which facilitates the analysis of SO effects on metal ligand bond orders and such.<sup>187</sup> An interface has been developed recently in a Python code project named *eXatomic*,<sup>188</sup> to create NBO inputs from the SO-RASSI type 1 density matrices. Initial applications of this functionality have focused on the chemical bonding in actinide (An) complexes in the context of core-excitation spectroscopy.<sup>189,190</sup>

X-ray absorption near-edge structure (XANES) exper-
iments on actinide complexes are pursued to obtain in-
sight into the metal-ligand bonding or metal oxidation
state. The interpretation of XANES spectra in terms
of An-ligand chemical bonding should be accompanied
by theoretical calculations. The accurate modelling of
XANES of An complexes requires the calculation of mul-
tireference correlated valence and core-hole wave func-
tions that account for the core-hole interactions with the
valence shell and the SO coupling (SOC). The chemical
bonding information from the complex spin-orbit cou-
pled core-hole states probed by the XANES experiments
can then be obtained in terms of real NOs of the type
mentioned above, or localised combinations such as those
generated by the NBO algorithms. A *compact* description
for the different electronic states may be best achieved
by describing each state with its own optimal set of or-
bitals. This leads to the concept of orbital relaxation,
for instance upon an electronic excitation, which may al-
ternatively be viewed as the impact of the differential

electron correlation between the two states.

RASSCF XANES calculations with [Open]Molcas have been reported for a variety of transition metal complexes<sup>178,179,191–193</sup> (see also Section III D 3). The applicability of the RASSCF/RASSI protocol in modelling XANES spectra of An complexes has only recently been demonstrated.<sup>189,190</sup> Computational details are available in these references and the SI of the present article. Valence and different groups of core-hole states are ideally converged in separate, state-averaged RASSCF calculations by either exploiting symmetry and/or the HEXS keyword. The RASSI module then produces spin-orbit coupled manifolds and transition dipoles (and other multipoles) needed to generate the XANES spectrum.

The *eXatomic* [Open]Molcas–NBO interface was used in Ref. ¹⁸⁹ to rationalise the measured high-resolution An M₄.₅ XANES spectra of AnO₂⁺ (An = U, Np and Pu) molecules¹⁹⁴ in terms of the An–ligand bonding. The spectra (see Figure 3 in Ref. 189) feature a main peak assigned to transitions from An 3d mainly to non-bonding (φ, δ) 5f orbitals, and a satellite at high energy due to excitations from 3d into ligand–An 5f σ<sup>*</sup><sub>u</sub> orbitals. In the ground states, the AnO₂⁺ show increasing 5f–σ covalent bonding when going from U to Np to Pu. A concomitant stronger σ<sub>u</sub>/σ<sup>*</sup><sub>u</sub> splitting is then expected to lead to an increasing off-set of the σ<sup>*</sup><sub>u</sub> satellites from the main peaks when going from U to Pu. However, the experiments and the OpenMolcas calculations show the opposite trend. NBO analyses of the dipole-intense SO core-excited states, using the *eXatomic* [Open]Molcas–NBO interface, revealed a strong σ/σ<sup>*</sup> orbital relaxation upon excitation into the σ<sup>*</sup><sub>u</sub> orbital, such that the σ An–ligand covalency is much reduced in the core-excited state compared to the GS. These results highlight the problem of interpreting XANES and other types of absorption spectra in terms of ground-state orbitals only. The usage of *eXatomic* is illustrated in the SI.

E. Final considerations for large calculations

Implementation of multiconfigrational methods re-
quires large computational resources. In the SI we col-
lected known limitations of the different modules in
[Open]Molcas and estimated the time required for a typ-
ical run.

1. Cost of calculating spectra

Compared to ground-state calculations, the inclusion
of large numbers of excited states leads to significant
increase in cost, roughly proportional to the number
of states, even with recent improvements in the CI
algorithm.⁸⁸ As an example, for X-ray calculations of
iron(III) hexacyanide with a ten-orbital valence space,
each combination of irreducible representation and spin
multiplicity takes around 300 core hours. More than two
thirds of the time is spent on the PT2 step, while the cost of the RASSI calculation is almost negligible.¹⁹² The “packaging” of states into independent, task-farmable RASSCF+PT2 calculations therefore constitute a trivial but efficient parallelization of the simulation protocol, unavailable for methods that, unlike RASSI, rely on a shared state averaging for all the states.

## 2. Techniques for cost-effective calculations

In order to mitigate the costs of computing energies and properties by means of the highly correlated methods available in OpenMolcas, a number of techniques have been implemented.¹⁹⁵⁻¹⁹⁸ In particular, the use of the so-called frozen natural orbital (FNO) approximation¹⁹⁵,¹⁹⁶,¹⁹⁹ is a very robust technique for reducing the computational costs of any correlated calculation (e.g., CC methods, or the PT2 step of the various multiconfigurational methods). With the help of FNO, the number of virtual orbitals required to compute the dynamical correlation energy for ground and excited states can be reduced considerably—and controllably—leading to a massive reduction in the number of wave function amplitudes to be determined. In practice, this means that a calculation for which the cost is dominated by for example the PT2 step, can be effectively sped up, often by up to one order of magnitude, without sacrificing accuracy of the computed excitation energies. The keyword FNOC implements this option within CASPT2, whereas FNOM in the MBPT2 input is used when trying to benefit from this cost-effective alternative in a subsequent coupled-cluster calculation (CHCC and CHT3 modules).

## IV. CONCLUSIONS AND OUTLOOK ON FUTURE DEVELOPMENT

The change of the development model in MOLCAS code (from a code made by a closed community to an opensource project) unleashed a rapid development of new features and simplified the integration with other computational codes.

The [Open]Molcas open-source code is in many aspects a legacy code. For it to adapt to future soft- and hardware environments it need substantial rewriting and modifications. In particular the use of a multitude of different files and old algorithms makes efficient parallelization impossible.

Near-future extensions of the capabilities of the QCMAQUIS module to further enable DMRG-related techniques are broad. Some work will focus on exploiting the current capabilities of [Open]Molcas to implement efficiently the time-dependent CASSCF (TD-CASSCF) algorithm²⁰⁰⁻²⁰². Thanks to the implementation of the time-dependent DMRG algorithm within the QCMAQUIS package²⁰³, active spaces with about 20–30 orbitals could be targeted. This would pave the route

toward *ab initio* simulations of time-resolved spectra of complex molecules, including metal complexes. Moreover, we plan to extend the traditional SA-RASSCF-based algorithms designed to simulate X-ray absorption spectra¹⁹² to DMRG within the QCMAQUIS setting. In this way, L-edge absorption spectra of large metal complexes could be simulated by combining the SA formulation of DMRG⁴⁴ to optimize a large number of excited states, and the MPSSI algorithm⁴⁷ to properly account for spin-orbit effects. Further developments of QCMAQUIS in OpenMolcas will also concentrate on the improvement of (parallelized) traditional multi-configurational methods within the MPS/MPO framework including i) DMRG-SCF, ii) DMRG-NEVPT2 and DMRG-CASPT2, and iii) DMRG-MC-srDFT²⁰⁴. It is also planned to streamline and clear the use of external libraries such that it will be possible to enhance the performance of QCMAQUIS by a factor of more than ten.

New approaches for analyzing the most important orbitals or even most important individual excitations will be developed and polished. The possibility to verify an active space selection by the usage of DMRG codes in combination with graphical user interfaces (e.g. AUTO-CAS) will convert multiconfigurational calculations into automated routine and will allow to integrate these calculations into existing automated workflows.

Development of the DYNAMIX module will include an improvement of the sampling procedure to generate initial conditions for molecular dynamics simulations. Further, the trajectory surface hopping algorithm will be extended to be used with the CASPT2 method.

The [Open]Molcas environment is constantly evolving, as demonstrated by the varied possible features here described. A common denominator, and an excellence of the program suite, is definitely its capacity of treating complex molecular systems at the highest possible level of theory, both as a stand-alone and as interfaced to other computational codes. Furthermore, constant development ensures the possibility of employing [Open]Molcas to reproduce, rationalize and predict experimental data obtained employing the latest spectroscopical techniques.

## V. SUPPLEMENTARY INFORMATION

The supplementary information contains several sections and illustrates the short description of features, described in the article. It also contains ready to use examples.

The structure of the supplementary information:

• S1. Overview of main computational codes

– Short description of computational modules

– Flowchart

– Description of Molcas input language (EMIL)

• S2. Frenkel excitonic Hamiltonian
• S3. Computation of BCHL excitation energies

• S4. X-ray absorption calculations

• S5. ARM protocol

• S6. Examples of DMRG calculations

• S7. Tools to select active space

• S8. Excited states of lanthanides in solids

• S9. Example of magnetic properties calculation

• S10. Molecular dynamics calculations

• S11. RASSCF/RASPT2/RASSI methodology for heavy elements

• S12. SINGLE_ANISO calculations

• S13. Example of applying the semi-*ab initio* approach: fragmentation + POLY_ANISO

ACKNOWLEDGMENTS

F. A. acknowledges financial support by the EU-H2020 research and innovation programme under grant agreement No 654360 within the framework of the NFFA-Europe Transnational Access Activity. Part of this work has been performed thanks to computer resources provided by CINECA under the Project HPC-EUROPA3 (INFRAIA-2016-1-730897), with the support of the EC Research Innovation Action of the H2020 Programme. D.-C. S. and J. A. acknowledge support from the U.S. Department of Energy, Office of Basic Energy Sciences, Heavy Element Chemistry program, under grant DE-SC0001136. S. B. acknowledges support by the Swiss National Science Foundation (grant number P2SKP2_184034). A. B. is grateful for support by ETH Zurich (ETH Fellowship No. FEL-49 18-1), M. R. acknowledges support by the Swiss National Science Foundation (Project No. 200021_182400). L. D. V., L. P.-G., and M. Ol. acknowledge a MIUR (Ministero dell'Istruzione, dell'Università e della Ricerca) grant "Dipartimento di Eccellenza 2018 - 2022". M. Ol. acknowledges a NSF grant number CHE-CLP-1710191. M. D. and M. L. acknowledges support from the Olle Engkvist Foundation. E. D. L. and V. V. acknowledge computational resources provided by SNIC through LUNARC and NSC. T. B. P. acknowledges support from the Research Council of Norway through its Centres of Excellence scheme, project number 262695, and through Research Grant number 240698. K. P. acknowledges financial support from KU Leuven through the grant #C14/15/052. L. S. acknowledges financial support from Ministerio de Economía y Competitividad, Spain (Dirección General de Investigación y Gestión del Plan Nacional de I+D+i, grant MAT2017-83553-P). J. S.-M. acknowledges support from the EU-H2020 Marie Curie Actions (*AttoDNA*, FP8-MSCA-IF, grant no 747662). I. S.

gratefully acknowledges funding by the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (Grant No. 678169 PhotoMutant). L. U. and X. G. gratefully acknowledge scientific grants R-143-000-A80-114, R-143-000-A65-133 from the National University of Singapore. Computational resources of the NSCC (ASPIRE-1, grant 11001278) were used for this study.

VI. DATA AVAILABILITY STATEMENT

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

VII. REFERENCES

<footnote n="1">¹Wikipedia contributors, “List of quantum chemistry and solid-state physics software — Wikipedia, the free encyclopedia,” (2020), [Online; accessed 17-January-2020].</footnote>

<footnote n="2">²G. Karlström, R. Lindh, P.-Å. Malmqvist, B. O. Roos, U. Ryde, V. Veryazov, P.-O. Widmark, M. Cossi, B. Schimmelpfennig, P. Neogrady, and L. Seijo, “MOLCAS: a program package for computational chemistry,” Comput. Mater. Sci. **28**, 222–239 (2003).</footnote>

<footnote n="3">³F. Aquilante, L. De Vico, N. Ferré, G. Ghigo, P.-Å. Malmqvist, P. Neogrády, T. B. Pedersen, M. Pitoňák, M. Reiher, B. O. Roos, L. Serrano-Andrés, M. Urban, V. Veryazov, and R. Lindh, “MOLCAS 7: The Next Generation,” J. Comput. Chem. **31**, 224–247 (2010).</footnote>

<footnote n="4">⁴F. Aquilante, J. Autschbach, R. K. Carlson, L. F. Chibotaru, M. G. Delcey, L. De Vico, I. Fdez. Galván, N. Ferré, L. M. Frutos, L. Gagliardi, M. Garavelli, A. Giussani, C. E. Hoyer, G. Li Manni, H. Lischka, D. Ma, P. Å. Malmqvist, T. Müller, A. Nenov, M. Olivucci, T. B. Pedersen, D. Peng, F. Plasser, B. Pritchard, M. Reiher, I. Rivalta, I. Schapiro, J. Segarra-Martí, M. Stenrup, D. G. Truhlar, L. Ungur, A. Valentini, S. Vancoillie, V. Veryazov, V. P. Vysotskiy, O. Weingart, F. Zapata, and R. Lindh, “Molcas 8: New Capabilities for Multiconfigurational Quantum Chemical Calculations Across the Periodic Table,” J. Comput. Chem. **37**, 506–541 (2016).</footnote>

<footnote n="5">⁵I. Fdez. Galván, M. Vacher, A. Alavi, C. Angeli, F. Aquilante, J. Autschbach, J. J. Bao, S. I. Bokarev, N. A. Bogdanov, R. K. Carlson, L. F. Chibotaru, J. Creutzberg, N. Dattani, M. G. Delcey, S. S. Dong, A. Dreuw, L. Freitag, L. M. Frutos, L. Gagliardi, F. Gendron, A. Giussani, L. González, G. Grell, M. Guo, C. E. Hoyer, M. Johansson, S. Keller, S. Knecht, G. Kovačević, E. Källman, G. Li Manni, M. Lundberg, Y. Ma, S. Mai, J. P. Malhado, P. Å. Malmqvist, P. Marquetand, S. A. Mewes, J. Norell, M. Olivucci, M. Oppel, Q. M. Phung, K. Pierloot, F. Plasser, M. Reiher, A. M. Sand, I. Schapiro, P. Sharma, C. J. Stein, L. K. Sørensen, D. G. Truhlar, M. Ugandi, L. Ungur, A. Valentini, S. Vancoillie, V. Veryazov, O. Weser, T. A. Wesolowski, P.-O. Widmark, S. Wouters, A. Zech, J. P. Zobel, and R. Lindh, “Openmolcas: From source code to insight,” Journal of Chemical Theory and Computation **15**, 5925–5964 (2019).</footnote>

<footnote n="6">⁶V. Veryazov, P.-O. Widmark, L. Serrano-Andrés, R. Lindh, and B. O. Roos, “MOLCAS as a development platform for quantum chemistry softwares,” Int. J. Quantum Chem. **100**, 626–653 (2004).</footnote>

<footnote n="7">⁷G. Kovačević and V. Veryazov, “Luscus: molecular viewer and editor for MOLCAS,” Journal of Cheminformatics **7**, 16 (2015).</footnote>
<footnote n="8">⁸F. Aquilante, T. B. Pedersen, V. Veryazov, and R. Lindh, “Molcas—a software for multiconfigurational quantum chemistry calculations,” Wiley Interdiscip. Rev.: Comput. Mol. Sci. **3**, 143–149 (2013).</footnote>

<footnote n="9">⁹Molcas contributors, “Molcas manual,” (2019), [Online; accessed 17-January-2020], https://www.molcas.org/documentation/manual/.</footnote>

<footnote n="10">¹⁰F. Aquilante, M. G. Delcey, T. B. Pedersen, I. Fdez. Galván, and R. Lindh, “Inner projection techniques for the low-cost handling of two-electron integrals in quantum chemistry,” Mol. Phys. **115**, 2052–2064 (2017).</footnote>

<footnote n="11">¹¹F. Aquilante, L. Boman, J. Boström, H. Koch, T. B. Pedersen, A. Sánchez de Merás, and R. Lindh, “Cholesky decomposition techniques in electronic structure theory,” in *Challenges and Advances in Computational Chemistry and Physics*, Vol. 13, edited by M. G. Papadopoulos, R. Zalesny, P. G. Mezey, and J. Leszczynski (Springer, 2011) pp. 301–344.</footnote>

<footnote n="12">¹²S. Vancoillie, M. G. Delcey, R. Lindh, V. Vysotskiy, P.-Å. Malmqvist, and V. Veryazov, “Parallelization of a multiconfigurational perturbation theory,” Journal of computational chemistry **34**, 1937–1948 (2013).</footnote>

<footnote n="13">¹³J. Almlöf, “The molecule integral program,” (1974), University of Stockholm.</footnote>

<footnote n="14">¹⁴R. Lindh, U. Ryu, and B. Liu, “The reduced multiplication scheme of the rys quadrature and new recurrence relations for auxiliary function based two-electron integral evaluation,” The Journal of chemical physics **95**, 5889–5897 (1991).</footnote>

<footnote n="15">¹⁵R. Lindh, “The reduced multiplication scheme of the rys-gauss quadrature for 1st order integral derivatives,” Theoretica chimica acta **85**, 423–440 (1993).</footnote>

<footnote n="16">¹⁶A. Bernhardsson, R. Lindh, J. Olsen, and M. Fulscher, “A direct implementation of the second-order derivatives of multiconfigurational scf energies and an analysis of the preconditioning in the associated response equation,” Molecular Physics **96**, 617–628 (1999).</footnote>

<footnote n="17">¹⁷L. K. Sørensen, E. Kieri, S. Srivastav, M. Lundberg, and R. Lindh, “Implementation of a semiclassical light-matter interaction using the gauss-hermite quadrature: A simple alternative to the multipole expansion,” Phys. Rev. A **99**, 013419 (2019).</footnote>

<footnote n="18">¹⁸R. Lindh, J. W. Krogh, M. Schütz, and K. Hirao, “Semidirect parallel self-consistent field: the load balancing problem in the input/output intensive self-consistent field iterations,” Theoretical Chemistry Accounts **110**, 156–164 (2003).</footnote>

<footnote n="19">¹⁹T. B. Pedersen, F. Aquilante, and R. Lindh, “Density fitting with auxiliary basis sets from cholesky decompositions,” Theor. Chem. Acc. **124**, 1–10 (2009).</footnote>

<footnote n="20">²⁰F. Aquilante, L. Gagliardi, T. B. Pedersen, and R. Lindh, “Atomic cholesky decompositions: A route to unbiased auxiliary basis sets for density fitting approximation with tunable accuracy and efficiency,” J. Chem. Phys. **130**, 154107 (2009).</footnote>

<footnote n="21">²¹I. Fdez. Galván, M. G. Delcey, T. B. Pedersen, F. Aquilante, and R. Lindh, “Analytical state-average complete-active-space self-consistent field nonadiabatic coupling vectors: Implementation with density-fitted two-electron integrals and application to conical intersections,” J. Chem. Theory Comput. **12**, 3636–3653 (2016).</footnote>

<footnote n="22">²²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).</footnote>

<footnote n="23">²³P. Merlot, T. Kjærgaard, T. Helgaker, R. Lindh, F. Aquilante, S. Reine, and T. B. Pedersen, “Attractive electron-electron interactions within robust local fitting approximations,” J. Comput. Chem. **34**, 1486–1496 (2013).</footnote>

<footnote n="24">²⁴L. N. Wirz, S. S. Reine, and T. B. Pedersen, “On Resolution-of-the-Identity Electron Repulsion Integral Approximations and Variational Stability,” J. Chem. Theory Comput. **13**, 4897–4906 (2017).</footnote>

<footnote n="25">²⁵B. O. Roos, “The complete active space self-consistent field method and its applications in electronic structure calculations,” in *Advances in Chemical Physics; Ab Initio Methods in Quantum Chemistry - II*, edited by K. P. Lawley (John Wiley & Sons Ltd., Chichester, England, 1987) Chap. 69, p. 399.</footnote>

<footnote n="26">²⁶J. Olsen, B. O. Roos, P. Jorgensen, and H. J. A. Jensen, “Determinant based configuration interaction algorithms for complete and restricted configuration interaction spaces,” J. Chem. Phys. **89**, 2185–2192 (1988).</footnote>

<footnote n="27">²⁷D. Ma, G. Li Manni, and L. Gagliardi, “The generalized active space concept in multiconfigurational self-consistent field methods,” J. Chem. Phys. **135**, 044128 (2011).</footnote>

<footnote n="28">²⁸P.-Å. Malmqvist, A. Rendell, and B. O. Roos, “The restricted active space self-consistent-field method, implemented with a split graph unitary group approach,” J. Phys. Chem. **94**, 5477–5482 (1990).</footnote>

<footnote n="29">²⁹S. Battaglia and R. Lindh, “Extended dynamically weighted CASPT2: the best of two worlds,” J. Chem. Theory Comput. (2020).</footnote>

<footnote n="30">³⁰P.-Å. Malmqvist, “Calculation of transformation density matrices by nonunitary orbital transformations,” Int. J. Quantum Chem. **30**, 479 (1986).</footnote>

<footnote n="31">³¹P. Å. Malmqvist and B. O. Roos, “The CASSCF state interaction method,” Chem. Phys. Letters **155**, 189–194 (1989).</footnote>

<footnote n="32">³²P.-Å. Malmqvist, B. O. Roos, and B. Schimmelpfennig, “The Restricted Active Space (RAS) State Interaction Approach with Spin-Orbit Coupling,” Chem. Phys. Letters **357**, 230–240 (2002).</footnote>

<footnote n="33">³³P.-Å. Malmqvist and V. Veryazov, “The binatural orbitals of electronic transitions,” Molecular Physics **110**, 2455–2464 (2012).</footnote>

<footnote n="34">³⁴S. Keller, M. Dolfi, M. Troyer, and M. Reiher, “An efficient matrix product operator representation of the quantum-chemical hamiltonian,” J. Chem. Phys. **143**, 244118 (2015).</footnote>

<footnote n="35">³⁵S. Keller and M. Reiher, “Spin-adapted matrix product states and operators,” J. Chem. Phys. **144**, 134101 (2016).</footnote>

<footnote n="36">³⁶S. Knecht, E. D. Hedegaard, S. Keller, A. Kovyrshin, Y. Ma, A. Muolo, C. J. Stein, and M. Reiher, “New approaches for ab initio calculations of molecules with strong electron correlation,” Chimia **70**, 244–251 (2016).</footnote>

<footnote n="37">³⁷S. R. White, “Density matrix formulation for quantum renormalization groups,” Phys. Rev. Lett. **69**, 2863–2866 (1992).</footnote>

<footnote n="38">³⁸S. R. White, “Density-matrix algorithms for quantum renormalization groups,” Phys. Rev. B **48**, 10345–10356 (1993).</footnote>

<footnote n="39">³⁹A. Baiardi and M. Reiher, “The density matrix renormalization group in chemistry and molecular physics: Recent developments and new challenges,” J. Chem. Phys. **152**, 040903 (2020).</footnote>

<footnote n="40">⁴⁰S. Battaglia, S. Keller, and S. Knecht, “An efficient relativistic density-matrix renormalization group implementation in a matrix-product operator formulation,” J. Chem. Theory Comput. **14**, 2353–2369 (2018).</footnote>

<footnote n="41">⁴¹Y. Ma, S. Knecht, S. Keller, and M. Reiher, “Second-order self-consistent-field density-matrix renormalization group,” J. Chem. Theory Comput. **13**, 2533–2549 (2017).</footnote>

<footnote n="42">⁴²T. Dresselhaus, J. Neugebauer, S. Knecht, S. Keller, Y. Ma, and M. Reiher, “Self-consistent embedding of density-matrix renormalization group wavefunctions in a density functional environment,” J. Chem. Phys. **142**, 044111 (2015).</footnote>

<footnote n="43">⁴³Y. Ma, S. Knecht, and M. Reiher, “Multiconfigurational effects in theoretical resonance Raman spectra,” Chem. Phys. Chem. **18**, 384–393 (2017).</footnote>

<footnote n="44">⁴⁴L. Freitag, Y. Ma, A. Baiardi, S. Knecht, and M. Reiher, “Approximate analytical gradients and nonadiabatic couplings for the state-average density matrix renormalization group self-consistent-field method,” J. Chem. Theory Comput. **15**, 6724–6737 (2019).</footnote>

<footnote n="45">⁴⁵L. Freitag, S. Knecht, C. Angeli, and M. Reiher, “Multireference perturbation theory with Cholesky decomposition for the density matrix renormalization group,” J. Chem. Theory Comput. **13**, 451–459 (2017).</footnote>
<footnote n="46">⁴⁶P. Sharma, V. Bernales, S. Knecht, D. G. Truhlar, and L. Gagliardi, “Density matrix renormalization group pair-density functional theory (DMRG-PDFT): Singlet-triplet gaps in polyacenes and polyacetylenes,” Chem. Sci. **10**, 1716 (2019).</footnote>

<footnote n="47">⁴⁷S. Knecht, S. Keller, J. Autschbach, and M. Reiher, “A nonorthogonal state-interaction approach for matrix product state wave functions,” J. Chem. Theory Comput. **12**, 5881–5894 (2016).</footnote>

<footnote n="48">⁴⁸S. Dolgov, B. Khoromskij, I. Oseledets, and D. Savostyanov, “Computation of extreme eigenvalues in higher dimensions using block tensor train format,” Comput. Phys. Comm. **185**, 1207–1216 (2014).</footnote>

<footnote n="49">⁴⁹J. Ståring, A. Bernhardsson, and R. Lindh, “Analytical gradients of a state average mcsf state and a state average diagnostic,” Mol. Phys. **99**, 103–114 (2001).</footnote>

<footnote n="50">⁵⁰C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger, and J.-P. Malrieu, “Introduction of n-electron valence states for multireference perturbation theory,” J. Chem. Phys. **114**, 10252 (2001).</footnote>

<footnote n="51">⁵¹C. Angeli, R. Cimiraglia, and J.-P. Malrieu, “N-electron valence state perturbation theory: A fast implementation of the strongly contracted variant,” Chem. Phys. Lett. **350**, 297–305 (2001).</footnote>

<footnote n="52">⁵²C. Angeli, R. Cimiraglia, and J.-P. Malrieu, “N-electron valence state perturbation theory: A spinless formulation and an efficient implementation of the strongly contracted and of the partially contracted variants,” J. Chem. Phys. **117**, 9138–9153 (2002).</footnote>

<footnote n="53">⁵³C. Angeli, S. Borini, M. Cestari, and R. Cimiraglia, “A quasidgenerate formulation of the second order n-electron valence state perturbation theory approach,” J. Chem. Phys. **121**, 4043–4049 (2004).</footnote>

<footnote n="54">⁵⁴T. Husch, L. Freitag, and M. Reiher, “Calculation of ligand dissociation energies in large transition-metal complexes,” J. Chem. Theory Comput. **14**, 2456–2468 (2018).</footnote>

<footnote n="55">⁵⁵S. Wouters, W. Poelmans, P. W. Ayers, and D. Van Neck, “CheMPS2: A free open-source spin-adapted implementation of the density matrix renormalization group for ab initio quantum chemistry,” Comput. Phys. Commun. **185**, 1501–1514 (2014).</footnote>

<footnote n="56">⁵⁶S. Wouters, V. Van Speybroeck, and D. Van Neck, “DMRG-CASPT2 study of the longitudinal static second hyperpolarizability of all-trans polyenes,” J. Chem. Phys. **145**, 054120 (2016).</footnote>

<footnote n="57">⁵⁷T. Yanai, Y. Kurashige, W. Mizukami, J. Chalupský, T. N. Lan, and M. Saitow, “Density matrix renormalization group for ab initio calculations and associated dynamic correlation methods: A review of theory and applications,” Int. J. Quantum Chem. **115**, 283–299 (2015).</footnote>

<footnote n="58">⁵⁸Q. M. Phung, A. Domingo, and K. Pierloot, “Dinuclear Iron(II) spin-crossover compounds: A theoretical study,” Chem. Eur. J. **24**, 5183–5190 (2018).</footnote>

<footnote n="59">⁵⁹K. Pierloot, Q. M. Phung, and A. Domingo, “Spin state energetics in first-row transition metal complexes: contribution of (3s3p) correlation and its description by second-order perturbation theory,” J. Chem. Theory Comput. **13**, 537–553 (2017).</footnote>

<footnote n="60">⁶⁰Q. M. Phung and K. Pierloot, “The dioxygen adducts of iron and manganese porphyrins: electronic structure and binding energy,” Phys. Chem. Chem. Phys. **20**, 17009–17019 (2018).</footnote>

<footnote n="61">⁶¹Q. M. Phung and K. Pierloot, “Low-lying electronic states in chloro-ligated Iron(IV)-Oxo porphyrin as a model for compound I, studied with second-order perturbation theory based on density matrix renormalization group,” J. Chem. Theory Comput. (2019).</footnote>

<footnote n="62">⁶²Q. M. Phung, S. Wouters, and K. Pierloot, “Cumulant approximated second-order perturbation theory based on the density matrix renormalization group for transition metal complexes: a benchmark study,” J. Chem. Theory Comput. **12**, 4352–4361 (2016).</footnote>

<footnote n="63">⁶³L. N. Pham and M. T. Nguyen, “Another look at photoelectron spectra of the anion $Cr_2O_2$: Multireference character and energetic degeneracy,” J. Chem. Theory Comput. **14**, 4833–4843 (2018).</footnote>

<footnote n="64">⁶⁴Q. M. Phung and K. Pierloot, “Electronic structure of N-bridged high-valent Diiron-Oxo,” Chem. Eur. J. **25**, 12491–12496 (2019).</footnote>

<footnote n="65">⁶⁵Y. Kurashige, J. Chalupský, T. N. Lan, and T. Yanai, “Complete active space second-order perturbation theory with cumulant approximation for extended active-space wavefunction from density matrix renormalization group,” J. Chem. Phys. **141**, 174111 (2014).</footnote>

<footnote n="66">⁶⁶B. O. Roos, R. Lindh, P. Å. Malmqvist, V. Veryazov, and P.-O. Widmark, *Multiconfigurational Quantum Chemistry* (Wiley & Sons, Hoboken, 2016) pp. 1–224.</footnote>

<footnote n="67">⁶⁷V. Veryazov, P. Å. Malmqvist, and B. O. Roos, “How to select active space for multiconfigurational quantum chemistry?” Inter. J. Quantum Chem. **111**, 3329–3338 (2011).</footnote>

<footnote n="68">⁶⁸C. Stein and M. Reiher, “SCINE autoCAS,” (2018), https://scine.ethz.ch/download/autocas (accessed Dec 5th 2019).</footnote>

<footnote n="69">⁶⁹C. J. Stein and M. Reiher, “autoCAS: a program for fully automated multi-configurational calculations,” J. Comput. Chem. **40**, 2216–2226 (2019).</footnote>

<footnote n="70">⁷⁰C. J. Stein and M. Reiher, “Automated selection of active orbital spaces,” J. Chem. Theory Comput. **12**, 1760–1771 (2016).</footnote>

<footnote n="71">⁷¹C. J. Stein and M. Reiher, “Automated identification of relevant frontier orbitals for chemical compounds and processes,” Chimia **71**, 170–176 (2017).</footnote>

<footnote n="72">⁷²C. J. Stein, V. von Burg, and M. Reiher, “The delicate balance of static and dynamic electron correlation,” J. Chem. Theory Comput. **12**, 3764–3773 (2016).</footnote>

<footnote n="73">⁷³B. A. Heß, C. M. Marian, U. Wahlgren, and O. Gropen, “A mean-field spin-orbit method applicable to correlated wavefunctions,” Chem. Phys. Letters **251**, 365–371 (1996).</footnote>

<footnote n="74">⁷⁴S. Vancoillie, P.-Å. Malmqvist, and K. Pierloot, “Calculation of EPR g Tensors for Transition-Metal Complexes Based on Multiconfigurational Perturbation Theory (CASPT2),” Chem. Phys. **8**, 1803–1815 (2007).</footnote>

<footnote n="75">⁷⁵B. A. Hess, “Relativistic electronic-structure calculations employing a 2-component no-pair formalism with external-field projection operators,” Phys. Rev. A **33**, 3742–3748 (1986).</footnote>

<footnote n="76">⁷⁶B. O. Roos and P.-Å. Malmqvist, “Relativistic quantum chemistry: the multiconfigurational approach,” Phys. Chem. Chem. Phys. **6**, 2919–2927 (2004).</footnote>

<footnote n="77">⁷⁷B. O. Roos, P.-Å. Malmqvist, and L. Gagliardi, “Exploring the actinide-actinide bond: Theoretical studies of the chemical bond in $Ac_2$, $Th_2$, $Pa_2$, and $U_2$,” J. Am. Chem. Soc. **128**, 17000–17006 (2006).</footnote>

<footnote n="78">⁷⁸J. T. Lyon, L. Andrews, P.-Å. Malmqvist, B. O. Roos, T. Wang, and B. E. Bursten, “Infrared spectrum and bonding in uranium methylidene dihydride, $CH_2=UH_2$,” Inorg. Chem. **46**, 4917–4925 (2007).</footnote>

<footnote n="79">⁷⁹P.-Å. Malmqvist, K. Pierloot, A. R. M. Shahi, C. J. Cramer, and L. Gagliardi, “The restricted active space followed by second-order perturbation theory method: Theory and application to the study of $CuO_2$ and $Cu_2O_2$ systems,” J. Chem. Phys. **128**, 204109 (2008).</footnote>

<footnote n="80">⁸⁰V. Sauri, L. Serrano-Andrés, A. R. M. Shahi, L. Gagliardi, S. Vancoillie, and K. Pierloot, “Multiconfigurational second-order perturbation theory restricted active space (RASPT2) method for electronic excited states: a benchmark study,” J. Chem. Theory Comput. **7**, 153–168 (2011).</footnote>

<footnote n="81">⁸¹S. Vancoillie, P.-Å. Malmqvist, and V. Veryazov, “Potential energy surface of the chromium dimer re-re-revisited with multiconfigurational perturbation theory,” J. Chem. Theory Comput. **12**, 1647–1655 (2016).</footnote>

<footnote n="82">⁸²J. S. Anderson, A. T. Gallagher, J. A. Mason, and T. D. Harris, “A five-coordinate heme dioxygen adduct isolated within a metal–organic framework,” J. Am. Chem. Soc. **136**, 16489–16492 (2014).</footnote>

<footnote n="83">⁸³A. T. Gallagher, M. L. Kelty, J. G. Park, J. S. Anderson, J. A. Mason, J. P. Walsh, S. L. Collins, and T. D. Harris, “Dioxygen binding at a four-coordinate cobaltous porphyrin site in a metal–organic framework: structural, EPR, and $O_2$ adsorption</footnote>
analysis," Inorg. Chem. Front. 3, 536-540 (2016).

⁸⁴G. N. Schrauzer and L.-P. Lee, "The molecular and electronic structure of vitamin B₁₂*, cobaloximes(II), and related compounds," J. Am. Chem. Soc. **90**, 6541-6543 (1968).

⁸⁵J. H. Zagal and M. T. Koper, "Reactivity descriptors for the activity of molecular MN4 catalysts for the oxygen reduction reaction," Angew. Chem. Int. Ed. **55**, 14510-14521 (2016).

⁸⁶J. E. Linard, P. E. Ellis Jr, J. R. Budge, R. D. Jones, and F. Basolo, "Oxygenation of iron (II) and cobalt (II) "capped" porphyrins," J. Am. Chem. Soc. **102**, 1896-1904 (1980).

⁸⁷A. Nyvang, P. A. Malmqvist, and V. Veryazov, "A generally contracted gaussian basis set for photochemistry and photo-physics in the range 0-2 keV," (2017).

⁸⁸M. G. Delcey, L. K. Sørensen, M. Vacher, R. C. Couto, and M. Lundberg, "Efficient calculations of a large number of highly excited states for multiconfigurational wavefunctions," Journal of Computational Chemistry **40**, 1789-1799 (2019).

⁸⁹L. K. Sørensen, M. Guo, R. Lindh, and M. Lundberg, "Applications to metal K pre-edges of transition metal dimers illustrate the approximate origin independence for the intensities in the length representation," Mol. Phys. **115**, 174-189 (2017).

⁹⁰R. J. Cogdell, A. Gall, and J. Köhler, "The architecture and function of the light-harvesting apparatus of purple bacteria: from single molecules to in vivo membranes," Q. Rev. Biophys. **39**, 227-324 (2006).

⁹¹A. Anda, T. Hansen, and L. De Vico, "Multireference Excitation Energies for Bacteriochlorophylls A within Light Harvesting System 2," J. Chem. Theory Comput. **12**, 1305-1313 (2016).

⁹²L. De Vico, A. Anda, V. A. Osipov, A. Ø. Madsen, and T. Hansen, "Macrocycle Ring Deformation as the Secondary Design Principle for Light-Harvesting Complexes," Proc. Natl. Acad. Sci. U. S. A. **115**, E9051-E9057 (2018).

⁹³A. Anda, L. De Vico, and T. Hansen, "Intermolecular Modes between LH2 Bacteriochlorophylls and Protein Residues: The Effect on the Excitation Energies," J. Phys. Chem. B **121**, 5499-5508 (2017).

⁹⁴B. O. Roos, R. Lindh, P.-Å. Malmqvist, V. Veryazov, and P.-O. Widmark, "New relativistic ANO basis sets for transition metal atoms," J. Phys. Chem. A **109**, 6575-6579 (2005).

⁹⁵A. Anda, T. Hansen, and L. De Vico, "Qy and Qx Absorption Bands for Bacteriochlorophyll a Molecules from LH2 and LH3," J. Phys. Chem. A **123**, 5283-5292 (2019).

⁹⁶L. Seijo and Z. Barandiarán, "Ab initio calculations on excited states of lanthanide containing materials," in *Handbook on the Physics and Chemistry of Rare Earths*, Vol. 50, edited by J. C. Bünzli and V. K. Pecharsky (Elsevier, Amsterdam, 2016) Chap. 285, pp. 65-89.

⁹⁷Z. Barandiarán and L. Seijo, "The ab initio model potential representation of the crystalline environment, theoretical study at the local distortion on nacl:cu⁺," J. Chem. Phys. **89**, 5739-5746 (1988).

⁹⁸M. de Jong, D. Biner, K. Krämer, Z. Barandiarán, L. Seijo, and A. Meijerink, "New insights in the 4f¹²⁵d¹ excited states of Tm²⁺ through excited state excitation spectroscopy," J. Phys. Chem. Lett. **7**, 2730-2734 (2016).

⁹⁹R. B. Hughes-Currie, K. V. Ivanovskikh, J. P. R. Wells, M. F. Reid, R. A. Gordon, L. Seijo, and Z. Barandiarán, "X-ray excitation triggers ytterbium anomalous emission in CaF₂:Yb but not in SrF₂:Yb," J. Phys. Chem. Lett. **8**, 1175-1178 (2017).

¹⁰⁰Z. Barandiarán, M. Bettinelli, and L. Seijo, "Color control of Pr³⁺ luminescence by electron-hole recombination energy transfer in CaTiO₃ and CaZrO₃," J. Phys. Chem. Lett. **8**, 3095-3100 (2017).

¹⁰¹J. Joos, L. Seijo, and Z. Barandiarán, "Direct evidence of intervalence charge-transfer states of Eu-doped luminescent materials," J. Phys. Chem. Lett. **10**, 1581-1586 (2019).

¹⁰²L. Seijo, S. P. Feofilov, and Z. Barandiarán, "Fine tuning the Cr³⁺ R₁-line by controlling Pauli antisymmetry strength," J. Phys. Chem. Lett. **10**, 3176-3180 (2019).

¹⁰³A. Abragam and B. Bleaney, *Electron Paramagnetic Resonance of Transition Ions* (Clarendon Press, Oxford, 1970).

¹⁰⁴L. F. Chibotaru and L. Ungur, "Ab initio calculation of anisotropic magnetic properties of complexes. I. Unique definition of pseudospin Hamiltonians and their derivation," Journal of Chemical Physics **137**, 064112 (2012).

¹⁰⁵L. F. Chibotaru and L. Ungur, "Negative *g* factors, berry phases, and magnetic properties of complexes," Phys. Rev. Lett. **109**, 246403 (2012).

¹⁰⁶M. H. L. Pryce, "Sign of *g* in Magnetic Resonance, and the Sign of the Quadrupole Moment of Np²³⁷," Phys. Rev. Lett. **3**, 375-375 (1959).

¹⁰⁷L. Ungur and L. F. Chibotaru, "Ab initio crystal field for lanthanides," Chemistry – A European Journal **23**, 3708–3718 (2017).

¹⁰⁸C. Rudowicz and C. Y. Chung, "The generalization of the extended stevens operators to higher ranks and spins, and a systematic review of the tables of the tensor operators and their matrix elements," Journal of Physics: Condensed Matter **16**, 5825-5847 (2004).

¹⁰⁹N. Iwahara and L. F. Chibotaru, "Exchange interaction between *j* multiplets," Phys. Rev. B **91**, 174438 (2015).

¹¹⁰L. Ungur and L. F. Chibotaru, "Magnetic anisotropy in the excited states of low symmetry lanthanide complexes," Phys. Chem. Chem. Phys. **13**, 20086-20090 (2011).

¹¹¹L. Ungur and L. F. Chibotaru, "Strategies toward high-temperature lanthanide-based single-molecule magnets," Inorganic Chemistry **55**, 10043-10056 (2016).

¹¹²K. L. M. Harriman, D. Errulat, and M. Murugesu, "Magnetic axiality: Design principles from molecules to materials," Trends in Chemistry **1**, 425-439 (2019).

¹¹³L. Ungur, J. J. Le Roy, I. Korobkov, M. Murugesu, and L. F. Chibotaru, "Fine-tuning the local symmetry to attain record blocking temperature and magnetic remanence in a single-ion magnet," Angewandte Chemie International Edition **53**, 4413-4417 (2014).

¹¹⁴S. T. Liddle and J. van Slageren, "Improving *f*-element single molecule magnets," Chem. Soc. Rev. **44**, 6655-6669 (2015).

¹¹⁵C. A. P. Goodwin, F. Ortu, D. Reta, N. F. Chilton, and D. P. Mills, "Molecular magnetic hysteresis at 60 kelvin in dysprosocenium," Nature **548**, 439-442 (2017).

¹¹⁶Y.-C. Chen, J.-L. Liu, W. Wernsdorfer, D. Liu, L. F. Chibotaru, X.-M. Chen, and M.-L. Tong, "Hyperfine-Interaction-Driven Suppression of Quantum Tunneling at Zero Field in a Holmium(III) Single-Ion Magnet," Angewandte Chemie International Edition **56**, 4996-5000 (2017).

¹¹⁷M. Murugesu and R. L. A. Layfield, eds., *Lanthanides and Actinides in Molecular Magnetism* (Wiley-VCH, 2015).

¹¹⁸Y.-C. Chen, J.-L. Liu, L. Ungur, J. Liu, Q.-W. Li, L.-F. Wang, Z.-P. Ni, L. F. Chibotaru, X.-M. Chen, and M.-L. Tong, "Symmetry-Supported Magnetic Blocking at 20 K in Pentagonal Bipyramidal Dy(III) Single-Ion Magnets," Journal of the American Chemical Society **138**, 2829-2837 (2016).

¹¹⁹J. Tang and P. Zhang, *Lanthanide Single Molecule Magnets* (Springer, 2015).

¹²⁰F.-S. Guo, B. M. Day, Y.-C. Chen, M.-L. Tong, A. Mansikkamäki, and R. A. Layfield, "Magnetic hysteresis up to 80 kelvin in a dysprosium metallocene single-molecule magnet," Science **362**, 1400-1403 (2018), https://science.sciencemag.org/content/362/6421/1400.full.pdf.

¹²¹L. Ungur and L. F. Chibotaru, "Computational modelling of the magnetic properties of lanthanide compounds," in *Lanthanides and Actinides in Molecular Magnetism* (John Wiley and Sons, Ltd, 2015) Chap. 6, pp. 153-184.

¹²²L. F. Chibotaru, L. Ungur, and A. Soncini, "The origin of non-magnetic kramers doublets in the ground state of dysprosium triangles: Evidence for a toroidal magnetic moment," Angewandte Chemie International Edition **47**, 4126-4129 (2008).

¹²³M. E. Lines, "Orbital angular momentum in the theory of para-magnetic clusters," The Journal of Chemical Physics **55**, 2977-
2984 (1971).

124K. Meyer, D. J. Mindiola, T. A. Baker, W. M. Davis, and C. C. Cummins, "Uranium hexakisamido complexes," Angew. Chem. Int. Ed. **39**, 3063-3066 (2000).

125J. 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. Journal of Quantum Chem. **114**, 321-327 (2014).

126A. M. Sand, C. E. Hoyer, K. Sharkas, K. M. Kidder, R. Lindh, D. G. Truhlar, and L. Gagliardi, "Analytic gradients for complete active space pair-density functional theory," J. Chem. Theory Comput. **14**, 126-138 (2017).

127F. Aquilante, R. Lindh, and T. B. Pedersen, "Analytic derivatives for the Cholesky representation of the two-electron integrals," J. Chem. Phys. **129**, 034106 (2008).

128J. Boström, F. Aquilante, T. B. Pedersen, and R. Lindh, "Analytical gradients of Hartree-Fock exchange with density fitting approximations," J. Chem. Theory Comput. **9**, 204-212 (2013).

129M. G. Delcey, T. B. Pedersen, F. Aquilante, and R. Lindh, "Analytical gradients of the state-average complete active space self-consistent field method with density fitting," J. Chem. Phys. **143**, 044110 (2015).

130R. Lindh, A. Bernhardsson, and M. Schütz, "Force-constant weighted redundant coordinates in molecular geometry optimizations," Chem. Phys. Letters **303**, 567-575 (1999).

131L. De Vico, M. Olivucci, and R. Lindh, "New General Tools for Constrained Geometry Optimizations," J. Chem. Theory Comput. **1**, 1029-1037 (2005).

132M. Stenrup, R. Lindh, and I. Fdez. Galván, "Constrained numerical gradients and composite gradients: Practical tools for geometry optimization and potential energy surface navigation," J. Comput. Chem. **36**, 1698-1708 (2015).

133F. Häse, I. Fdez. Galván, A. Aspuru-Guzik, R. Lindh, and M. Vacher, "How machine learning can assist the interpretation of ab initio molecular dynamics simulations and conceptual understanding of chemistry," Chem. Sci. **10**, 2298-2307 (2019).

134M. Vacher, A. Brakestad, H. O. Karlsson, I. Fdez. Galván, and R. Lindh, "Dynamical insights into the decomposition of 1,2-dioxetane," Journal of Chemical Theory and Computation **13**, 2448-2457 (2017).

135M. Vacher, P. Farahani, A. Valentini, L. M. Frutos, H. O. Karlsson, I. Fdez. Galván, and R. Lindh, "How do methyl groups enhance the triplet chemiexcitation yield of dioxetane?" The Journal of Physical Chemistry Letters **8**, 3790-3794 (2017).

136A. Valentini, D. Rivero, F. Zapata, C. García-Iriepa, M. Marazzi, R. Palmeiro, I. Fdez. Galván, D. Sampedro, M. Olivucci, and L. M. Frutos, "Optomechanical control of quantum yield in trans-cis ultrafast photoisomerization of a retinal chromophore model," Angew. Chem. Int. Ed. **56**, 3842-3846 (2017).

137M. Manathunga, X. Yang, H. L. Luk, S. Gozem, L. M. Frutos, A. Valentini, N. Ferré, and M. Olivucci, "Probing the Photo-dynamics of Rhodopsins with Reduced Retinal Chromophores," Journal of Chemical Theory and Computation **12**, 839-850 (2016).

138G. Marchand, J. Eng, I. Schapiro, A. Valentini, L. M. Frutos, E. Pieri, M. Olivucci, J. Léonard, and E. Gindensperger, "Directionality of Double-Bond Photoisomerization Dynamics Induced by a Single Stereogenic Center," The Journal of Physical Chemistry Letters **6**, 599-604 (2015).

139C. García-Iriepa, M. Marazzi, F. Zapata, A. Valentini, D. Sampedro, and L. M. Frutos, "Chiral hydrogen bond environment providing unidirectional rotation in photoactive molecular motors," Journal of Physical Chemistry Letters **4**, 1389-1396 (2013).

140C. Schnedermann, X. Yang, M. Liebel, K. M. Spillane, J. Lugtenburg, I. Fernández, A. Valentini, I. Schapiro, M. Olivucci, P. Kukura, and R. A. Mathies, "Evidence for a vibrational phase-dependent isotope effect on the photochemistry of vision," Nature Chemistry **10**, 449455 (2018).

141O. Schalk, J. Galiana, T. Geng, T. L. Larsson, R. D. Thomas, I. Fdez. Galvn, T. Hansson, and M. Vacher, "Competition between ring-puckering and ring-opening excited state reactions exemplified on 5h-furan-2-one and derivatives," The Journal of Chemical Physics **152**, 064301 (2020).

142J. C. Tully, "Molecular dynamics with electronic transitions," J. Chem. Phys. **93**, 1061-1071 (1990).

143N. Ferré and J. G. Ángyán, "Approximate electrostatic interaction operator for qm/mm calculations," Chem. Phys. Lett. **356**, 331339 (2002).

144V. A. Borin, C. Wiebeler, and I. Schapiro, "A qm/mm study of the initial excited state dynamics of green-absorbing proteins," Faraday Discuss. **207**, 137-152 (2018).

145M. Richter, P. Marquetand, J. González-Vázquez, I. Sola, and L. González, "SHARC: Ab initio molecular dynamics with surface hopping in the adiabatic representation including arbitrary couplings," J. Chem. Theory Comput. **7**, 1253-1258 (2011).

146J. Norell, M. Odelius, and M. Vacher, "Ultrafast dynamics of photo-excited 2-thiopyridone: Theoretical insights into triplet state population and proton transfer pathways," Structural Dynamics **7**, 024101 (2020).

147O. Weingart, A. Nenov, P. Altoè, I. Rivalta, J. Segarra-Martí, I. Dokukina, and M. Garavelli, "COBRAMM 2.0 — a software interface for tailoring molecular electronic structure calculations and running nanoscale (QM/MM) simulations," Journal of Molecular Modeling **24**, 271 (2018).

148G. A. Worth, K. Giri, G. W. Richings, I. Burghardt, M. H. Beck, A. Jäckle, and H.-D. Meyer, "The QUANTICS Package, Version 1.1, University of Birmingham, Birmingham, U.K." (2015).

149G. W. Richings, I. Polyak, K. E. Spinlove, G. A. Worth, I. Burghardt, and B. Lasorne, "Quantum dynamics simulations using gaussian wavepackets: the vMCG method," Int. Rev. Phys. Chem. **34**, 269-308 (2015).

150F. Melaccio, M. d. C. Marín, A. Valentini, F. Montisci, S. Rinaldi, M. Cherubini, X. Yang, Y. Kato, M. Stenrup, Y. Orozco-González, N. Ferré, H. L. Luk, H. Kandori, and M. Olivucci, "Toward Automatic Rhodopsin Modeling as a Tool for High-Throughput Computational Photobiology," J. Chem. Theory Comput. **12**, 6020-6034 (2016).

151L. Pedraza-González, L. De Vico, M. d. C. Marín, F. Fanelli, and M. Olivucci, "a-ARM: Automatic Rhodopsin Modeling with Chromophore Cavity Generation, Ionization State Selection and External Counter-ion Placement," J. Chem. Theory Comput. **15**, 3134-3152 (2019).

152D.-E. Nilsson, "Photoreceptor evolution: Ancient siblings serve different tasks," Curr. Biol. **15**, R94-R96 (2005).

153O. P. Ernst, D. T. Lodowski, M. Elstner, P. Hegemann, L. S. Brown, and H. Kandori, "Microbial and animal rhodopsins: Structures, functions, and molecular mechanisms," Chem. Rev. **114**, 126-163 (2013).

154I. Gushchin, V. Shevchenko, V. Polovinkin, V. Borshchevskiy, P. Buslaev, E. Bamberg, and V. Gordeliy, "Structure of the light-driven sodium pump KR2 and its implications for optogenetics," The FEBS J. **283**, 1232-1238 (2016).

155D. Agathangelou, Y. Orozco-Gonzalez, M. d. C. Marín, P. Roy, J. Brazard, H. Kandori, K.-H. Jung, J. Léonard, T. Buckup, N. Ferré, M. Olivucci, and S. Haacke, "Effect of point mutations on the ultrafast photo-isomerization of *Anabaena* sensory rhodopsin," Faraday Discuss. **207**, 55-75 (2018).

156M. d. C. Marín, D. Agathangelou, Y. Orozco-Gonzalez, A. Valentini, Y. Kato, R. Abe-Yoshizumi, H. Kandori, A. Choi, K.-H. Jung, S. Haacke, and M. Olivucci, "Fluorescence enhancement of a microbial rhodopsin via electronic reprogramming," J. Am. Chem. Soc. **141**, 262-271 (2019).

157K. Inoue, M. d. C. Marín, S. Tomida, R. Nakamura, Y. Nakajima, M. Olivucci, and H. Kandori, "Red-shifting mutation of light-driven sodium-pump rhodopsin," Nat. Commun. **10**, 1993 (2019).

158S. Gholami, L. Pedraza-González, X. Yang, A. A. Granovsky, I. N. Ioffe, and M. Olivucci, "Multi-state multi-configuration
quantum chemical computation of the two-photon absorption spectra of bovine rhodopsin," J. Phys. Chem. Lett. **10**, 6293-6300 (2019).

¹⁵⁹L. Zhang and J. Hermans, "Hydrophilicity of cavities in proteins," Proteins: Structure, Function, and Bioinformatics **24**, 433–438 (1996).

¹⁶⁰S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apos-
tolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. Van Der Spoel,
B. Hess, and E. Lindahl, "Gromacs 4.5: A high-throughput and
highly parallel open source molecular simulation toolkit," Bioin-
formatics **29**, 845–854 (2013).

¹⁶¹T. Andruniów, N. Ferré, and M. Olivucci, "Structure, initial excited-state relaxation, and energy storage of rhodopsin resolved at the multi-configurational perturbation theory level," Proc. Natl. Acad. Sci. U. S. A. **101**, 17908–17913 (2004).

¹⁶²J. A. Rackers, Z. Wang, C. Lu, M. L. Laury, L. Lagardére, M. J. Schnieders, J.-P. Piquemal, P. Ren, and J. W. Ponder, "Tinker 8: Software tools for molecular design," J. Chem. Theory Comput. **14**, 5273–5289 (2018).

¹⁶³E. Pieri, V. Ledentu, M. Sahlin, F. Dehez, M. Olivucci, and N. Ferré, “CpHMD-Then-QM/MM Identification of the Amino Acids Responsible for the Anabaena Sensory Rhodopsin pH-Dependent Electronic Absorption Spectrum,” J. Chem. Theory Comput. **15**, 4535–4546 (2019).

¹⁶⁴L. Pedraza-González, M. D. C. Marín, A. N. Jorge, T. D. Ruck, X. Yang, A. Valentini, M. Olivucci, and L. De Vico, “Web-ARM: a Web-Based Interface for the Automatic Construction of QM/MM Models of Rhodopsins,” J. Chem. Inf. Model. **60**, 1481–1493 (2020).

¹⁶⁵F. Segatta, L. Cupellini, S. Jurinovich, S. Mukamel, M. Dapor, S. Taioli, M. Garavelli, and B. Mennucci, “A Quantum Chemical Interpretation of Two-Dimensional Electronic Spectroscopy of Light-Harvesting Complexes,” J. Am. Chem. Soc. **139**, 7558–7567 (2017).

¹⁶⁶A. Nenov, R. Borrego-Varillas, A. Oriana, L. Ganzer, F. Segatta, I. Conti, J. Segarra-Martí, J. Omachi, M. Dapor, S. Taioli, C. Manzoni, S. Mukamel, G. Cerullo, and M. Garavelli, “Uv-light-induced vibrational coherences: The key to understand kasha rule violation in trans -azobenzene,” Journal of Physical Chemistry Letters **9**, 1534–1541 (2018).

¹⁶⁷F. Segatta, A. Nenov, S. Orlandi, A. Arcioni, S. Mukamel, and M. Garavelli, "Exploring the capabilities of optical pump x-ray probe NEXAFS spectroscopy to track photo-induced dynamics mediated by conical intersections," Faraday Discussions (2020), 10.1039/c9fd00073a.

¹⁶⁸R. Borrego-Varillas, D. Teles-Ferreira, A. Nenov, I. Conti, L. Ganzer, C. Manzoni, M. Garavelli, A. Maria De Paula, and G. Cerullo, "Observation of the sub-100 femtosecond population of a dark state in a thiobase mediating intersystem crossing," Journal of the American Chemical Society **140**, 16087–16093 (2018).

¹⁶⁹A. Picchiotti, A. Nenov, A. Giussani, V. I. Prokhorenko, R. J. D. Miller, S. Mukamel, and M. Garavelli, "Pyrene, a test case for deep-ultraviolet molecular photophysics," Journal of Physical Chemistry Letters **10**, 3481–3487 (2019).

¹⁷⁰D. C. Teles-Ferreira, I. Conti, R. Borrego-Varillas, N. A., L. Van Stokkum, I. H. M. adn Ganzer, C. Manzoni, A. M. de Paula, G. Cerullo, and M. Garavelli, “A unified experimental/theoretical description of the ultrafast photophysics of single and double thionated uracils,” Chemistry – A European Journal **26**, 336–343 (2020).

¹⁷¹A. Nenov, F. Segatta, A. Bruner, S. Mukamel, and M. Garavelli, “X-ray linear and non-linear spectroscopy of the esca molecule,” Journal of Chemical Physics **151** (2019), 10.1063/1.5116699.

¹⁷²J. Segarra-Martí, S. Mukamel, M. Garavelli, A. Nenov, and I. Rivalta, “Towards accurate simulation of two-dimensional electronic spectroscopy,” Topics in Current Chemistry **376** (2018).

¹⁷³A. Nenov, I. Conti, R. Borrego-Varillas, G. Cerullo, and M. Gar-
avelli, “Linear absorption spectra of solvated thiouracils resolved

at the hybrid RASPT2/MM level,” Chemical Physics **515**, 643–653 (2018).

¹⁷⁴D. Abramavicius, B. Palmieri, D. Voronine, F. Šanda, and S. Mukamel, “Coherent multidimensional optical spectroscopy of excitons in molecular aggregates; quasiparticle versus supermolecule perspectives,” Chemical Reviews **109**, 2350–2408 (2009).

¹⁷⁵R. Borrego-Varillas, A. Nenov, L. Ganzer, A. Oriana, C. Man-
zoni, A. Tolomelli, I. Rivalta, S. Mukamel, M. Garavelli, and
G. Cerullo, "Two-dimensional uv spectroscopy: A new insight
into the structure and dynamics of biomolecules," Chemical Sci-
ence **10**, 9907–9921 (2019).

¹⁷⁶T. Wolf, R. Myhre, J. Cryan, S. Coriani, R. Squibb, A. Battistoni, N. Berrah, C. Bostedt, P. Bucksbaum, G. Coslovich, R. Feifel, K. Gaffney, J. Grilj, T. Martinez, S. Miyabe, S. Moeller, M. Mucke, A. Natan, R. Obaid, T. Osipov, O. Plekan, S. Wang, H. Koch, and M. Ghr, “Probing ultrafast ππ\*-nπ\* internal conversion in organic chromophores via k-edge resonant absorption,” Nature Communications **8** (2017), 10.1038/s41467-017-00069-7.

¹⁷⁷J. Norell, R. M. Jay, M. Hantschmann, S. Eckert, M. Guo, K. J. Gaffney, P. Wernet, M. Lundberg, A. Föhlisch, and M. Odelius, “Fingerprints of electronic, spin and structural dynamics from resonant inelastic soft X-ray scattering in transient photo-chemical species,” Phys. Chem. Chem. Phys. **20**, 7243–7253 (2018).

¹⁷⁸I. Josefsson, K. Kunnus, S. Schreck, A. Föhlisch, F. de Groot, P. Wernet, and M. Odelius, “Ab initio calculations of X-ray spectra: Atomic multiplet and molecular orbital effects in a multiconfigurational SCF approach to the L-edge spectra of transition metal complexes,” J. Phys. Chem. Lett. **3**, 3565–3570 (2012).

¹⁷⁹S. I. Bokarev, M. Dantz, E. Suljoti, O. Kühn, and E. F. Aziz, “State-dependent electron delocalization dynamics at the solute-solvent interface: Soft-X-ray absorption spectroscopy and ab initio calculations,” Phys. Rev. Lett. **111**, 083002–083007 (2013).

¹⁸⁰R. V. Pinjari, M. G. Delcey, M. Guo, M. Odelius, and M. Lundberg, “Restricted active space calculations of L-edge X-ray absorption spectra: From molecular orbitals to multiplet states,” J. Chem. Phys. **141**, 124116 (2014).

¹⁸¹P. Wernet, K. Kunnus, I. Josefsson, I. Rajkovic, W. Quevedo, M. Beye, S. Schreck, S. Gruebel, M. Scholz, D. Nordlund, W. Zhang, R. W. Hartsock, W. F. Schlotter, J. J. Turner, B. Kennedy, F. Hennies, F. M. F. de Groot, K. J. Gaffney, S. Techert, M. Odelius, and A. Föhlisch, “Orbital-specific mapping of the ligand exchange dynamics of Fe(CO)₅ in solution,” Nature **520**, 78–81 (2015).

¹⁸²M. Guo, E. Källman, L. K. Sørensen, M. G. Delcey, R. V. Pin-
jari, and M. Lundberg, “Molecular orbital simulations of metal
1s2p resonant inelastic x-ray scattering,” J. Phys. Chem. A **120**,
5848–5855 (2016).

¹⁸³F. Weinhold, "Natural bond orbital methods," in *Encyclopedia of computational chemistry*, edited by P. von Ragué Schleyer (John Wiley & Sons, Chichester, 1998) pp. 1792–1811.

¹⁸⁴E. D. Glendening, C. R. Landis, and F. Weinhold, “NBO 6.0: Natural bond orbital analysis program,” J. Comput. Chem. **34**, 1429–1437 (2013).

¹⁸⁵F. Gendron, D. Páez-Hernández, F.-P. Notter, B. Pritchard, H. Bolvin, and J. Autschbach, “Magnetic properties and electronic structure of neptunyl<sup>VI</sup> complexes: Wavefunctions, orbitals, and crystal-field models,” Chem. Eur. J. **20**, 7994–8011 (2014).

¹⁸⁶J. Autschbach, “Orbitals for analyzing bonding and magnetism of heavy-metal complexes,” Comments Inorg. Chem. **36**, 215–244 (2016).

¹⁸⁷F. Gendron, B. Le Guennic, and J. Autschbach, “Magnetic properties and electronic structures of Ar₃U<sup>IV</sup>-L complexes with Ar = C₅(CH₃)₄H⁻ or C₅H₅⁻ and L = CH₃, NO, and Cl,” Inorg. Chem. **53**, 13174–13187 (2014).
<footnote n="188">¹⁸⁸T. J. Duignan and A. Marchenko, “Exatomic: A unified platform for computational chemists,” https://github.com/exa-analytics/exatomic, doi:10.5281/zenodo.60053.</footnote>

<footnote n="189">¹⁸⁹D.-C. Sergentu, J. T. Duignan, and J. Autschbach, “Ab initio study of covalency in the ground versus core-excited states and X-ray absorption spectra of actinide complexes,” J. Phys. Chem. Lett. **9**, 5583–5591 (2018).</footnote>

<footnote n="190">¹⁹⁰G. Ganguly, D.-C. Sergentu, and J. Autschbach, “Ab initio analysis of metal-ligand bonding in An(COT)₂, An = Th, U, in their ground- and core-excited states,” Chem. Eur. J. (2020), 10.1002/chem.201904166.</footnote>

<footnote n="191">¹⁹¹P. Wernet, K. Kunnus, S. Schreck, W. Quevedo, R. Kurian, S. Techert, F. M. F. de Groot, M. Odelius, and A. Föhlisch, “Dissecting local atomic and intermolecular interactions of transition-metal ions in solution with selective X-ray spectroscopy,” J. Phys. Chem. Lett. **3**, 3448–3453 (2012).</footnote>

<footnote n="192">¹⁹²R. V. Pinjari, M. G. Delcey, M. Guo, M. Odelius, and M. Lundberg, “Cost and sensitivity of restricted active-space calculations of metal L-edge X-ray absorption spectra,” J. Comput. Chem. **37**, 477–486 (2016).</footnote>

<footnote n="193">¹⁹³M. Preuße, S. I. Bokarev, S. G. Aziz, and O. Kühn, “Towards an ab initio theory for metal L-edge soft X-ray spectroscopy of molecular aggregates,” Struct. Dyn. **3**, 062601 (2016).</footnote>

<footnote n="194">¹⁹⁴T. Vitova, I. Pidchenko, D. Fellhauer, P. S. Bagus, Y. Joly, T. Pruessmann, S. Bahl, E. Gonzalez-Robles, J. Rothe, M. Altmaier, M. A. Denecke, and H. Geckeis, “The role of the 5f valence orbitals of early actinides in chemical bonding,” Nat. Commun. **8**, 16053 (2017).</footnote>

<footnote n="195">¹⁹⁵J. Segarra-Martí, M. Garavelli, and F. Aquilante, “Converging many-body correlation energies by means of sequence extrapolation,” J. Chem. Phys. **148**, 034107 (2018).</footnote>

<footnote n="196">¹⁹⁶J. Segarra-Martí, M. Garavelli, and F. Aquilante, “Multiconfigurational second-order perturbation theory with frozen natural orbitals extended to the treatment of photochemical problems,” J. Chem. Theory Comput. **11**, 3772–3784 (2015).</footnote>

<footnote n="197">¹⁹⁷A. Zech, F. Aquilante, and T. A. Wesolowski, “Orthogonality of embedded wave functions for different states in frozen-density embedding theory,” J. Chem. Phys. **143**, 164106 (2015).</footnote>

<footnote n="198">¹⁹⁸P. Farahani, D. Roca-Sanjuán, and F. Aquilante, “A two-scale approach to electron correlation in multiconfigurational perturbation theory,” J. Comput. Chem. **35**, 1609–1617 (2014).</footnote>

<footnote n="199">¹⁹⁹F. Aquilante, T. K. Todorova, L. Gagliardi, T. B. Pedersen, and B. O. Roos, “Systematic truncation of the virtual space in multiconfigurational perturbation theory,” J. Chem. Phys. **131**, 034113 (2009).</footnote>

<footnote n="200">²⁰⁰M. Nest, T. Klamroth, and P. Saalfrank, “The multiconfiguration time-dependent HartreeFock method for quantum chemical calculations,” J. Chem. Phys. **122**, 124102 (2005).</footnote>

<footnote n="201">²⁰¹T. Sato and K. L. Ishikawa, “Time-dependent complete-active-space self-consistent-field method for multielectron dynamics in intense laser fields,” Phys. Rev. A **88**, 023402 (2013).</footnote>

<footnote n="202">²⁰²H. Miyagi and L. B. Madsen, “Time-dependent restricted-active-space self-consistent-field theory for laser-driven many-electron dynamics. II. Extended formulation and numerical analysis,” Phys. Rev. A **89**, 063416 (2014).</footnote>

<footnote n="203">²⁰³A. Baiardi and M. Reiher, “Large-scale quantum-dynamics with matrix product states,” J. Chem. Theory Comput. **15**, 3481–3498 (2019).</footnote>

<footnote n="204">²⁰⁴E. D. Hedegård, S. Knecht, J. S. Kielberg, H. J. A. Jensen, and M. Reiher, “Density matrix renormalization group with efficient dynamical electron correlation through range separation,” J. Chem. Phys. **142**, 224108 (2015).</footnote>