Non-Adiabatic Ab Initio Molecular Dynamics with Floating Occupation

Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

pubs.acs.org/JCTC

Nonadiabatic Ab Initio Molecular Dynamics with the Floating Occupation Molecular Orbital-Complete Active Space Configuration Interaction Method Daniel Hollas,† Lukás ̌ Šištík,† Edward G. Hohenstein,‡,§ Todd J. Martínez,∥,⊥ and Petr Slavíček*,†,# †

Department of Physical Chemistry, University of Chemistry and Technology, Prague, Technická 5, 16628 Prague 6, Czech Republic Department of Chemistry and Biochemistry, The City College of New York, New York, New York 10031, United States § PhD Program in Chemistry, The Graduate Center of the City University of New York, New York, New York 10016, United States ∥ Department of Chemistry and the PULSE Institute, Stanford University, Stanford, California 94305, United States ⊥ SLAC National Accelerator Laboratory, Menlo Park, California 94025, United States # J. Heyrovský Institute of Physical Chemistry, Dolejškova 3, 18223 Prague 8, Czech Republic ‡

S Supporting Information *

ABSTRACT: We show that the floating occupation molecular orbital complete active space configuration interaction (FOMO-CASCI) method is a promising alternative to the widely used complete active space self-consistent field (CASSCF) method in direct nonadiabatic dynamics simulations. We have simulated photodynamics of three archetypal molecules in photodynamics: ethylene, methaniminium cation, and malonaldehyde. We compared the time evolution of electronic populations and reaction mechanisms as revealed by the FOMO-CASCI and CASSCF approaches. Generally, the two approaches provide similar results. Some dynamical differences are observed, but these can be traced back to energetically minor differences in the potential energy surfaces. We suggest that the FOMO-CASCI method represents, due to its efficiency and stability, a promising approach for direct ab initio dynamics in the excited state.

■

(TDDFT), or coupled cluster (CC) based approaches.8 Although these methods typically provide a reasonable description of the Franck−Condon region, they struggle to reliably model bond breaking, biradical intermediates, or crossing between potential energy surfaces.9 Only relatively recently have single-reference methods been adapted for applications in photodynamics,10 and new approaches to circumvent some of the problems with single reference methods have been proposed (e.g., using references with different spin9,11−14 or particle number15). However, the most reliable approaches remain the multireference methods such as the complete active space selfconsistent field (CASSCF) and their correlated variants such as multireference configuration interaction (MRCI), second order multireference perturbation theory (CASPT2), or semiempirical methods including dynamical correlation built upon the CASSCF wave function.16 In practice, multireference methods can suffer from high computational cost, sensitivity to active space selection, and instabilities in the calculations due to

INTRODUCTION Since their appearance in the 1960s, molecular dynamical simulations have gradually gained an uncontested position in modeling dynamical and thermodynamical properties of molecules and materials. Because of the rapid development of computer technology, ab initio molecular dynamics has become a standard tool for simulations of reactive chemical processes. Such simulations are especially important in photodynamics where chemical intuition often fails. A further challenge in photodynamical simulations is the ubiquitous breakdown of the Born−Oppenheimer approximation. In last 25 years, various methods for simulations of nonadiabatic processes have been developed. These methods include both semiclassical approaches, such as Ehrenfest dynamics1 or surface hopping (SH)2 schemes, as well as more elaborate wavepacket-based schemes such as full multiple spawning (FMS)3,4 or variational multiconfiguration Gaussian wavepackets (vMCG).5,6 The major limitation of photodynamics simulations remains an accurate description of the underlying potential energy surfaces.7 Various electronic structure methods have been used to model excited-state potential energy surfaces. These include single-reference methods, such as restricted open-shell Kohn− Sham (ROKS), time-dependent density functional theory © XXXX American Chemical Society

Received: September 12, 2017 Published: December 5, 2017 A

DOI: 10.1021/acs.jctc.7b00958 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

where ni is the occupation number of molecular orbital i and Cιμ are orbital coefficients. In the FON-HF method, the occupation number of a given orbital changes continuously according to the energy of the orbital. For example, the orbital occupation numbers can be given by a Fermi−Dirac distribution function36

changes in the active space orbitals that can lead to energy discontinuities along the trajectory. In the CASSCF method, the wave function is expanded as a linear combination of configuration state functions (CSFs):17 CSFs

ψ CASSCF =

∑ ckϕk(C)

ni =

(1)

k

