# Time evolution of natural orbitals in ab initio molecular dynamics

Alejandro Rivero Santamaría, Mario Piris

► To cite this version:

Alejandro Rivero Santamaría, Mario Piris. Time evolution of natural orbitals in ab initio molecular dynamics. The Journal of Chemical Physics, 2024, 160 (7), 10.1063/5.0188491. hal-04491144

HAL Id: hal-04491144

https://hal.science/hal-04491144v1

Submitted on 8 Mar 2024

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Time evolution of natural orbitals in *ab initio* molecular dynamics

Alejandro Rivero Santamaría*

Univ. Lille, CNRS, UMR 8523 – PhLAM – Physique
des Lasers Atomes et Molécules, F-59000 Lille, France

Mario Piris†

Donostia International Physics Center (DIPC), 20018 Donostia, Euskadi,
Spain; Euskal Herriko Unibertsitatea (UPV/EHU), PK 1072, 20080 Donostia, Euskadi,
Spain; and Basque Foundation for Science (IKERBASQUE), 48009 Bilbao, Euskadi, Spain

(Dated: January 10, 2024)
Abstract

This work combines for the first time *ab initio* molecular dynamics (AIMD) within the Born-Oppenheimer approximation, with a global natural orbital functional (GNOF), an approximate functional of the one-particle reduced density matrix. The most prominent feature of GNOF-AIMD is the ability to display the real-time evolution of natural orbitals, providing detailed information on the time-dependent electronic structure of complex systems and processes, including reactive collisions. The quartet ground-state reaction $N(^4S) + H_2(^1Σ) → NH(^3Σ) + H(^2S)$ is taken as validation test. Collision energy influences on integral cross sections for different initial ro-vibrational states of $H_2$ and rotational-state distributions of NH product are discussed, showing a good agreement with previous high-quality theoretical results.

One-particle reduced density matrix (1RDM) functional theory[1-3] is an alternative formalism to both density functional and wavefunction based methods. A pragmatic approach results in approximate functionals of the 1RDM in its diagonal form, that is, the use of natural orbitals (NOs) and its occupation numbers (ONs) as the fundamental variables, which define a natural orbital functional (NOF).[4] An important issue is that the approximate NOFs continue to depend on the two-particle RDM (2RDM),[5] so it is necessary to consider their functional N-representability.[6, 7] The NOF theory is currently an active research field,[8-31] which can already be applied to large molecular systems of general chemical interest[32, 33] using open-source software like DoNOF.[34, 35] Exhaustive reviews of NOF-based methods can be found elsewhere.[36-38]

Knowledge of the analytical derivatives with respect to nuclear coordinates in NOF theory[39, 40] allows routine calculation of molecular structures and related properties. This ability to calculate gradients analytically also makes NOF-based *ab initio* molecular dynamics (AIMD)[41] simulations feasible. The basic idea is to calculate the forces acting on the nuclei from the calculations of the electronic structure as the molecular dynamics trajectory is generated. In this way we avoid determining in advance the potentials of interatomic interactions that can have serious drawbacks in chemically complex systems where the bond pattern changes qualitatively in the course of dynamics.

In many cases, the electrons respond almost instantaneously to nuclear motion. In such situations, the Born-Oppenheimer (BO) approximation allows us to decouple the electronic and nu-

<footnote>* alejandro.rivero@univ-lille.fr</footnote>

<footnote>† mario.piris@ehu.eus</footnote>
clear problems, defining a potential energy surface (PES). In this work, we focus on the gas-phase
dynamics of molecules in their electronic ground state (GS), where the BO approximation is gen-
erally accurate. There are other problems, for example, when the dynamics starts in an excited
electronic state, where strong couplings between two or more PESs occur, requiring quantum
treatment of both nuclei and electrons. Non-adiabatic molecular dynamics[42, 43] is the method
of choice to model these processes, but that is outside the scope of this work.

Nowadays, the most reliable approaches to deal with electron correlation in AIMD are multiref-
erence methods, such as the complete active space self-consistent field (CASSCF)[44, 45] and its
combination with second-order perturbation theory (CASPT2).[46–49] Unfortunately, multirefer-
ence methods can suffer from high computational cost, sensitivity to active space selection, and
instabilities in calculations due to changes in active space orbitals that can lead to energy disconti-
nuities along the trajectory. Consequently, few simulations using CASPT2 nuclear gradients have
been reported to perform non-adiabatic AIMD for reduced molecular models [50, 51] and, so far,
there are no examples of convergent CASPT2-AIMD calculations in reactive collisions.

