{'fraction_non_alphanumeric': 0.033129543952412424, 'fraction_numerical': 0.013384005287508261, 'mean_word_length': 4.696235294117647, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Photo-induced reactions from efficient molecular dynamics with electronic transitions using the FIREBALL local-orbital density functional theory formalism

The computational simulation of photo-induced processes in large molecular systems is a very challenging problem. Here, we present a detailed description of our implementation of a molecular dynamics with electronic transitions algorithm within the local-orbital density functional theory code Fireball, suitable for the computational study of these problems. Our methodology enables simulating photo-induced reaction mechanisms over hundreds of trajectories; therefore, large statistically significant ensembles can be calculated to accurately represent a reaction profile. As an example of the application of this approach, we report results on the [2+2] cycloaddition of ethylene with maleic anhydride and on the [2+2] photo-induced polymerization reaction of two C60 molecules. We identify different deactivation channels of the initial electron excitation, depending on the time of the electronic transition from LUMO to HOMO, and the character of the HOMO after the transition.


I. INTRODUCTION
One of the greatest challenges in simulation methodology is to accurately and efficiently interpret photo-induced reaction mechanisms which are critically important in many chemical processes, such as vision, photodamage, UV absorption, photosynthesis, photography, etc. In Born-Oppenheimer (or adiabatic) molecular dynamics -the nuclei follow classical trajectories defined on a single potential energy surface (PES); this PES corresponds to the ground state electronic energy of the system, for frozen atomic configurations. To accurately simulate photo-induced reactions, the excited potential energy surface must be appropriately accessed beyond the Born-Oppenheimer approximation. Additionally, the ability of the electron to transition between the ground PES and the excited PES at the point of the conical intersection between the two PESs is very important in accurately portraying the transition that will occur in photo-induced reactions.
Broadly speaking, two types of non-adiabatic methods dominate the literature: Ehrenfest Dynamics (ED) methods and surface hopping (SH) methods. [1][2][3] In the ED approach, the nuclei move classically on a single effective PES obtained by averaging over all the adiabatic states involved; however, in SH methods the classical degrees of freedom (nuclei coordinates) evolve on single adiabatic surfaces, 1 and make probabilistic hops from one PES to another. Surface-Hopping methods present key advantages over ED methods, because the system is always on a certain pure adiabatic state, not in an average non-physical state. Moreover, ED calculations sometimes lead to nonphysical features, such as energetically inaccessible levels populated or violation of microscopic reversibility. 4 Photo-induced processes in large molecular systems are computationally intensive to simulate due to the complexities in modelling the vast search space associated with photo-activated reaction mechanisms. Also, the search for conical intersections is quite more complex than geometry optimization. [5][6][7][8] Photo-induced polymerization reactions, such as the photo-induced [2+2] cycloaddition reactions of two C 60 molecules, would easily require large simulation systems containing hundreds of atoms. Biological photo-induced reactions are even more demanding in computational resources. 9 Additionally, photo-induced reaction mechanisms are probability events.
Hence, to accurately simulate reaction mechanisms, hundreds, if not thousands, of trajectories are required to obtain the proper statistical information. Thus, efficient algorithms are needed, at least, at an exploratory level. There is a definitive need for simulation tools that can, at minimum, obtain sufficiently accurate results that will efficiently yield significant insight in photo-induced reaction mechanisms without overburdening computational resources.
Tully's molecular dynamics with electronic transitions (MDET) algorithm 10 employes a surface hopping strategy 11 to switch between different adiabatic states. In Tully's MDET approach, the system evolves on one of the PES corresponding to a particular adiabatic state E A of the system.
Using the idea of ensemble of trajectories, the transitions between different PES are determined through a probabilistic approach; in particular, the transitions between adiabatic electronic states are determined using the fewest-switches algorithm that considers the minimum number of nonadiabatic transitions that are compatible with the correct statistical distribution of state populations at all times. The MDET method is probably the most used and successful method to deal with problems in which coupling between electrons and ions is required.
As such, we have focused primarily on implementing a MDET algorithm within FIREBALL, an efficient local-orbital DFT MD technique. 25 Our approach to the nonadiabatic MD simulaton is extremely efficient because we pre-calculate and store the integrals needed during the MD simulation (e.g. the integrals for the calculation of the hamiltonian matrix elements). 26,27 During the simulations, interactions are loaded into memory from the pre-computed data files and no additional time is needed for on-the-fly integration. We have previously demonstrated a successful application of the FIREBALL DFT-MDET in the photo-isomerization processes of stilbene and azobenzene. 19 The remainder of this paper is organized as follows. In Section II, we present a detailed description of our implementation of DFT-MDET within the FIREBALL code. In Section III, we present results for two different photo-induced cycloaddition reactions: the [2+2] cycloaddition of ethylene with maleic anhydride to form cyclobutane-1-2-dicarboxilic anhydride and the [2+2] photo-induced polymerization reaction of two C 60 molecules. Finally, in section IV, we summarize the main conclusions.

II. DESCRIPTION OF THE METHOD
In DFT-MDET simulations the atoms follow classical trajectories {R α (t)} on one of the PESs corresponding to a particular adiabatic state E A . The coupling between the nuclear motion and the electronic quantum state is taken into account by means of probabilistic hops between different adiabatic PESs (E A → E B ). The probability for these hops is determined from the time evolution of the electronic states and their non-adiabatic coupling vectors (see below). In the following discussion, we present details about our implementation of the MDET algorithm within the FIREBALL local-orbital DFT code. Step 1: Solution of the stationary Kohn-Sham Hamiltonian. In this first step we obtain the self-consistent (SCF) solution of the stationary Kohn-Sham (KS) Hamiltonian for the current atomic configuration and adiabatic state: where R represents the atomic positions {R α (t)} at time t and ψ j (r; R) are the single-particle adiabatic KS orbitals (at time t). The different adiabatic states, E A , are defined by means of the corresponding set of occupied KS adiabatic orbitals or, in other words, by the occupancies f k of the ψ k orbitals, with f k = 1 (occupied) or 0 (empty). Thus, the KS Hamiltonian in eq.
(1) is defined using the electronic configuration for the current PES (e.g. the electron density is ρ = k f k |ψ k | 2 ); for this purpose, we perform a constrained DFT calculation in which the level occupancy is imposed for a particular PES to maintain the electronic configuration during the SCF cycle. In order to keep track of the identity of the adiabatic KS states at different time steps we calculate the overlaps ψ k (t − ∆t)|ψ m (t) .
Step 2: Calculation of the non-adiabatic couplings. The time evolution of the electronic state is calculated using time-dependent KS theory. [28][29][30] For this purpose, we expand the time-dependent KS orbitals, ϕ i , in terms of the adiabatic KS states and evolve them using the time-dependent, Schrödinger-like, KS equations: This yields the equations of motion for the coefficients c ik (t) in eq. (2): In this equation the coupling between the classical motion of the atoms and the time evolution of the electronic quantum state is reflected in the non-adiabatic coupling (NAC) term: where V α = ∂R α /∂t is the atomic velocity of atom α and d α km are the non-adiabatic coupling vectors (NACV) between adiabatic KS states k and m: The NACVs d α km are the central quantity for non-adiabatic molecular dynamics. The atomic motion induces changes in the adiabatic states populations |c ik (t)| 2 through the coupling of the atomic velocities and the NACVs, see eq. (4). Recently, we have derived an expression to calculate the NACVs d α km in a basis set of local-orbitals on the fly along the MD simulation in a practical and computationally efficient way. 31 Alternatively, the NAC for the time evolution (eq. (4)) of the coefficients c ik (t) can be directly calculated using the following numerical derivative 32 : which gives the value of the NAC at time (t − ∆t/2).
Step 3: Time evolution of the electronic states. Once the NACs have been calculated, we propagate the coefficients c ik from time (t − ∆t) to time t using a small time step ∆t el and a 4-th order Runge-Kutta numerical scheme to integrate eq. (4). Typically ∆t el ∼ ∆t/500. This is the inner loop for the propagation of the electronic states for each step, ∆t, in the atomic motion.
Step 4: Surface hopping. At each time step ∆t, a transition between adiabatic PESs E A → E B may occur through a probabilistic "hop" using the fewest-switches algorithm. 10 The MDET algorithm incorporates the following switching probability g jk (t) to define the transition between where the coefficient b kj is defined as: The transitions are determined using a Monte-Carlo approach following the probabilities g jk (t).
When a transition takes place, the occupancies of the adiabatic KS orbitals ψ j and ψ k are modified: f j : 1 → 0; f k : 0 → 1, and the stationary KS Hamiltonian, eq. (1), is re-calculated for the new electronic configuration (see step 1). In this case, energy conservation is imposed by re-scaling velocities along the direction of the NACV, eq. (6)). 33,34 Hops to higher energy PESs are rejected if insufficient kinetic energy is present. Velocity rescaling and hop rejection ensures detailed balance between transitions up and down in energy. If the switch does not occur, the simulation continues on the current PES.
Step 5: Motion of the atoms. Finally, the atomic positions are updated {R α (t)} → {R α (t + ∆t)} using the forces corresponding to the current PES, E A : where F α is the force acting on atom α. The calculation goes back to step 1, using the new atomic positions. This scheme is repeated until the simulation is finished.
We now describe the parameters used during simulations of photo-cycloaddition (PCA) reactions discussed in the next section. The optimized atomic structure in the ground state is used as the initial molecular configuration for the reactant species. Atomic positions of the reactant species were relaxed with convergence criteria for the total energy and forces of 10 −6 eV and 0.05 eV/Å respectively to obtain the optimized structures. All calculations were performed using local density approximation (LDA). The optimized numerical atomic-like orbitals (NAOs) were confined to regions limited by the corresponding cutoff radii r c . 26,35 The following r c for the different To gain more insight into the reaction mechanisms, we performed, in total, 50 independent DFT-MDET simulations, each with total time of 2 ps; in each simulation, a single electron was initially promoted from the HOMO to the LUMO state of the molecular complex ( Each simulation took approximately 16 hours on single CPU Intel Xeon E5420 @ 2.50GHz which demonstrates the computational efficiency of our FIREBALL-MDET code. All simulations were initiated from an identical configuration -both planar molecules were oriented parallel and the carbon-carbon double bond (C 3 = C 4 ) of ethylene was positioned above the carbon-carbon double bond (C 1 = C 2 ) of maleic anhydride (see Fig. 2B top). The initial distance between the molecular planes of the anhydride and ethylene was 3Å. The velocities were initialized with a normal distribution corresponding to an initial temperature of 500 K. In order to show clearly the stochastic nature of the algorithm, the same initial positions and velocities were applied for each trajectory. Four frontier orbitals (HOMO-1, HOMO, LUMO, LUMO+1) were propagated in time and thus the electron transitions were allowed only amongst these states.
We were able to identify three different deactivation channels of the initial electron excitation A direct addition reaction of reactants was observed in 12 cases. The detailed analysis of the reaction pathway leading to direct formation of cyclobutane diacid anhydride is presented in Fig. 3.  Fig. 3 b). The time evolution of the energy spectrum for the four frontiers molecular orbitals with their temporal occupancies is shown in Fig. 3 b). Complementary, Fig. 3 c) shows the PES corresponding to the ground S 0 and first excited S 1 state along the simulation; the cyan dots indicate the actual PES for the simulation. The grey line represents the evolution of the NACVs between the HOMO and the LUMO orbitals with its maximum occurring near the conical intersection on reffiganh-1a c).
To characterize the pathway of the reaction process, we traced a time evolution of a characteristic distance d reac and force F reac acting between reactants. The characteristic distance d reac measures the distance from the center of ethylene carbon-carbon bond (C 1 = C 2 , see Fig. 2) to the center of anhydride carbon-carbon bond (C 3 = C 4 ). We also evaluated the forces acting between the two molecules as F reac = F ethylene − F anhydride projected along the vector from the center of ethylene carbon-carbon bond to the center of anhydride carbon-carbon bond ( d reac ).
The distance d reac and force F reac between reactants are plotted in Fig. 3  We should note that, in this particular case, the switching mechanism does not occur at the distance with the NACV maximum, but later due to the stochastic nature of the fewest switches  Fig. 6 and Fig. 7. In this case the system typically finalizes (22 cases from 50 total) in the decoupling of the molecular reactants but we detected 3 cases where the intermediate case produces the complete addition reaction (see Fig. 6). One typical example of the transition state, which evolves into the full dissociation is presented in Fig. 7. At the beginning of this transition state reaction profile, the time evolution of the system is the same as in the previous dissociative cases (e.g. see As a result, the ethylene molecule dissociates from the anhydride as seen from the increase of the distance d reac plotted on Fig. 7d). For the rest of the simulation, the system is in the ground state and the HOMO orbital is localized on the ethylene molecule.