Two sets of coefficients are then optimized: orbital coefficients C and configuration expansion coefficients ck. In photodynamics, often a single set of orbital coefficients is used for all electronic states (state-averaged CASSCF). The CASSCF method represents a full configuration interaction (FCI) method within a preselected active space, i.e., CSFs corresponding to all excitations within an active space are considered. The active space is typically rather small, and the method is thus best at describing so-called “static” correlation.18 The CASSCF method can be viewed as a multiconfigurational analogue of the Hartree−Fock method. Unfortunately, the computational cost of the CASSCF method grows exponentially with the active space. Furthermore, the highly nonlinear orbital optimization increases the computational cost and can lead to multiple local minima. Avoiding the orbital optimization step significantly simplifies the problem. This is achieved in the complete active space configuration interaction (CASCI) method, where only the configuration expansion coefficients are optimized. Indeed, there have been many attempts to develop an accurate CASCI method based on preprocessed orbitals. Hartree−Fock orbitals are usually a poor choice for the subsequent CASCI calculation, whereas Kohn−Sham orbitals are somewhat better. Several other approaches have been proposed, such as natural orbitals generated by the unrestricted HF method (UNO orbitals),19,20 high multiplicity natural orbitals (HMNO),21 improved virtual orbitals (IVO),22 configuration-averaged HF orbitals,23 or natural orbitals generated from configuration interaction singles (CISNO).24 Here, we advocate the use of floating occupation molecular orbitals (FOMO) from a finite-temperature self-consistent field calculation (this is often called fractional occupation number or FON-SCF).25−27 In connection with the CASCI approach, these orbitals were first used in a semiempirical context,25,28−31 but the method was later extended to ab initio wave functions.26,27 More recently, analytical ab initio FOMOCASCI gradients27 and nonadiabatic coupling vectors32 were formulated and implemented in the GPU-based code TeraChem.33,34 The GPU implementation of the FOMOCASCI method has enabled ab initio molecular dynamics simulations of photoexcited systems in the condensed phase including hundreds of atoms.35 In this paper, we further test the use of FOMO-CASCI for nonadiabatic surface hopping (SH) simulations and benchmark the method against CASSCF results for several archetypal systems. Let us begin with a summary of the method. The total ensemble HF (FON-HF) energy can be expressed as E HF =

∑ Pμνhμν+ ∑ μν

μνλσ

fi (ε) =

i

(4)

1 β

2 −[(ε − εi)2 /2β 2] e π

(5)

where the occupation number of orbital i is obtained by the integration of the distribution function up to the Fermi level ni =

εf

∫−∞ fi (ε) dε

(6)

These two approaches provide essentially identical results.26 The broadening coefficient β, which is analogous to an electronic temperature, is usually set to a constant value. In the context of CASCI, it can also be incorporated into the optimization cycle as a variational parameter.37 The FOMO orbital optimization closely resembles the usual HF procedure with the set of occupation numbers updated after each SCF iteration. At convergence, the self-consistent orbitals are used as a basis for optimization of the configuration interaction wave function in FOMO-CASCI. In this work, we test the applicability of the FOMO-CASCI approach as an electronic structure method for computational photodynamics. First, we briefly summarize the technical details of our simulations. Then, we focus on the photodynamics of several archetypal reactions, including the cis−trans isomerization of ethylene, the three-state problem of the methyliminium cation, and finally excited-state proton transfer in the malonaldehyde molecule.

■

COMPUTATIONAL DETAILS We have modeled the nonadiabatic dynamical processes using a surface hopping (SH) scheme. The approach is based on a semiclassical ansatz, where atomic nuclei evolve classically on a single adiabatic potential energy surface; however, the particle is characterized by coefficients ck(t) defining the electronic population of the particles in the kth electronic state. The coefficients ck(t) evolve in the adiabatic representation as d c k (t ) = −iℏEk (t ) − dt

∑ cl(t )hckl(R(t ))·v c(t ) l

(7)

hckl(R(t))

where is the nonadiabatic coupling vector between states k and l in geometry R(t) and vc(t) is the nuclear velocity vector. The fewest switches algorithm2 was employed to calculate the probability for jumping between electronic states. The decoherence correction of Granucci and Persico was applied at every time step with the decoherence parameter α set to its recommended value of 0.1 atomic units.38,39 The goal of this paper is to evaluate the FOMO-CASCI method for the purpose of nonadiabatic simulations. For comparison, we have also simulated the photodynamical

⎛ ⎞ 1 PμνPλσ ⎜(μν|λσ ) − (μσ |λν)⎟ ⎝ ⎠ 2