NOFs offer an intermediate cost between multireference methods and common density func-
tionals. In fact, approximate NOFs have been shown to be more accurate than their electron
density-dependent counterparts and to have better scaling with respect to the number of basic func-
tions than wavefunction-based methods for systems with a large amount of strong non-dynamic
correlation. Interestingly, NOFs further corroborate the motivation behind the floating occupa-
tion molecular orbital complete active space configuration interaction (FOMO-CASCI) method, a
promising alternative to CASSCF in dynamics simulations.[52, 53]

The aim of this article is to present for the first time a BO AIMD based on an approximate NOF,
that is, the nuclei will propagate according to the classical equations of motion, in an adiabatic PES
obtained by solving in each time step the quantum mechanical electronic structure problem using
a NOF. The price to pay is that the correlation lengths and relaxation times that are accessible are
smaller than those available through standard molecular dynamics, but we can handle chemically
complex systems, avoid the dimensionality bottleneck that arises when PESs are calculated in
advance,[41] and see the real-time evolution of NOs and their ONs during complex dynamics,
obtaining detailed information about the time-dependent electronic structure of such processes.

As a functional, we will use the recently proposed global NOF (GNOF)[54] for electronic
systems with any value of spin regardless of external potential. It has been demonstrated[55]
that GNOF provides a good balance between static and dynamic electronic correlations leading to
accurate total energies while preserving spin, even for systems with a highly multi-configurational character. In addition, GNOF has proven to be successful in dealing with aromaticity[56] and charge delocalization error.[57] GNOF correlates all electrons into all available orbitals for a given basis set; which today is not possible for large systems with current wavefunction-based methods. Regarding the dynamics, the NOs vary along the trajectory, adapting at each time step to the most favorable interactions of the corresponding nuclei configuration. In the following, we briefly describe GNOF.

We consider a mixed quantum state (multiplet) of an N-electron molecule with N<sub>I</sub> spin-unpaired electrons, N<sub>II</sub> = N - N<sub>I</sub> spin-paired, and total spin S. We focus on the state of highest multiplicity: 2S + 1 = N<sub>I</sub> + 1.[58] For the whole ensemble, the expected value of Ŝ<sub>z</sub> is zero; so the spin-restricted theory is adopted, i.e., a single set of orbitals is used for α and β spins: φ<sub>p</sub><sup>α</sup>(<b>r</b>) = φ<sub>p</sub><sup>β</sup>(<b>r</b>) = φ<sub>p</sub>(<b>r</b>). Accordingly, all spatial orbitals are double occupied in the ensemble, so that occupancies for particles with α and β spins are equal: n<sub>p</sub><sup>α</sup> = n<sub>p</sub><sup>β</sup> = n<sub>p</sub>.

The orbital space is divided in turn into two subspaces: Ω = Ω<sub>I</sub> ⊕ Ω<sub>II</sub>. Ω<sub>II</sub> is composed of N<sub>II</sub>/2 mutually disjoint subspaces Ω<sub>g</sub>, each of which contains one orbital |g⟩ with g ≤ N<sub>II</sub>/2, and N<sub>g</sub> orbitals |p⟩ with p > N<sub>II</sub>/2, namely,

$$
\Omega_g = \{|g\rangle, |p_1\rangle, |p_2\rangle, \dots, |p_{N_g}\rangle\}, \quad \Omega_g \in \Omega_{II}. \qquad (1)
$$

Taking into account the spin, the total occupancy for a given subspace Ωg is 2, which is reflected
in the following sum rule:

$$
\sum_{p \in \Omega_g} n_p = n_g + \sum_{i=1}^{N_g} n_{p_i} = 1, \quad g = 1, 2, \dots, \frac{N_{II}}{2}. \qquad (2)
$$

In general, N<sub>g</sub> can be different for each subspace as long as it describes the electron pair well. For convenience, we usually take it the same for all subspaces Ω<sub>g</sub> ∈ Ω<sub>II</sub>. The maximum possible value of N<sub>g</sub> is determined by the number of basis functions (N<sub>B</sub>). From (2), it follows that