B. Polymerization of two C 60 molecules via cycloaddition
Here, as a second example of the application of our FIREBALL-MDET code, we present simulations of the polymerization process of two free-standing C 60 molecules via a 2+2 PCA reaction.
This reaction was first reported by Rao et al. in 1993. 37 The basic principles for this case are similar to our discussion in the previous section, see Fig. 2 Accordingly, the HOMO of the molecular complex C 60 -C 60 has an anti-bonding character and the bonding orbital is associated to the LUMO+1 (see Fig. 8C). We should note that the LUMO and LUMO+1 orbitals are nearly degenerate with an energy difference of 0.08 eV in the initial configuration. Nevertheless the LUMO+1 shows a strong bonding character with pronounced electron densities localized between the frontiers atoms of the two C 60 molecules.
Here we consider the case where two C 60 molecules are oriented by 56/65 edges ( see Fig. 8 B) against each other. 38 In the first step, we examine the PESs as function of distance between the two C 60 molecules in the ground S 0 , first excited S 1 , and second excited S 2 states for distances in the range from 1.3 to 3.0Å, with incremental changes in the distance of 0.05Å. In this calculations we fixed z and y coordinates of the four interacting carbon atoms forming the 56/65 edge see Fig. 8 B), i.e. those involved in the cycloaddition chemical bond, and the rest of the atomic positions were allowed to relax. The calculated PES S 0 , S 1 and S 2 are plotted in Fig. 8 A). We can identify two minima on the S 0 PES: (i) a global one located at the distance of ≈ 2.69Å, which corresponds to two weakly interacting C 60 molecules; (ii) and the second, a local minimum found at the distance of ≈ 1.56Å. This second minimum corresponds to two chemically bonded C 60 molecules. In between these two energy minima, there is an energy barrier of ≈ 0.6 eV located at a distance of 1.75Å, which avoids a spontaneous transition from the energetically higher chemically bound state to the weakly bound state. At this distance there is a conical intersection between S 0 and S 1 .
In our DFT-MDET simulations, we initially placed two C 60 molecules at a distance of 2.0Å (vertical line in Fig. 8 A) i.e. 0.25Å away from the conical intersection, and one electron was promoted from the HOMO to the LUMO+1 state. During the simulation, we let the system evolve in time for 0.5 ps at a constant temperature of 100K, with initial velocities chosen as discussed in the previous section. Electron transitions are allowed between three orbitals -HOMO, LUMO and LUMO+1. Our simulations identify two scenarios: (i) a PCA reaction forming a C 60 bonded dimer; or (ii) a complete decoupling of the two C 60 molecules towards the weakly bound minimum. In total, we carried 20 trajectories and 3 of the simulations yielded the polymerization process; 17 simulations lead to the separation of the C 60 molecules. Fig. 9 shows the time evolution for a case where electron de-excitation happens via the 2+2 cycloaddition reaction. Fig. 9 a) represents the HOMO, LUMO and LUMO+1 KS orbitals at three different snapshots corresponding to the vertical black lines in Fig. 9 b). Fig. 9 b) shows the time evolution of the HOMO, LUMO and LUMO+1 KS eigenvalues. Fig. 9 c) shows the PESs corresponding to the ground S 0 , first S 1 and second S 2 excited states throughout the simulation.
The cyan circles identify on which PES the system is located at a given time during the simulation.
The solid/dashed grey line represents the evolution of the NACVs between the LUMO and HOMO and between the LUMO+1 and LUMO, respectively. Finally, Fig. 9 d) displays the evolution of the force F reac acting between the two molecules and their distance d reac , which is defined as an average distance between the frontiers atoms on the 56/65 edges of each molecule. Similarly, the force F reac is the difference between the sum of forces acting on each C 60 molecule projected along the vector d reac .
Initially, the S 1 and S 2 states are almost degenerate, see Fig. 9 c). Very rapidly in the simulation, the LUMO and LUMO+1 energy states show a strong mixing and a large NACV, and an electronic hopping between LUMO+1 and LUMO occurs (at t = 4.75 fs). The system evolves now in the An example of the non-reactive pathway is presented in Fig. 10. Initially, the system evolves very similarly to the previous case where the system appraoches a conical intersection between S 1 and S 2 ; due to the probabilistic nature of the fewest switches method, the electronic transition from LUMO+1 to LUMO takes place now a little later, at t = 5.75 fs, see Fig. 10 b). The system evolves then in the S 1 state, moving towards a S 0 -S 1 conical intersection, but the electron hopping does not occur now. Thus, the systems stays in the S 1 PES, with a repulsive force F reac acting between both molecules; the distance d reac oscillates (Fig. 10 d), due to the attraction associated with the bonding LUMO (occupied with one electron), see Fig. 10 a). During the second oscillation cycle, when d reac decreases again the S 1 and S 0 PESs (and the LUMO and HOMO energy levels) approach each other again and the electron hop takes place now, at time t = 100 fs (Fig. 10 b,c).
The HOMO (occupied by two electrons) presents an anti-bonding character now (see Fig. 10 a), t = 110 fs), enhancing the repulsive force F reac . Consequently the distance d reac increases and the system goes to the weakly interacting ground state (Fig. 8) and the PCA reaction does not occur.