where hμν and (μν|λσ) represent one- and two-electron integrals, respectively, and Pμν is the density matrix

∑ niCiμCiν

1+e

where εi and εf represent energy of orbital i and Fermi energy, respectively, kB is the Boltzmann constant, and T is the thermodynamic temperature. The factor of 2 in the numerator is related to the starting molecular orbital occupation; here, we consider closed-shell singlet wave functions. In this work, we used an alternative Gaussian broadening scheme

(2)

Pμν =

2 (εi − εf )/ kBT

(3) B


Article


Figure 1. Potential energy surface scans along two main reaction coordinates for the ground and first excited states of ethylene. (A) Geometric distortions characterizing the reaction coordinates (all other coordinates were kept fixed). Potential energy scans of the ground and first excited states as well as the contour plot of the excitation energies at each geometry of the scan (rightmost panel) are shown in B and C.

Unless otherwise noted, the energy conservation issue occurred in less than 10% of trajectories. The initial conditions (positions and momenta) for surface hopping simulations were generated by sampling the Wigner transform of the initial harmonic vibrational wave function. The normal modes and frequencies were calculated at the MP2/631++G** level in the MOLPRO package. For each molecule and method, one hundred surface hopping trajectories were carried out with a time step of 5 au (∼0.121 fs). The total length of the simulation was set to 200 fs for ethylene, 100 fs for methaniminium, and 300 fs for malonaldehyde. The integration of the classical equations of motion was performed with the velocity Verlet algorithm. We estimated the standard error (SE) of the electronic state populations via

processes using the corresponding CASSCF electronic wave function. For both methods to be compared directly, the dynamical runs were started from the same set of structures and from the same set of random number seeds. The active space is denoted as (X,Y), where X is the number of active electrons and Y is the number of active orbitals. The fractionally occupied orbitals in the FON-HF procedure were restricted to the CASCI active space as required by the implementation of the analytical gradients in TeraChem.32 The surface hopping calculations were performed with our in-house code ABIN.40 The FOMO-CASCI energies, gradients, and nonadiabatic coupling vectors were calculated at each step using the TeraChem program.33,34 These values were passed via a new message passing interface (MPI) code connecting ABIN (or any other MD program) and TeraChem. The interface is analogous to the one developed for QM/MM simulations with the Amber package.41,42 The initial guesses for the wave function, floating occupation numbers, and CI vectors were passed from the previous step. The state-averaged (SA) CASSCF energies, forces, and nonadiabatic couplings were calculated using the MOLPRO program, 43 which was connected with ABIN via a file-based interface. For performance reasons, the nonadiabatic coupling vectors were calculated only between states that were sufficiently close in energy44 (energy difference below 3 eV) and that had significant electronic populations above 0.001. During the simulations, the total energy conservation was monitored, and the simulation was stopped when energy was not conserved due to changes in the active space. We have verified that the procedure does not affect the statistical significance of the results presented here.

SE = z

1 p (t )(1 − pi (t )) N i

(8)

where N is the number of trajectories and pi is the population of state i at time t. We used z = 1.96 corresponding to the 95% confidence level. This formula is an approximate estimate for a confidence interval of a binomial (or multinomial) distribution, where we treat each trajectory as an independent trial at each time interval.

■

RESULTS AND DISCUSSION In this section, we compare the results of TSH-MD simulations with FOMO-CASCI and CASSCF wave functions. We test the FOMO-CASCI approach for the ethylene molecule, methaniminium cation, and malonaldehyde molecule. C


Article


Figure 2. Electronic populations (A) and structural evolution (B) of the ethylene molecule excited into the first excited state (S1) at t = 0. The populations were modeled with a surface hopping scheme using the SA-3-CASSCF(2,2) and FOMO-CASCI(2,2) wave functions with the 6-31G* basis set. Only the first two states were considered in the dynamical simulations. The broadening parameter β for FOMO-CASCI was set to 0.2 au.