$$
2 \sum_{p \in \Omega_{II}} n_p = 2 \sum_{g=1}^{N_{II}/2} \left( n_g + \sum_{i=1}^{N_g} n_{p_i} \right) = N_{II}. \qquad (3)
$$

On the other hand, Ω<sub>I</sub> is composed of N<sub>I</sub> mutually disjoint subspaces, each of which contains only one orbital |g⟩. This spatial orbital is individually occupied, n<sub>g</sub><sup>α</sup> + n<sub>g</sub><sup>β</sup> = 2n<sub>g</sub> = 1, but we do not know whether the electron has α or β spin. This leads to

$$
2 \sum_{p \in \Omega_I} n_p = 2 \sum_{g=N_{II}/2+1}^{N_{\Omega}} n_g = N_I. \qquad (4)
$$
In Eq. (4), $N_Ω = N_{II}/2 + N_I$ denotes the total number of suspaces in $Ω$. Note that the orbitals are arranged as follows: from 1 to $N_{II}/2$, there are the strongly occupied orbitals belonging to $Ω_{II}$, followed by the single-occupied orbitals belonging to $Ω_I$ ranging from $N_{II}/2 + 1$ to $N_Ω$. After that, there are the weakly occupied orbitals of each $Ω_g \in Ω_{II}$. Taking into account Eqs. (3) and (4), the trace of the 1RDM is verified to be equal to the number of electrons:

$$2 \sum_{p \in \Omega} n_p = 2 \sum_{p \in \Omega_{II}} n_p + 2 \sum_{p \in \Omega_I} n_p = N_{II} + N_I = N. \quad (5)$$

It is essential to note that orbitals undergo changes throughout the optimization process to find the most favorable orbital interactions. As a result, the orbitals are not static during the optimization process; they adapt to the specific problem.

A reconstruction functional for the 2RDM in terms of the ONs leads to the following electronic energy (GNOF):

$$E_{el} = E^{intra} + E_{HF}^{inter} + E_{sta}^{inter} + E_{dyn}^{inter} \quad (6)$$

The intra-pair component is formed by summing the energies $E_g$ of electron pairs with opposite spins and the one-electron energies of unpaired electrons, specifically,

$$E^{intra} = \sum_{g=1}^{N_{II}/2} E_g + \sum_{g=N_{II}/2+1}^{N_Ω} H_{gg} \quad (7)$$

$$E_g = 2 \sum_{p \in \Omega_g} n_p H_{pp} + \sum_{q,p \in \Omega_g} \Pi(n_q, n_p) L_{pq} \quad (8)$$

where $H_{pp}$ are the diagonal one-electron matrix elements of the kinetic energy and external potential operators, whereas $L_{pq} = \langle pp|qq \rangle$ are the exchange-time-inversion integrals.[59] The matrix elements $\Pi(n_q, n_p) = c(n_q)c(n_p)$, where $c(n_p)$ is defined by the square root of the ONs according to the following rule:

$$c(n_p) = \begin{cases} \sqrt{n_p}, & p \le \frac{N_{II}}{2}, \\ -\sqrt{n_p}, & p > \frac{N_{II}}{2} \end{cases}, \quad p \in \Omega_g \in \Omega_{II} \quad (9)$$

that is, the phase factor of $c_p$ is chosen to be +1 for the strongly occupied orbital of a given subspace $Ω_g$, and -1 otherwise. The inter-subspace HF term is

$$E_{HF}^{inter} = \sum_{p,q=1}^{N_B} n_q n_p (2J_{pq} - K_{pq}) \quad (10)$$
where $J_{pq} = \langle pq | pq \rangle$ and $K_{pq} = \langle pq | qp \rangle$ are the Coulomb and exchange integrals, respectively. The prime in the summation indicates that only the inter-subspace terms are taking into account. The inter-subspace static component is written as

$$ E_{sta}^{inter} = - \left( \sum_{p=1}^{N_{\Omega}} \sum_{q=N_{\Omega}+1}^{N_{B}} + \sum_{p=N_{\Omega}+1}^{N_{B}} \sum_{q=1}^{N_{\Omega}} + \sum_{p,q=N_{\Omega}+1}^{N_{B}} \right)' \Phi_{q} \Phi_{p} L_{pq} - \frac{1}{2} \left( \sum_{p=1}^{N_{II}/2} \sum_{q=N_{II}/2+1}^{N_{\Omega}} + \sum_{p=N_{II}/2+1}^{N_{\Omega}} \sum_{q=1}^{N_{II}/2} \right)' \Phi_{q} \Phi_{p} L_{pq} - \frac{1}{4} \sum_{p,q=N_{II}/2+1}^{N_{\Omega}} K_{pq} \quad (11) $$