IV. CONCLUSIONS
The computational simulation of photo-induced processes in large molecular systems is a very challenging problem. Firstly, to properly simulate photo-induced reactions the PESs corresponding to excited states must be appropriately accessed; secondly, understanding the mechanisms of these processes requires the exploration of complex configurational spaces and the localization of conical intersections; finally, photo-induced reactions are probability events, that require the simulation of hundreds of trajectories to obtain the statistical information for the analysis of the reaction profiles.
Here, we have presented a detailed description of our implementation of a DFT-MDET technique within our local-orbital code FIREBALL, suitable for the computational study of these problems. The MDET technique consists of two nested time loops: an outer loop for the motion the atoms and an inner loop for the propagation of the electronic states. The atoms move in adiabatic PESs, that are calculated using constrained DFT calculations for the different electronic configurations. The time evolution of the electronic states is calculated using time-dependent KS theory.
The coupling of the atomic motion and the time evolution of the electronic states is taken into account by means of probabilistic hops between different PESs, using the fewest switches algorithm.
When hops take place, energy conservation is imposed re-scaling the velocities along the direction of the non-adiabatic coupling vectors, that are calculated on the fly along the molecular dynamics simulation 31 .
As an example of the application of this DFT-MDET approach, we have presented some simulations for two different photo-induced cycloaddition reactions: the reaction of two small organic molecules (maleic anhydride and ethylene, forming cyclobutane-1-2-dicarboxilic anhydride) and the polymerization reaction of two C 60 molecules. In the first case, an electron is initially promoted from the HOMO to the LUMO; the LUMO presents a bonding character, inducing an attractive force between both molecules. At a time of ∼ 60 fs, the system goes through a S 1 -S 0 conical intersection, with a large value of the NACV between HOMO and LUMO. We identify three different deactivation channels of the initial electron excitation, depending on the time of the electronic transition from LUMO to HOMO, and the character of the HOMO after the transition. In particular, the cycloaddition reaction is established when the HOMO is a bonding state after the transition.
In the case of the larger C 60 molecules, the S 1 and S 2 states are initially almost degenerate and we promote one electron from the HOMO to the LUMO+1. Very rapidly in the simulations (∼ 5 fs) the electron hops from the LUMO+1 to the LUMO and the system evolves in the S 1 state towards a S 1 -S 0 conical intersection. In this part of the simulations there is an overall repulsive force between both molecules and an attractive force between the frontier atoms on each molecule due to the excited electron in the bonding LUMO state; thus, the molecules tend to elongate. We identify two different scenarios, depending on the time of the electronic transition from LUMO to HOMO and the character of the HOMO state after the transition: (i) a PCA reaction forming a bonded C 60 dimer; and (ii) the motion of the system towards the weakly coupled ground state. to the center of maleic anhydride carbon-carbon bond (C 3 − C 4 ) (red) and forces acting between reactants (blue).