Photodynamics of Ethylene: A Two-State Problem. Ethylene is an archetypal system of photoisomerization dynamics.45 It has been investigated many times with both theoretical31,46−53 and experimental46,47,54−58 techniques. It was shown that ethylene relaxes to the ground state through two types of conical intersections (CIs).47 The first of these is located near a twisted and pyramidalized structure, and the second involves an ethylidene structure (CH3CH) accessed after H migration. Experimentally, two main channels are observed on a longer time scale: (1) hydrogen molecule formation or (2) dissociation of two hydrogen atoms.59,60 In the present calculations, we have used an active space of two electrons distributed in two active orbitals and the 6-31G* basis set. Only the lowest two singlet states were considered in the dynamics simulations. However, to avoid root-flipping problems, the CASSCF calculations were performed using state averaging over three states, i.e., SA-3-CASSCF(2,2). Despite its simplicity, this active space is known to capture the essential features of ethylene photodynamics.4,45,46,61 The first valence excitation in ethylene is of ππ* character and is found experimentally at 7.7 eV.54,62 The CASSCF method significantly overestimates the first excitation energy, placing it around 10 eV,26 and this value is also reproduced by the FOMO-CASCI calculations. Nevertheless, the shape of the potential energy surfaces is similar to what would be obtained with a higher quality wave function.4,63 As discussed above, two major reaction coordinates dominate ethylene photodynamics in the first excited state: twisting of the CH2 group about the CC bond and pyramidalization of the CH2 group (see Figure 1A). After photoexcitation, the molecule rotates around the nominal double bond, and after

some time, pyramidalization begins at one of the C atoms. This leads to a point of degeneracy (conical intersection) between the excited and ground state. Here, the electronic population can be efficiently funneled into the vibrationally hot electronic ground state. Figure 1 shows a cut through the potential energy surfaces of the ground and first excited states along the two reaction coordinates. At first sight, there is no major difference between the energetics as described by CASSCF and FOMO-CASCI. The largest discrepancy is observed in the region with the pyramidalization angle larger than 60 degrees and the twisting angle around its ground state value. However, this region is rarely accessed by the molecule after excitation to the first excited state. In Figure 2A, the excited-state populations of the ethylene molecule calculated with CASSCF and FOMO-CASCI wave functions are depicted. There is essentially no difference between these two methods in terms of the lifetime of the S1 state. The S1 population decay can be reliably fitted by a delayed exponential function e−(t−t0)/τ (the population is considered to be 1 until t0). The fitted parameters t0 and τ are 24 and 46 fs for the CASSCF dynamics and 29 and 46 fs for FOMO-CASCI dynamics. The predicted mechanism of the deactivation is the same for both CASSCF and FOMO-CASCI methods. The dominant mechanism is twisting and pyramidalization, followed by population transfer. According to our SA-3-CASSCF(2,2) simulations, ∼63% of the ethylene molecules relax to the ground state via twisted-pyramidalized conical intersections and the remaining 37% via the ethylidene-like conical intersection that involves hydrogen transfer. The FOMO-CASCI simuD


Article


Figure 3B shows the branching among various reaction channels. Dissociation of the C−N bond takes place shortly after the system quenches into the ground state and represents the only significant reaction channel. Both the CASSCF and FOMO-CASCI methods predict that C−N bond rupture occurs ∼60% of the time. We can conclude that there are no significant differences between the CASSCF and FOMOCASCI results for this system. Malonaldehyde: Excited-State Proton Transfer. Malonaldehyde (MA) is one of the simplest molecules where proton transfer is observed in both the ground and excited states. Figure 4A shows the reactant and product in the excited-state proton transfer (which are indistinguishable in malonaldehyde but not in larger related molecules exhibiting excited-state intramolecular proton transfer).66,67 The excited-state proton transfer in MA is particularly interesting because of the role of a 3-state conical intersection in the photodynamics.68−70 The first two excitations in malonaldehyde are of nπ* and ππ* character. The calculated excitation energies are listed in Table 2. Although both the FOMO-CASCI and CASSCF methods predict the energy of the first excited state close to the reference CASPT2 calculation, the energy of the second excited state is significantly affected by dynamical correlation. The active space used in our calculations consists of four electrons in six orbitals, and the CASSCF orbitals were optimized for the three lowest states. Although a slightly smaller (4,4) active space with only four active orbitals was used in previous CASSCF studies,69,70 we found this active space unstable for FOMO-CASCI due to rotations between π* and σ* orbitals. The FON-HF equations have at least two different local minima for the (4,4) active space, and it was difficult to ensure convergence to one of these minima uniformly. In contrast, the active space (4,6) was remarkably stable, and 94 of 100 trajectories reached the full 300 fs simulation. With CASSCF using the (4,6) active space, we only succeeded in bringing 70 trajectories to the total simulation time. Our code does not implement the full set of safeguards to ensure active space orbital continuity that is normally used in ab initio multiple spawning dynamics61 so it is possible that this could be improved. However, the important point is that FOMO-CASCI is able to provide a continuous set of orbitals without any special care. We further discuss the stability of FOMO-CASCI and CASSCF trajectories below. For malonaldehyde, we tested two different broadening β parameters for the FOMO-CASCI calculations. On the basis of the comparison to the CASSCF results for the Franck−Condon point, β = 0.3 is a reasonable choice and was used in a previous study.26 In fact, the FOMO-CASCI method with this β value provides good energetics for different important points on the malonaldehyde PES (see Figure SI1). However, we found that photodynamical simulations with β set to 0.2 are in somewhat better agreement with respect to the CASSCF results. We can conclude that focusing on a few specific points of the PES is not always the best guiding principle while also stressing the potential for sensitivity of the results to fine details of the energetics. Figure 5 depicts the electronic population for the three electronic states of malonaldehyde. The MA molecule quenches rapidly from the S2 state (originally of ππ* character) into the S1 state (nπ* character in the Franck−Condon region). The S2 lifetimes from CASSCF and FOMO-CASCI are around 75 and 30−50 fs, respectively. The subsequent population depletion of the S1 state is slower for CASSCF and FOMO-CASCI with β =