where $\Phi_p = \sqrt{n_p h_p}$ with the hole $h_p = 1 - n_p$. Finally, the inter-subspace dynamic energy is

$$ E_{dyn}^{inter} = \sum_{p,q=1}^{N_B} [n_q^d n_p^d + \Pi(n_q^d, n_p^d)] (1 - \delta_{q\Omega_{II}^b} \delta_{p\Omega_{II}^b}) L_{pq} \quad (12) $$

In Eq. (12), $\Omega_{II}^b$ denotes the subspace composed of orbitals below the level $N_{II}/2$, and the dynamic part of the ON $n_p$ is defined as

$$ n_p^d = n_p \cdot e^{-\left(\frac{h_g}{h_c}\right)^2}, \quad p \in \Omega_g, \quad g = 1, 2, \dots, \frac{N_{II}}{2} \quad (13) $$

with $h_c = 0.02\sqrt{2}$. $n_p^d$ is in accordance with the Pulay's criterion that establishes an occupancy deviation of approximately 0.01 with respect to 1 or 0 for a NO to contribute to the dynamic correlation.

In the BO approximation, the total energy of the molecule can be cast as $E = E_{nuc} + E_{el}$ with $E_{nuc} = \sum_{A<B} [Z_A Z_B / R_{AB}]$. $Z_A$ represents the atomic number of nucleus A, and $R_{AB}$ is the distance between nuclei A and B. Considering that all NOs $\phi_i$ are expanded in a fixed atomic basis set, $\phi_i = \sum_v C_{vi} \zeta_v$, the derivative of the total energy with respect to the coordinate $x$ of nucleus A is given by [40]

$$ \frac{dE}{dx_A} = \frac{\partial E_{nuc}}{\partial x_A} + \sum_{\mu\nu} \Gamma_{\mu\nu} \frac{\partial H_{\mu\nu}}{\partial x_A} + \sum_{\mu\eta\nu\delta} D_{\mu\eta\nu\delta} \frac{\partial \langle \mu\eta | v\delta \rangle}{\partial x_A} - \sum_{\mu\nu} \lambda_{\mu\nu} \frac{\partial S_{\mu\nu}}{\partial x_A}. \quad (14) $$

where $\Gamma_{\mu\nu}$ and $D_{\mu\eta\nu\delta}$ are the 1RDM and 2RDM, respectively, $S_{\mu\nu} = \langle \mu | v \rangle$ is the overlap matrix, and $\lambda_{\mu\nu}$ are the Lagrange multipliers obtained from RDMs, all in the atomic orbital representation. The first term of Eq. (14) is the derivative of the nuclear energy, the second represents the negative Hellmann-Feynman force, while the third contains the explicit derivatives of two-electron integrals. The last term, known as the density force, arises from the implicit dependence of $C_{vi}$ on
<figure><img src="image_2.png" /><figcaption>FIG. 2. Time evolution of the two strongly occupied natural orbitals involved in the bond pattern change during the collision.</figcaption></figure>

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

geometry. The implicit dependence of ONs on geometry does not contribute to analytic gradients since $E_{el}$ is stationary with respect to variations in all of the ONs.[39]

All derivatives in Eq. (14) have an explicit dependence on the nuclear coordinate $x_A$, so the force acting on each nucleus A ($F_A = -\nabla_A E$) can be obtained by a single static evaluation at each time step for the fixed nuclear positions at that instant. Consequently, we can calculate the
trajectories of the nuclei according to the classical equations of motion, but taking into account
the quantization of the reactants, a procedure known as quasiclassical trajectory (QCT) method.
It is worth noting that the QCT method does not take the tunneling effect into account, so it can
produce inaccurate results near the threshold energy.

NOF-based QCT calculations can be performed using the new molecular dynamics module im-
plemented in DoNOF,[34] which allows the calculation of nuclear trajectories by determining “on
the fly” the forces using NOF gradients (14). Beeman’s algorithm[60] is used to numerically inte-
grate Newton’s equations of motion, whereas the initial conditions are obtained using a standard
Monte Carlo sampling procedure.[61]

