Time evolution of natural orbitals in ab initio molecular dynamics

This work combines for the ﬁrst 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( 4 S) + H 2 ( 1 Σ ) → NH( 3 Σ ) + H( 2 S) is taken as validation test. Collision energy inﬂuences 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 deﬁne a natural orbital functional (NOF). 4 An important issue is that the approximate NOFs

One-particle reduced density matrix (1RDM) functional theory [1][2][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). 4An 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 Exhaustive reviews of NOF-based methods can be found elsewhere.[36][37][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 nuclear 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 generally accurate.There are other problems, for example, when the dynamics starts in an excited electronic state, where strong couplings between two or a) Electronic mail: alejandro.rivero@univ-lille.frb) Electronic mail: mario.piris@ehu.eusmore 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.
7][48][49] Unfortunately, multireference 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 discontinuities 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 functionals.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 functions than wavefunction-based methods for systems with a large amount of strong non-dynamic correlation.Interestingly, NOFs further corroborate the motivation behind the floating occupation molecular orbital complete active space configuration interaction (FOMO-CASCI) method, a promising alternative to CASSCF in dynamics simulations. 52,53he 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 multiconfigurational character.In addition, GNOF has proven to be successful in dealing with aromaticity 56 and charge delocalization error. 57GNOF 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 Nelectron molecule with N I spin-unpaired electrons, N II = N − N I spin-paired, and total spin S. We focus on the state of highest multiplicity: 2S + 1 = N I + 1. 58 For the whole ensemble, the expected value of Ŝz is zero; so the spin-restricted theory is adopted, i.e., a single set of orbitals is used for α and β spins, which results in doubly occupied spatial orbitals.
The orbital space is divided in turn into two subspaces: subspaces Ω g , each of which contains one orbital |g with g ≤ N II /2, and N g orbitals |p with p > N II /2.The total occupancy for a given subspace Ω g ∈ Ω II is 2, and N g corresponds to the maximum possible value determined by the number of basis functions (N B ).On the other hand, Ω I is composed of N I mutually disjoint subspaces, each of which contains only one orbital g with 2n g = 1.This orbital is individually occupied, but we do not know whether the electron has α or β spin.
A reconstruction functional for the 2RDM in terms of the ONs leads to the following electronic energy (GNOF): The intra component is formed by the sum of the pair and single-electron energies, namely where N Ω =N II /2 + N I , H pp are the diagonal one-electron matrix elements of the kinetic energy and external potential operators, and Π (n q , n p ) = √ n q n p δ qΩ a δ pΩ a − δ qg − δ pg .Ω a denotes the subspace composed of orbitals above the level N Ω , whereas J pq = pq|pq and L pq = pp|qq are the Coulomb and exchange-time-inversion integrals, 59 respectively.The inter-subspace HF term is where K pq = pq|qp are the exchange integrals.The prime in the summation indicates that only the inter-subspace terms are taking into account.The inter-subspace static component is written as where Φ p = n p h p with the hole h p = 1 − n p .Finally, the inter-subspace dynamic energy is ) In Eq. ( 6), Ω b II denotes the subspace composed of orbitals below the level N II /2, and the dynamic part of the ON n p is defined as with h c = 0.02 √ 2. n d p 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 Z A represents the atomic number of nucleus A, and R AB is the distance between nuclei A and B. Considering that all NOs φ i are expanded in a fixed atomic basis set, φ i = ∑ υ C υi ζ υ , the derivative of the total energy with respect to the coordinate x of nucleus A is given by 40 where Γ µυ and D µηυδ are the 1RDM and 2RDM, respectively, S µυ = µ|υ is the overlap matrix, and λ µυ are the Lagrange multipliers obtained from RDMs, all in the atomic orbital representation.The first term of Eq. ( 8) 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 υi on 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. 39ll derivatives in Eq. ( 8) have an explicit dependence on the nuclear coordinate x A , so the force acting on each nucleus A (F A = −∇ 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 implemented in DoNOF, 34 which allows the calculation of nuclear trajectories by determining "on the fly" the forces using NOF gradients (8).Beeman's algorithm 60 is used to numerically integrate Newton's equations of motion, whereas the initial conditions are obtained using a standard Monte Carlo sampling procedure. 61he N( 4 S) + H 2 ( 1 Σ) → NH( 3 Σ) + H( 2 S) 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][65][66][67][68][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. 70he initial separation between the nitrogen atom and the center-of-mass of H 2 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 T = 2.45 eV and H 2 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( 4 S) + H 2 ( 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 2 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 2 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 val-ues (b max ).Finally, 5000 trajectories were carried out using appropriate value of b max for each translational energy and rovibrational state (ν, J) of the H 2 diatom, which provided the ICS values through the following equation: where N r is the number of reactive trajectories, N t the total number of trajectories for which the impact parameter satisfies b ≤ 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][72][73] .Fig. 3 presents b max and σ as a function of E T 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 b max values were found to vary from 0.42 Å to 1.74 Å on increasing the collision energy.For the H 2 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 2 (J = 10,15), we observe that a decay then begins.Notice that the initial rotational excitation of H 2 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   Regarding the final ro-vibrational energy distribution of the product when the H 2 molecule is initially in its ground state (ν = 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 (ν = 0, J = 0) → 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 prefect 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][2][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.

FIG. 2 .
FIG. 2. Time evolution of the two strongly occupied natural orbitals involved in the bond pattern change during the collision.

FIG. 3 .
FIG.3.Maximum impact parameter (b max ) and integral cross section (σ ) as a function of translational energy (E T ) for different H 2 rovibrational states (ν, J).Dashed lines correspond to values reported in Ref.68 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 b max .