lations provide similar results; the relaxation dynamics are split between the twisted-pyramidalized (75%) and ethylidene (25%) conical intersections. Figure 2B shows the reaction products formed within the simulation time. The main photochemical channel observed for the ethylene molecule is hydrogen atom transfer (the formation of CH3CH). In some 10% of our simulations, we observed the formation of molecular hydrogen. Other channels, like C−C bond breaking or single hydrogen dissociation were negligible for both methods. The CASSCF and FOMO-CASCI simulations thus provide essentially the same results within statistical accuracy. Methaniminium: A Three-State Problem. The methaniminium cation (CH2NH+2 , methyleniminium) is arguably the simplest model to investigate three state photodynamics. The ion also represents an example of a strongly polarized π bond. As such, it has served as a convenient model in photodynamical simulations.44,64,65 The first two excitations are of σπ* and ππ* character. The first excited state has very small oscillator strength, whereas the second excited state is optically bright. The excitation energies of methaniminium cation are shown in Table 1. The FOMO-CASCI method predicts slightly higher Table 1. Vertical Excitation Energies of the Methaniminium Cation Calculated Using the (12/8) Active Space with the 6-31G* Basis Set and the Structure Optimized at the MP2/ aug-cc-pVDZ Level energy/eV excitation character

CASSCF

FOMO-CASCI (β = 0.2 au)

MR-CISD+Qa

σπ* ππ*

9.17 10.06

9.27 10.34

8.35 9.17

a

MR-CISD+Q/SA-9-[CAS(6/4)+AUX(4)]/d-aug-cc-pVDZ level, taken from ref 64.

excitation energies than those of the CASSCF method. Dynamical correlation further shifts down these energies by ∼0.9 eV.64 The reaction coordinates and the potential energy surface shapes should be similar to the ethylene case. However, the twisting coordinate results in a S1/S0 crossing near the 90° angle without pyramidalization. The simulation protocol is similar to the one used above for ethylene, but the molecule was initially excited to the bright S2 state. An active space consisting of 12 electrons in 8 orbitals was used.44 The CASSCF orbitals were state-averaged over the three lowest states. The broadening parameter β = 0.2 au was used for FOMO-CASCI simulations. We note that, in the CASSCF dynamics, around 20% of the trajectories ended prematurely due to problems with energy conservation. Similar issues with the same active space were encountered in previous TSH-MD studies.44,64,65 Upon irradiation, the methaniminium molecule quickly quenches into the S1 state (Figure 3A). This ultrafast quenching cannot be reliably fitted to an exponential function. For both FOMO-CASCI and CASSCF, the S2 population drops below 20% within 20 fs. This is in agreement with previous surface hopping results with the MR-CISD method.64 The population of the S1 state peaks around 10 fs for both the CASSCF and FOMO-CASCI electronic wave functions. Even the long-time dynamics are quite similar for both approaches. E


Article


Figure 3. Electronic populations (A) and structural evolution (B) of the methaniminium cation initially excited into the second excited state (S2). The populations were modeled with a surface hopping scheme using the SA-3-CASSCF(12,8) and FOMO-CASCI(12,8) wave functions with the 6-31G* basis set. The broadening parameter β for FOMO-CASCI was set to 0.2 au.