The N(^4S) + H₂(^1Σ) → NH(^3Σ) + H(^2S) reaction using the cc-pVDZ basis set[62] has been taken as a validation test for the GNOF-based BO AIMD. This reaction is important in the decomposition of ingredients in solid propellants used for rockets,[63] so it has been the subject of several high-quality theoretical studies[64–69] due to the experimental difficulty in preparing N atoms. These studies have shown that the reaction occurs via an abstraction mechanism dominated by the quartet GS, and presents a forward experimental barrier of 1.4 ± 0.3 eV.[70]

The initial separation between the nitrogen atom and the center-of-mass of H₂ was set at 6 Å,
and each trajectory was integrated until the separation between the final fragments was greater
than 6 Å. A time step of 0.1 fs was used, which yields to a conservation of the total energy with an
average error of 0.004 eV. Typical profiles of the kinetic, potential and total energies for a reactive
trajectory with a translational energy E<sub>T</sub> = 2.45 eV and H<sub>2</sub> at the GS can be seen in Fig. 1, while
in the supplementary material a movie of this trajectory can be found. From here we can conclude
that the collision occurs mainly in the time range from 20 fs to 40 fs.

In the N(^4S) + H₂(^1Σ) reaction, six are the strongly occupied NOs, namely three with occupancy
close to 2 and three singly occupied responsible for the quartet state. During reactive dynamics,
the lowest energy orbitals correspond to the 1s and 2s atomic orbitals of N. These NOs start
from the isolated atom and remain in the NH radical without significant changes. Similarly, two
of the 2p atomic orbitals of N undergo small transformations during the collision and continue
to maintain their character in the final NH, occupying directions perpendicular to the bonding
direction. Consequently, two NOs are responsible for the change of the bond pattern during the
collision, whose temporal evolutions are represented in Fig. 2 by specific snapshots. Thus, we
have a first NO that begins as a σ "ss" bonding orbital of the H₂ singlet with ON = 1.97 and
transforms into the σ "sp" bonding orbital of the NH triplet with ON = 1.96, while the other
individually occupied NO transforms from a 2p atomic orbital of the N into the 1s atomic orbital of the H.

Reaction probabilities and integral cross sections (ICSs) were calculated for translational energies up to 5.0 eV, and five different initial ro-vibrational states (ν, J) of H₂ molecule, namely the GS (0,0), the first two excited vibrational states (1,0) and (2,0), and two excited rotational states (0,10) and (0,15). Reduced batches of 1000 trajectories were run to determine a first estimate of the maximum impact parameter values (bmax). Finally, 5000 trajectories were carried out using appropriate value of bmax for each translational energy and ro-vibrational state (ν, J) of the H₂ diatom, which provided the ICS values through the following equation:

$$ \sigma \approx \pi b_{\max}^2 \frac{N_r}{N_t}, \qquad (15) $$

where N_r is the number of reactive trajectories, N_t the total number of trajectories for which the impact parameter satisfies $b \le b_{max}$, and $P_r = N_r/N_t$ is the total reaction probability. The number of propagated trajectories ensures a Monte Carlo statistical error of less than 5% for the ICS. Additionally, the final energy distribution of the product was obtained after a standard semi-classical determination of the vibrational, and rotational energies[71–73].

Fig. 3 presents bmax and σ as a function of E<sub>T</sub> for the five ro-vibrational states studied. For comparison, values reported in Ref.[68] based on the highly accurate DMBE PES[66] calculated at the MRCI-FVCAS/aug-cc-pVQZ level of theory are included. The bmax values were found to vary from 0.42 Å to 1.74 Å on increasing the collision energy. For the H₂ GS, the figure shows that the ICS has a threshold energy of approximately 1.5 eV and gradually decreases with increasing ν and J to 0.7 and 0.97 eV, respectively, in perfect agreement with the forward barrier observed for the reaction. Above the threshold energy, the ICSs increase almost linearly up to around 4 eV, where a much slower growth begins. For the ICSs corresponding to initial rotationally excited H₂ (J = 10,15), we observe that a decay then begins. Notice that the initial rotational excitation of H₂ molecule cause a more pronounced linear increase in ICSs at lower collision energies, and that the increase is significant in reactivity with increasing reagent vibrational excitation, as expected in this endoergic triatomic reaction with a late barrier. The shape of the ICS curves presented in Fig. 3 is similar to what was observed in previous studies, however our ICS values are lower. This ICS underestimation is understandable considering the modest basis set used in this work that lead to lower values of bmax.

Regarding the final ro-vibrational energy distribution of the product when the H₂ molecule is
FIG. 3. Maximum impact parameter ($b_{max}$) and integral cross section ($\sigma$) as a function of translational energy ($E_T$) for different H$_2$ ro-vibrational states ($v, J$). Dashed lines correspond to values reported in Ref.[68]

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

FIG. 4. Rotational-state distribution of product NH for the reaction $N + H_2(v=0, J=0) \rightarrow NH(v', J') + H$.

<figure><img src="image_5.png" /></figure>
initially in its ground state ($v=0, J=0$) our simulations show that below 3 eV, only the vibrational ground state ($v'=0$) and the first vibrational level ($v'=1$) are accessible to the NH molecule after its formation. At the highest collision energy studied in this work (5 eV), the highest vibrational level populated was $v'=5$. These results indicate a poor transfer between the translational and vibrational energy modes during the collision. However, an efficient rotational excitation of the NH molecule for the studied range of collision energies was observed. For a collision energy of 5 eV, 35 rotational levels are accessible, while for 1.5 eV, the six first rotational levels were populated. We show in Fig. 4 the calculated rotational-state distributions of the product for the reaction $N + H_2(v=0, J=0) \rightarrow NH(v',J') + H$ for six collision energies. The figure clearly indicates that the rotational excitation increases with increasing $E_T$. The higher the collision energy of the reactants, the higher the rotational energies of the product obtained. This behaviour is in perfect agreement with previous theoretical works [65, 67] on the studied reaction. The results presented in figures 3 and 4 emphasize the sensitivity of our method to qualitatively and quantitatively fully captures the physics involve in the dynamics of complex reactions.

In summary, we have shown that GNOF-AIMD is a method of choice to investigate the evolution of complex electronic problems, particularly reactive collisions. According to the existence theorems [1–3] of the 1RDM functional, there is a one-to-one mapping between the ground-state 1RDM and the ground-state N-particle density matrix, so by observing the real-time evolution of the NOs along with their ONs, we are seeing in real time the evolution of the solution to the electronic problem dynamically. The unique NO representation is especially useful for viewing the real-time evolution of changes in bond patterns. Indeed, NOs vary along trajectories calculated on the fly, adapting at each time step to the most favorable interactions of the corresponding nuclei configuration. GNOF-AIMD also allows to study dynamics with any transformation in the spins and the number of electrons of the component subsystems, conserving the total spin of the entire system. Therefore, GNOF-AIMD opens a promising field of research: AIMD based on natural orbital functionals.

# I. SUPPLEMENTARY MATERIAL

See the supplemental material for a movie of a reactive trajectory with translational energy $E_T$ = 2.45 eV and $H_2$ in the ground state.

**Acknowledgements:** Support comes from Grant PID 2021-126714NB-I00 funded by MCIN/AEI/
10.13039/501100011033 and the Eusko Jaurlaritza (Ref.: IT1584-22). The authors thank for technical and human support provided by IZO-SGI SGIker of UPV/EHU and DIPC.

[1] T. L. Gilbert, Phys. Rev. B **12**, 2111 (1975).

[2] M. Levy, Proc. Natl. Acad. Sci. USA **76**, 6062 (1979).

[3] S. M. Valone, J. Chem. Phys. **73**, 1344 (1980).

[4] M. Piris, in *Reduced-Density-Matrix Mech. with Appl. to many-electron atoms Mol.*, Vol. 134, edited by D. A. Mazziotti (John Wiley and Sons, Hoboken, New Jersey, USA, 2007) Chap. 14, pp. 387–427.

[5] R. A. Donnelly, J. Chem. Phys. **71**, 2874 (1979).

[6] E. V. Ludeña, F. J. Torres, and C. Costa, J. Mod. Phys. **04**, 391 (2013).

[7] M. Piris, in *Many-body approaches Differ. scales a Tribut. to N. H. March Occas. his 90th Birthd.*, edited by G. G. N. Angilella and C. Amovilli (Springer, New York, 2018) Chap. 22, pp. 261–278.

[8] S. M. Sutter and K. J. H. Giesbertz, Phys. Rev. A **107**, 022210 (2023).