Figure 4. (A) Two equivalent malonaldehyde structures illustrating the proton transfer. (B) Important points on the excited potential energy surface of malonaldehyde: hydrogen transfer intersection (HTI) between S2 and S1 states, 3-state intersection (3SI), and minimum conical intersection (MECI) between S1/S0 states.

simulations performed here are consistent with previous ab initio multiple spawning (AIMS) calculations.69,70 On the other hand, much faster depopulation of the S1 state is seen for FOMO-CASCI simulations with β = 0.3. Depopulation of the S2 state proceeds via two competing channels: either at geometries close to the hydrogen transfer intersection (HTI) or via a conical intersection with a twisted CC bond (see Figure 4b). The twisted conical intersection is in fact a three-state conical intersection connecting the S2, S1, and S0 states69 (see Figure SI1). The HTI channel is the dominant one even though energetically it lies well above the twisted 3-state conical intersection. This finding is in line with previous simulations.69 However, the twisted conical intersection plays a much more important role in the CASSCF simulations, where it is responsible for the deactivation of

Table 2. Vertical Excitation Energies of Malonaldehyde Calculated Using the (4/6) Active Space with the 6-31G* Basis Set and the Structure Optimized at the MP2/6-31++G** Levela energy/eV excitation character

CASSCF



MSCASPT2

nπ* ππ*

4.03 6.04

4.72 6.21

4.23 6.24

3.55 4.51

a

A reference calculation was performed at the MS-CASPT2 level using the (4/6) active space with the aug-cc-pVDZ basis set.

0.2, and both methods end up with a population of 50% still residing in the first excited state at t = 300 fs. The CASSCF F


Article


Figure 5. Electronic populations for the malonaldehyde molecule excited into the second excited (S2) state at t = 0. The populations were modeled within a surface hopping scheme using the CASSCF and FOMO-CASCI methods with the (4/6) active space averaged over three states using the 6-31G* basis set. The broadening parameter β for FOMO-CASCI was set to 0.2 or 0.3 au.

Figure 6. Fraction of trajectories with a transferred hydrogen for each electronic state of malonaldehyde. The populations were calculated within the surface hopping scheme using the SA-3-CASSCF and FOMO-CASCI methods with the (4/6) active space averaged over three states using the 6-31G* basis set. The broadening parameter β for the FOMO-CASCI method was set to 0.2 or 0.3 au.

G


Article


Figure 7. Comparison of the energy conservation for the SA-CASSCF and FOMO-CASCI methods. Distributions of maximum total energy differences along surface hopping trajectories for ethene (left panels), methaniminium cation (middle pannels), and malonaldehyde (right panels). See the text for the simulation details.

Energy Conservation During Surface Hopping MD. Next, we discuss some more technical aspects of the FOMOCASCI and CASSCF surface hopping simulations. As discussed in the Introduction, the FOMO-CASCI method is expected to be more stable as it decouples the optimization of orbital coefficients from the configuration interaction coefficients. This inherent stability is nicely illustrated by comparing the energy conservation during the FOMO-CASCI- and CASSCF-based surface hopping trajectories. Figure 7 shows a histogram of the maximum total energy “jump” that took place between the consecutive time steps along each trajectory for each method and molecule. Using this metric, the majority of FOMO-CASCI simulations have energy conservation better than 0.05 eV, and energy jumps greater than 0.2 eV are not rare for the CASSCF method using the same active space, as implemented in our ABIN code interfaced to MolPro. The occurrence of large energy jumps in CASSCF is in part because our code does not take special care to follow a continuous set of active orbitals throughout the dynamics. We do use the active space orbitals from the previous time step as an initial guess but otherwise rely on the convergence algorithm in MolPro to provide a continuous set of active orbitals. The importance of this comparison between FOMO-CASCI and CASSCF is that it shows one does not need to take special care with orbital continuity when using FOMO-CASCI (and not because it provides any indication of the energy discontinuities, which will be obtained with a carefully implemented CASSCF dynamics code). Finally, we should note that energy jumps in CASSCF do not necessarily affect the underlying dynamics of the

approximately 25% of trajectories. For the FOMO-CASCI trajectories with both broadening parameters, the population transfer takes place almost exclusively via the HTI intersection. This probably explains a faster depopulation of the S2 state on the FOMO-CASCI potential energy surface because the HTI intersection is closer to the Franck−Condon point. This discrepancy may be caused by a slightly higher barrier along the reaction coordinate connecting the Franck−Condon point and the twisted conical intersection. The barrier is almost negligible (

Non-Adiabatic Ab Initio Molecular Dynamics with Floating Occupation

Recommend Documents