[9] S. Di Sabatino, J. Koskelo, J. A. Berger, and P. Romaniello, Phys. Rev. B **107**, 035111 (2023).

[10] J. Cioslowski, C. Schilling, and R. Schilling, J. Chem. Phys. **158**, 084106 (2023).

[11] J. F. H. Lew-Yee and J. M. del Campo, J. Chem. Phys. **157**, 104113 (2022).

[12] S. Di Sabatino, J. Koskelo, J. A. Berger, and P. Romaniello, Phys. Rev. B **105**, 235123 (2022).

[13] B. Senjean, S. Yalouz, N. Nakatani, and E. Fromager, Phys. Rev. A **106**, 032203 (2022).

[14] S. Bousiadi, N. I. Gidopoulos, and N. N. Lathiotakis, Phys. Chem. Chem. Phys. **24**, 19279 (2022).

[15] J. Wang and E. J. Baerends, Phys. Rev. Lett. **128**, 013001 (2022).

[16] R. Schade, C. Bauer, K. Tamoev, L. Mazur, C. Plessl, and T. D. Kühne, Phys. Rev. Res. **4**, 033160 (2022).

[17] M. Rodríguez-Mayorga, K. J. H. Giesbertz, and L. Visscher, SciPost Chem. **1**, 004 (2022).

[18] J. Liebert, F. Castillo, J.-P. Labbé, and C. Schilling, J. Chem. Theory Comput. **18**, 124 (2022).

[19] Y.-F. Yao, W.-H. Fang, and N. Q. Su, J. Phys. Chem. Lett. **12**, 6788 (2021).

[20] D. Gibney, J. N. Boyn, and D. A. Mazziotti, J. Phys. Chem. Lett. **12**, 385 (2021).

[21] J. Liebert and C. Schilling, Phys. Rev. Res. **3**, 013282 (2021).

[22] J. Schmidt, M. Fadel, and C. L. Benavides-Riveros, Phys. Rev. Res. **3**, L032063 (2021).

[23] C. Schilling and S. Pittalis, Phys. Rev. Lett. **127**, 023001 (2021).

[24] J. M. Mercero, J. M. Ugalde, and M. Piris, Theor. Chem. Acc. **140**, 74 (2021).
[25] R. Quintero-Monsebaiz, L. I. Perea-Ramírez, M. Piris, and A. Vela, Phys. Chem. Chem. Phys. **23**, 2953 (2021).

[26] M. Rodríguez-Mayorga, I. Mitxelena, F. Bruneval, and M. Piris, J. Chem. Theory Comput. **17**, 7562 (2021).

[27] J. Cioslowski, J. Chem. Theory Comp. **16**, 1578 (2020).

[28] K. J. Giesbertz, Phys. Rev. A **102**, 052814 (2020).

[29] C. L. Benavides-Riveros, J. Wolff, M. A. Marques, and C. Schilling, Phys. Rev. Lett. **124**, 180603 (2020).

[30] I. Mitxelena and M. Piris, J. Phys.: Cond. Matter **32**, 17LT01 (2020).

[31] I. Mitxelena and M. Piris, J. Chem. Phys. **152**, 064108 (2020).

[32] Y. Lemke, J. Kussmann, and C. Ochsenfeld, J. Chem. Theory and Comp. **18**, 4229 (2022).

[33] J. F. H. Lew-Yee, J. M. del Campo, and M. Piris, J. Chem. Theory Comput. **19**, 211 (2023).

[34] M. Piris and I. Mitxelena, Comput. Phys. Commun. **259**, 107651 (2021), available at https://github.com/DoNOF/, documentation at https://donof.readthedocs.io/.

[35] J. F. H. Lew-Yee, M. Piris, and J. M. del Campo, J. Chem. Phys. **154**, 064102 (2021).

[36] K. Pernal and K. J. H. Giesbertz, Top Curr Chem **368**, 125 (2016).

[37] R. Schade, E. Kamil, and P. Blöchl, Eur. Phys. J. Spec. Top. **226**, 2677 (2017).

[38] I. Mitxelena, M. Piris, and J. M. Ugalde, in *State Art Mol. Electron. Struct. Comput. Correl. Methods, Basis Sets More*, Advances in Quantum Chemistry, Vol. 79, edited by P. Hoggan and U. Ancarani (Academic Press, 2019) Chap. 7, pp. 155–177.

[39] I. Mitxelena and M. Piris, J. Chem. Phys. **146**, 014102 (2017).

[40] I. Mitxelena and M. Piris, J. Chem. Phys. **153**, 044101 (2020).

[41] D. Marx and J. Hutter, *Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods* (Cambridge University Press, 2009) pp. 1–567.

[42] E. Tapavicza, G. D. Bellchambers, J. C. Vincent, and F. Furche, Phys. Chem. Chem. Phys. **15**, 18336 (2013).

[43] B. F. Curchod and T. J. Martínez, Chem. Rev. **118**, 3305 (2018).

[44] B. O. Roos, P. R. Taylor, and P. E. M. Sigbahn, Chem. Phys. **48**, 157 (1980).

[45] B. O. Roos, Adv. Chem. Phys. **69**, 399 (2007).

[46] B. O. Roos, P. Linse, P. E. M. Siegbahn, and M. R. A. Blomberg, Chem. Phys. **66**, 197 (1982).

[47] K. Andersson, P. A. Malmqvist, B. O. Roos, A. J. Sadlej, and K. Wolinski, J. Phys. Chem. **94**, 5483 (2010).
(1990).

[48] K. Andersson, P. A. Malmqvist, and B. O. Roos, J. Chem. Phys. **96**, 1218 (1992).

[49] P. Pulay, Int. J. Quantum Chem **111**, 3273 (2011).

[50] J. W. Park and T. Shiozaki, J. Chem. Theory Comput. **13**, 3676 (2017).

[51] I. Luzon, E. Livshits, K. Gope, R. Baer, and D. Strasser, J. Phys. Chem. Lett. **10**, 1361 (2019).

[52] P. Slavíček and T. J. Martínez, J. Chem. Phys. **132**, 234102 (2010).

[53] D. Hollas, L. Šištík, E. G. Hohenstein, T. J. Martínez, and P. Slavíček, J. Chem. Theory Comp. **14**, 339 (2018).

[54] M. Piris, Physical Review Letters **127**, 233001 (2021).

[55] I. Mitxelena and M. Piris, J. Chem. Phys. **156**, 214102 (2022).

[56] J. M. Mercero, R. Grande-Aztatzi, J. M. Ugalde, and M. Piris, Adv. Quantum Chem. **88**, 229 (2023).

[57] J. F. H. Lew-Yee, M. Piris, and J. M. Campo, J. Chem. Phys. **158**, 084110 (2023).

[58] M. Piris, Phys. Rev. A **100**, 32508 (2019).

[59] M. Piris, J. Math. Chem. **25**, 47 (1999).

[60] D. Beeman, J. Comput. Phys. **20**, 130 (1976).

[61] R. B. Bernstein, ed., *Atom-Molecule Collision Theory* (Springer US, 1979).

[62] T. H. Dunning and T. H. Dunning Jr., J. Chem. Phys. **90**, 1007 (1989).

[63] K. Prasad, R. A. Yetter, and M. D. Smooke, Combust. Sci. Technol. **124**, 35 (1997).

[64] S. Zhang and T. N. Truong, J. Chem. Phys. **113**, 6149 (2000).

[65] R. Z. Pascual, G. C. Schatz, G. Lendvay, and D. Troya, J. Phys. Chem. A **106**, 4125 (2002).

[66] L. A. Poveda and A. J. Varandas, Phys. Chem. Chem. Phys. **7**, 2867 (2005).

[67] X. Yu and Y. Yu, Integrated Ferroelectrics **147**, 67 (2013).

[68] Y. J. Yu and X. Y. Yu, Adv. Mater. Res. **924**, 227 (2014).

[69] J. Zhang, S. B. Gao, H. Wu, and Q. T. Meng, Chinese Phys. B **24**, 083104 (2015).

[70] M. Koshi, M. Yoshimura, K. Fukuda, H. Matsui, K. Saito, M. Watanabe, A. Imamura, and C. Chen, J. Chem. Phys. **93**, 8703 (1990).

[71] D. G. Truhlar and J. T. Muckerman, *Atom-molecule collision theory*, edited by R. B. Bernstein (Plenum, New York, Chapter 16, 1979).

[72] W. H. Miller, Adv. Chem. Phys **25**, 69 (1974).

[73] L. Bonnet and J.-C. Rayez, Chem. Phys. Lett. **277**, 183 (1997).