Diabatic States at Construction (DAC) through Generalized Singular

Oct 2, 2018 - ... is presented to transform adiabatic potential energy surfaces into a diabatic representation by generalized singular value decomposi...
0 downloads 0 Views 2MB Size
Letter Cite This: J. Phys. Chem. Lett. 2018, 9, 6038−6046

pubs.acs.org/JPCL

Diabatic States at Construction (DAC) through Generalized Singular Value Decomposition Meiyi Liu,† Xin Chen,† Adam Grofe,*,† and Jiali Gao*,†,‡ †

Laboratory of Theoretical and Computational Chemistry, Institute of Theoretical Chemistry, Jilin University, Changchun, Jilin Province 130023, China ‡ Department of Chemistry and Supercomputing Institute, University of Minnesota, Minneapolis, Minnesota 55455, United States

Downloaded via UNIV OF SUNDERLAND on October 4, 2018 at 22:30:04 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: A procedure, called generalized diabatic-at-construction (GDAC), is presented to transform adiabatic potential energy surfaces into a diabatic representation by generalized singular value decomposition. First, we use a set of localized, valence bond-like configuration state functions, called DAC, as the basis states. Then, the adiabatic ground and relevant excited states are determined using multistate density functional theory (MSDFT). GDAC differs in the opposite direction from traditional approaches based on adiabatic-to-diabatic transformation with certain property restraints. The method is illustrated with applications to a model first-order bond dissociation reaction of CH3OCH2Cl polarized by a solvent molecule, the ground- and first-excited-state potential energy surfaces near the minimum conical intersection for the ammonia dimer photodissociation, and the multiple avoided curve crossings in the dissociation of lithium hydride. The GDAC diabatization method may be useful for defining charge-localized states in studies of electron transfer and proton-coupled electron transfer reactions in proteins.

D

by comparison of the resulting adiabatic energies with experimental observables. The VB configurations that correspond to the asymptotic dissociation limits, whose valence characters are strictly maintained at all nuclear coordinates, are diabatic states by construction. The VB-like configurations in the DAC approach can in principle be used directly in dynamic simulations;7,39 however, they are typically nonorthogonal and span the entire active space in a multiconfiguration self-consistent-field (MCSCF) calculation such as ab initio VB (equivalent to complete active space SCF, CASSCF) and the present MSDFT.40 In practice, only a limited number of primary states (NP) are of interest in a photochemical reaction, involving mainly the ground and the first few excited states. In this case, it is convenient (but not essential) to employ the same number of equivalent diabatic states in dynamic calculations. This also gives rise to the chemically intuitive feature that two diabatic states cross when the adiabatic states interchange their valence characters, resulting in an avoided crossing or a conical intersection. In this Letter, we present a diabatization procedure, making use of generalized singular value decomposition (GSVD) to represent the adiabatic PES. The method is called the generalized diabatic-at-construction (GDAC), and it does not require additional computational efforts such as calculations of the

iabatic representation of potential energy surfaces (PESs) is widely used in studies of charge transfer and nonadiabatic dynamic simulations of photochemical reactions involving conical intersections or locally avoided crossings.1 It provides diabatic coupling elements of the Hamiltonian operator, which are scalar and smoothly varying with nuclear coordinates.2−9 On the other hand, the nonadiabatic coupling or derivative coupling vectors of the adiabatic PES are not smooth and exhibit singularities in regions near conical intersection.10,11 Consequently, numerous techniques have been developed for constructing diabatic states.4,8,9,12−36 However, diabatic states are not unique,37 and most methods involve an orthogonal transformation of the adiabatic surfaces with restraints by desired properties or minimization of the nonadiabatic coupling vectors. These approaches, known as diabatization methods, may be classified as the adiabatic-todiabatic (ATD) category because the adiabatic PESs are determined first and the diabatic states that are constructed through an ATD transformation are admixtures of the corresponding valence bond (VB) configurations of the asymptotic dissociation products. Recently, we described an alternative approach, called diabatic-at-construction (DAC),38 which starts with a set of VB-like configuration state functions as the basis states, based on which the adiabatic ground and excited surfaces are obtained through nonorthogonal configuration interaction in multistate density functional theory (MSDFT). The validity and the required number of VB-like configurations for the adiabatic states of interest can be verified © XXXX American Chemical Society

Received: August 11, 2018 Accepted: October 2, 2018 Published: October 2, 2018 6038

DOI: 10.1021/acs.jpclett.8b02472 J. Phys. Chem. Lett. 2018, 9, 6038−6046

Letter

The Journal of Physical Chemistry Letters nonadiabatic coupling vectors as an optimization restraint. The resulting GDAC states correspond to the closest representation of the adiabatic PES relative to the valence characters for the diabatic statesthe asymptotic dissociation limits for the reactions of interestin terms of the Frobenius norm. We illustrate this approach by three representative examples, including a simple two-state model in an SN1-like dissociation involving covalent and ionic curve crossing, potential energy curves of LiH with multiple states, and a conical intersection in photodissociation of an ammonia dimer. MSDFT is used to determine the ground and excited PESs for the illustrative cases.40,41 MSDFT follows a dynamic-thenstatic ansatz to treat electron correlation. First, (1) dynamic correlation is incorporated via block-localized Kohn−Sham density functional theory (BLKS-DFT) into the individual configurations in the active space.41 Then, (2) static correlation is treated by nonorthogonal configuration interaction, important for systems with strong correlation and neardegeneracy.40 Consequently, it does not require perturbation calculations to yield accurate results for ground and excited states. Furthermore, because dynamic correlation is included in the first place, the orbital space is also optimized such that MSDFT generally requires far fewer configurations than CASSCF. 40 Theoretical and computational details of MSDFT can be found in previous studies.40,42,43 The goal of this study is to seek a representation, i.e., diabatization, of NP adiabatic states {|ΘI⟩} in terms of the configurations corresponding to their asymptotic limits on the PES. This is accomplished by generalized SVD. Generalized Diabatic-at-Construction. The present GDAC method requires the coefficient matrices C and D of the NP adiabatic and diabatic reference states as input and gives the transformation matrix Ũ that defines the diabatic states. We first define and describe these matrices. Then, the GDAC method is formulated and solved as a mathematical problem of least-squares minimization. The wave functions and potential energies of the adiabatic ground and excited states are obtained from MSDFT calculations40,41

diabatic states. The procedure is called GDAC to distinguish the use of VB states directly, called DAC.38 Most ATD transformations,12−34 whether orbital- and property-dependent or -independent, direct or indirect, focus on finding the transformation X to obtain diabatic energies, but in GDAC, we emphasize optimization of the wave functions of diabatic states such that Φds will have the maximum resemblance of the valence characters that the diabatic states represent. The transformation matrix X becomes a natural consequence of the optimized Φds. This reflects the difference of the present DAC strategy (diabatic first followed by adiabatic),7,38 opposite from that of ATD methods (adiabatic first). The optimized diabatic states Φds are given in the MSDFT basis states by Φds = Ũ TΨ

where Ũ is the coefficient (transformation) matrix of the NP diabatic states expressed in the spin-adapted VB basis configurations Ψ (eq 1). A key ingredient in the GDAC method is to maximally maintain the “valence characters” of the diabatic states. We define the diabatic reference states. Diabatic reference states (D) are configurations defining the target valence characters that the diabatic states represent (see the Supporting Information). They are defined by the asymptotic dissociation limit of a photochemical reaction, denoted by |ΨDI (r → ∞)⟩, where r is a set of internal coordinates and r → ∞ is meant to indicate the fully separated reaction products in various electronic states {I = 1, ..., NP}. |Ψ ID(r → ∞)⟩ is an eigenfunction of the electronic Hamiltonian, only identical to the adiabatic states |ΘI⟩ at the equilibrium geometry of the final photochemical reaction products. Thus, the adiabatic states in the primary space are related to the diabatic reference states by {|ΘI⟩ → |ΨDI (r → ∞)⟩}. {|ΨDI (r → ∞)⟩} have specific localized (i.e., diabatic) valence characters, for example, biradical or ion-pair, π → π*, n → π*, etc. They can always be expressed in terms of Lewis structures, i.e., VB-like configurations, which could be just one of the DAC configurations in the active space of MSDFT or a combination of several such configurations, given by ΨD = DTΨ, to account for hybridization or electronic excitation at arbitrary coordinates of r.38 In comparison with standard stateaveraged CASSCF calculations, diabatic reference states correspond to adiabatic states at equilibrium geometries of the reaction products, but unlike CASSCF states, which become mixed with other states, diabatic reference states maintain their valence characters away from their equilibrium geometries. Diabatic reference states ΨD are essentially localized VB states at different geometries, and they are not associated with localized orbitals or any specific physical properties. Having defined diabatic reference states, given by the coefficient matrix D, the optimization of Φds = Ũ TΨ is translated into a purely mathematical problem of least-squares minimization, which can be easily solved by singular value decomposition (SVD). Given the coefficient matrices C and D for the adiabatic states and diabatic reference states, and their SVDs C = Ũ σṼ T and D = Q0sRT0 , both of rank NP, we seek to find the (diabatic) representation Ũ of C that minimizes the Frobenius norm ||Ũ − Q0||F2

NF

|ΘI (r)⟩ =

∑ CAI|ΨA(r)⟩ A=1

(1)

or, in matrix form Θad = CTΨ

(2)

where Θad = {|ΘI⟩; I = 1, ..., NP} and Ψ = {|ΨA⟩; A = 1, ..., NF} are, respectively, the adiabatic states and spin-adapted VB-like basis configurations arranged as column vectors, C is the NF × NP coefficient matrix with elements defined in eq 1, and NP and NF are the number of adiabatic states of interest (called primary states) and basis configurations (active space). Throughout this Letter, the symbols I, J, ... specify adiabatic states and A, B, ... denote diabatic or basis state functions. In the present GDAC method, we optimize the representation of the NP adiabatic states Θad by an equivalent set of diabatic states Φds that maximally maintain the valence characters of the photochemical reaction products Θad = XTΦds

(4)

(3)

where Φds denotes the GDAC diabatic states (ds) {|ΦB⟩; B = 1, ..., NP}, and X is an NP × NP unitary matrix that formally transforms the adiabatic states into an equivalent set of NP 6039

DOI: 10.1021/acs.jpclett.8b02472 J. Phys. Chem. Lett. 2018, 9, 6038−6046

Letter

The Journal of Physical Chemistry Letters Ũ = min ∥U − Q 0∥F 2 U C = Uσ VT

X ≡ σ Ṽ T = Q 0 TC (5)

However, in practice, it is more useful to transform the Hamiltonian matrix Hds in the orthogonal diabatic basis. Then, the diagonal elements of Hds correspond to the energies of the GDAC diabatic states, and the off-diagonal elements are diabatic couplings. Diagonalization of Hds yields eigenvalues that will be identical to the energies of the adiabatic states and the transformation matrix X, being the matrix of the eigenvectors. The GDAC computational procedure amounts to a simple orthogonal transformation, which may be summarized in the following three steps:

This amounts to a special case of a low-rank matrix approximation to an observed (target) matrix in that the ranks of both the “approximate” and target matrices are the same, e.g., the representation is exact, ensuring that the adiabatic potential energies are exactly reproduced in terms of the diabatic representation. This type of optimization has a wide range of applications, such as spectral decomposition, principal component analysis, and genome-scale analysis for cancer screening.44,45 Detailed solution and discussion of the minimization of eq 5 are given in the Supporting Information. Here, we present the main results. SVD is not unique (just as unitary transformation of molecular orbitals). However, there is a unique decomposition (Ũ σṼ T) that minimizes ||Ũ − Q0||F2. This implies that such a SVD is constrained, called generalized SVD (GSVD),46 denoted by the tilde over the decomposition matrices. In regular SVD, D = Q0sRT0 , both the left and right decomposition matrices are orthonormal (unitary in full dimension) (e.g., QT0 Q0 = RT0 R0 = I). In GSVD, C = Ũ σṼ T, they are orthonormal under the constraints W and M

Ṽ TWṼ = I

(6a)

Ũ TMŨ = I

(6b)

(1) Define D, the diabatic reference states that the diabatic states represent, and perform SVD to yield Q0. (2) Form PC = C(CTC)−1CT, and compute Ũ = PCQ0. Gram-Schmidt orthogonalize Ũ → Ũ ⊥, if the overlap ⟨Ψ|Ψ⟩ ≠ I, to yield an orthonormal Φds = Ũ T⊥Ψ. (3) Transform the MSDFT Hamiltonian matrix in the full NF basis of the active space into the NP GDAC diabatic basis, Hds = Ũ T⊥HMSDFTŨ ⊥, which yields the diabatic state energies (diagonal elements) and diabatic couplings (off-diagonal). If desired, the small NP × NP matrix Hds can be diagonalized to give the adiabatic state energies for validation of results.

Therefore, the minimization of the Frobenius norm ||Ũ − Q0||F2 in eq 5 is subjected to the constraints of eq 6. Rewriting the SVD as Ũ = CWṼ σ−1, we form the optimization object function h(W,λ) = ||CWṼ σ−1 − Q0||F2 + λ||Ṽ TWṼ − I||F2, where λ is a Lagrangian multiplier. Taking the derivatives h(W,λ)/∂W = 2CT(CWṼ σ−1 − Q0)σ−1Ṽ T + 2λṼ (Ṽ TWṼ − I) Ṽ T (see the Supporting Information)47 and solving for ∂h(W,λ)/∂W = 0, we find that the second term vanishes because of eq 6a, and CTCWṼ σ−1 − CTQ0 = 0 after multiplying WṼ σ from the right. Rearranging this result, we find WṼ σ −1 = (CTC)−1CTQ 0

Notice that it is not necessary to actually construct the constraints W and M in computation because analytical minimization results have been obtained via GSVD (eq 8). We note that the present GDAC method has some similarities at first glance to the block diagonalization approach introduced by Cederbaum and co-workers in the 1980s. 12−15,48 Mathematically, both methods are formulated as a minimization problem. In block diagonalization, the optimization focuses on the object min||X − 1||F2,13 the adiabatic−diabatic transformation, whereas in GDAC, we optimize the diabatic wave function (its coefficient matrix Ũ ) to be closest to the valence characters Q0, ||Ũ − Q0||F2. Because orthogonal projection has eigenvalues of 1, or rank NP, the minimum ||X − 1||F2 = 0 is ensured as a result. In practice, as Domcke and co-workers pointed out,14,15 it is typically necessary to transform adiabatic states through optimization of “diabatic” molecular orbitals.48 In GDAC, we do not work in orbital space, and the diabatic reference states are VB configuration state functions. Multistate density functional theory is used to obtain the adiabatic ground and excited state energies by diagonalizing the nonorthogonal configuration interaction Hamiltonian to yield the coefficient matrix C (eqs 1 and 2).40,41 We note that different VB-like basis configurations |ΨA(r)⟩ in MSDFT do not share a common set of orbitals because they are separately optimized and nonorthogonal. They are “VB-like” because we use one or a minimal number of spin-adapted block-localized determinants to represent a given Lewis structure, rather than using all possible VB configurations as in ab initio VB (which is equivalent in the nonorthogonal basis to CASSCF with orthogonal orbitals for a given active space). Each VB-like determinant is treated as a KS-DFT determinant, except that BLKS orbitals are used. The block-localization scheme effectively treats the full valence-bond resonance delocalization (or a local CAS) involving numerous determinants by a single BLKS configuration state function. Consequently, the number of configurations is significantly reduced locally in MSDFT

(7)

Then, we obtain the optimized coefficient matrix (eq 4) for the diabatic states Ũ = CWṼ σ −1 = C(CTC)−1CTQ 0 = PCQ 0 T

−1

(9)

(8)

T

where PC = C(C C) C is the orthogonal projection onto the adiabatic space. Recall that orthogonal projection is both idempotent (PC2 = PC) and symmetric (PTC = PC). Because the DAC basis states are generally nonorthogonal, the orthogonal projection has the form P̂ C = C(CTSC)−1CT, where S is the overlap matrix of the MSDFT basis states. In this case, the generalized idempotency condition is satisfied.41 If the adiabatic basis states are orthogonal, P̂ C is simply reduced to CCT. It is important to note that Ũ itself is in fact orthonormal with M = P C (eq 6): Ũ T Ũ = (P C Q 0 ) T (P C Q 0 ) = (PCQ0)T(PCPCQ0) = Ũ TMŨ = I. Equation 8 has the simple physical interpretation that the GDAC diabatic representation (Φds = Ũ TΨ) of the adiabatic states is a projection of the valence characters in the adiabatic space. This is the main result of this work. It is straightforward to obtain the adiabatic ↔ diabatic transformation matrix (see the Supporting Information) as follows 6040

DOI: 10.1021/acs.jpclett.8b02472 J. Phys. Chem. Lett. 2018, 9, 6038−6046

Letter

The Journal of Physical Chemistry Letters

to a determinant energy-weighted average of Kohn−Sham correlation for the two interacting states. Alternatively, an overlap integral weighting scheme can also be used.41,43 In most applications, both yield similar results, though eq 10 has an advantage for interactions with vanishing overlap integrals (for example, at distant separations between two fragments). Note that both exchange and static correlation (as a result of the CI) are partially included in the first term in eq 12, whereas only dynamic correlation (from the external space outside of the active space used in MSDFT) is present in the second term (eq 13).40,52 In the present context, VTDF AB [ΨA,ΨB,ρAB(r)] has the physical interpretation of dynamic correlation contribution to the diabatic couplings between diabatic states A and B. Illustrative examples of the GDAC method are presented below. We begin with a two-state model, exhibiting dual curve crossing in the C−Cl bond dissociation of H2O···CH3OCH2Cl in the gas phase. It also corresponds to a first-order nucleophilic substitution mechanism (SN1) in water. Because of the strongly electron-donating effect of the methoxy group along with solvent−molecule polarization, it is sufficient to lower the ionic state energy below the covalent curve in the intermediate distance between C and Cl, exhibiting two avoided curve crossings to level off as diradical products in the ground state. The PESs for the covalent ground and ionic states as a function of the C−Cl distance are determined using MSDFT with the PBE0 functional along with the aug-cc-pVDZ basis set. The main purpose here is to show the qualitative representation of the adiabatic PES by two diabatic states via a SVD with constraint. Of course, there is a simple, but arbitrary, rotation to make a 2 × 2 transformation in a twostate model,53,54 and the nonorthogonal VB basis states themselves constitute a “pure” diabatic representation by construction (although they are nonorthogonal). The GSVD yields a diabatic transformation of mixed states that have the closest resemblance to the VB states in terms of the Frobenius distance. A water molecule is located at 2.5 Å from the central carbon atom along the C−Cl bond direction of the substrate (MeO)CH2Cl, which is aligned with the z-axis. Then, the molecular complex is partitioned into three orbital blocks, respectively, according to the fragment (H2O), p orbitals of Cl, and all basis orbitals of CH3OCH2 plus the remainder basis orbitals on Cl. BLKS orbitals are written as linear combinations of the atomic orbital basis functions only from each individual orbital block. As a result, the total charge and electronic spin of each fragment block can be strictly constrained locally. With the above partition, we define three spin-adapted covalent configurations for the three sets of 2p2q and 3p1q basis orbitals and electrons of chlorine, where q = x, y, z, each in terms of two Slater determinants

relative to that for a full CAS space, with the trade-off for a “mean-field” representation of the local correlation. The latter (electron correlation within nondegenerate local fragments) is not a problem for KS-DFT to model systems with strong correlation. In fact, this is mostly compensated, if not completely, by including dynamic correlation via BLKS-DFT, resulting in improved accuracy over ab initio VB (or equivalently, CASSCF) that lacks dynamic correlation. The electron densities for the adiabatic ground and excited states are given by NP

ρI (r) =

∑ {CAI 2ρAms (r) + ∑ CAICBIρAB (r)} A=1

B≠A

(10)

where ρI(r) is the charge density of adiabatic state I and ρms A (r) and ρAB(r) are, respectively, the charge density for BLKS determinant A and the transition density involving configurations A and B.41 The MSDFT energy functional is given by40 NF

EIms[ρI (r)] =

∑ {CA 2EBLKS[ρAms (r)] A=1

∑ CAICBIETDF[ΨA, ΨB, ρAB (r)]}

+

B≠A BLKS

(11)

[ρms A (r)]

where E denotes the energy of configuration A, determined by using any standard Kohn−Sham exchangems correlation functional EKS xc [ρA (r)] with BLKS orbitals (PBE0 is used here throughout) and ETDF[ΨA,ΨB,ρAB(r)] is the transition density functional (TDF), a novel class of density functional that does not exist in the Kohn−Sham approximation of DFT.40,41 ETDF[ΨA,ΨB,ρAB(r)] is a functional of the configuration state functions |ΨA(r)⟩ and |ΨB(r)⟩ and the associated transition density. Its exact expression is not yet known, and it is approximated in applications, but under some special situations, its value can be exactly obtained, given a particular KS-DFT approximation used to describe the two determinant configurations A and B,42 such as the spin coupling of two unpaired electrons to yield a pair of singlet and triplet states. It would be of interest to develop a general purpose TDF in MSDFT in concert with the hybrid functional used in the BLKS-DFT for the diagonal energies. In the present application, the off-diagonal matrix element of the Hamiltonian that includes contributions from the transition density functional ETDF is computed by40,41 SD TDF HAB ≡ ETDF[ΨA , ΨB, ρAB (r)] = H AB + VAB

(12)

where H is the electronic Hamiltonian, and HSD AB = ⟨ΨA|H|ΨB⟩ is a nonorthogonal matrix element computed using BLKS Slater determinants (SD).39,49,50 The second term in eq 12 is the correlation energy due to the transition density and the BLKS energies and determinants,40,41,43,51 and it is determined by TDF VAB



SD H AB

EAHF

+

(EcBLKS[ρAms (r)] HF EB

1 |Ψq ,cov⟩ = Nq ,cov[Â {(Ω H2O)(ΩCφCCl )(φp̅ 3)} q



+

EcBLKS[ρBms (r)])

1 3 Â {(Ω H2O)(ΩCφCCl ̅ )(φpq)}]

(14)

In eq 14, Nq,cov is a normalization constant, the parentheses () specify fragment blocks, and the symbols ΩH2O and ΩC denote products of doubly occupied BLKS orbitals for water and the substrate fragment CH3OCH2Cl, except the pq orbitals on chlorine. φ1CCl specifies a fragment orbital on the methylene carbon along the z-axis pointing toward Cl, and φ3pq are

(13)

where, EHF and EBHF are the Hartree−Fock determinant A energies obtained by using the corresponding BLKS orbitals BLKS ms and EBLKS [ρms [ρB (r)] are DFT correlation c A (r)] and Ec energies for configurations A and B. Equation 13 corresponds 6041

DOI: 10.1021/acs.jpclett.8b02472 J. Phys. Chem. Lett. 2018, 9, 6038−6046

Letter

The Journal of Physical Chemistry Letters

ionic VB states essentially coincide with the adiabatic ground and excited states, respectively. The GSVD diabatic states show similar features as that of the pure VB states, but the ionic state has a quicker rise as a result of orthogonalization. Because of the large diabatic coupling at 3 Å, only the adiabatic PES is relevant, and for all practical purposes, only the outer avoided crossing is important for nonadiabatic dynamics. In this regard, the diabatic states and dual curve crossing in Figure 4 cannot be strictly considered as a realistic molecular analogue of the Tully model. Conical Intersection in (NH3)2. We next examine the PES associated with photodissociation of the ammonia dimer.57−60 Here, a discussion of the full PESs and the excited-state dynamics is beyond the scope of the present work; we focus on the “diabatization” of the adiabatic PES near the conical intersection along the H2N···HNH3 diradical separation pathway on the first excited state.59 For convenience, we have restricted the hydrogen atom that is transferred from the initially excited state (planar geometry) H2NH···NH3(Ã ) to the second ammonia (pyramidal) along the axis of the two heavy atoms, R(N−N).59 The full system is characterized by 14 spin-adapted VB, i.e., DAC, configurations from 21 determinants (Supporting Information), sufficient to describe both the adiabatic ground state and the first excited state. The active space of MSDFT is analogous to that of the NH3 monomer photodissociation, also using aug-cc-pVTZ and the PBE0 functionals,40 but the dissociating hydrogen atom is grouped into the second ammonia fragment. We note that if CASSCF was used, it would have required an active space involving 16 electrons (3 covalent bonds and 1 lone pair for each ammonia), consisting of a much larger number of determinants. MSDFT, as a density functional approximation, focuses on a minimal, valence-state representation that is directly associated with the adiabatic states of interesthere, the ground state and a valence-excited state. It is possible to obtain accurate results with a smaller number of configurations in the active space for two reasons. First, because the contributions from the vast majority of configurations in the corresponding CAS space are of dynamic correlation in nature as far as the two lowest adiabatic states are concerned, they are already included and can be successfully treated by blocklocalized KS-DFT (the dynamic-then-static ansatz). On the other hand, it is difficult to select or exclude any determinants in the CAS because of the use of delocalized molecular orbitals. Second, separate optimizations of BLKS orbitals in different individual (VB-like) basis configurations in MSDFT provide a fully compact (contracted) active space compared to a procedure employing a common set of orbitals such as CASSCF. In the latter case, a very large number of determinants must be used to effectively achieve the same effect of separate orbital optimization for a particular valence character important to the adiabatic PESs. At R(N−N) = 4.232 Å along the NH2 and NH4 diradical dissociation path on the excited state (Figure 2), a conical intersection between the ground (S0) and excited (S1) state surfaces is obtained, which corresponds to the minimum conical intersection point of the full surface at an energy of 4.83 eV above two separate ammonia molecules in the ground state. For comparison, an energy of about 4.9 eV has been estimated at the conical intersection at an N−N separation of ∼4.1 Å from state-averaged CASSCF/MRCI calculations with 16 electrons in 12 orbitals.58,59 In fact, all critical points on both the ground and excited PES for the ammonia dimer

products of orbitals of 2pq and 3pq character, with the latter having a single-electron occupation. Pairs of orbitals with and without a hat indicate α and β spin pairings to give an overall singlet configuration. The ionic configuration is represented by one determinant 0 |Ψion(C+Cl−)⟩ = Nion[Â {(Ω H2O)(ΩCφCCl )(Cl−)}]

(15)

Notice that we have kept all orbitals of the water molecule in a separate block without mixing with the substrate; there is no charge transfer nor bonding participation from the solvent molecule. Thus, the adiabatic potential energy profile follows strictly an SN1 (or covalent bond dissociation) mechanism. In principle, we could include a covalent configuration that combines the substrate fragment and the water molecule: Â {(Ω+H2OC)(Cl−)}. This, however, corresponds to a standard SN2 mechanism, different from the ionic and covalent curve crossing induced by solvation in the SN1 mechanism. The adiabatic and diabatic (VB and orthogonal) potential energy curves are shown in Figure 1. In the present model, the

Figure 1. Computed potential energy curves for the adiabatic ground (V1 in black) and excited (V2 in purple) states, covalent (dashed blue) and ionic (dashed yellow) VB configurations, and GDAC diabatic states with covalent (red) and ionic (green) dominant characters. The PBE0 functional with the aug-cc-pVDZ basis set is used in all calculations.

covalent and ionic diabatic states in dashed curves exhibit a surprising double curve crossing feature, at about 3.0 and 6.6 Å, respectively. The structures at different C−Cl distances were optimized using B3LYP/aug-cc-pVDZ for the adiabatic ground state. In the intermediate region between 3.0 and 6.6 Å, the energy of the ionic configuration dips below that of the covalent state, while the covalent configuration makes a dominant contribution to C−Cl bonding at shorter bond lengths and represents the asymptotic dissociation product. The dual crossing reaction path is in fact analogous to one of the most widely tested models proposed by Tully55 for nonadiabatic dynamics, except that the couplings are relatively weak and symmetric in that model.55,56 For CH3OCH2Cl, the electronic coupling at 3.0 Å is about 38 kcal/mol, strong enough to generate an energy gap of 74 kcal/mol between the adiabatic ground and excited states (Figure 1). In contrast, the coupling at 6.6 Å is just 0.1 kcal/mol, producing a weakly avoided curve crossing, switching, again, the valence character of the adiabatic ground state from ionic to covalent at longer C−Cl separations. Here, the effect of a single solvent molecule is not sufficient to fully lower the energy of the ionic state below that of covalent bond dissociation as a true SN1 process. Beyond the outer avoided curve crossing, the covalent and 6042

DOI: 10.1021/acs.jpclett.8b02472 J. Phys. Chem. Lett. 2018, 9, 6038−6046

Letter

The Journal of Physical Chemistry Letters

Figure 2. Ground and excited PESs for ammonia dimer about the conical intersection point along the minimum photodissociation path (N−N distance) and the H2N out-of-plane bending coordinates to yield the radicals H2N* and H4N* on the first excited state, and H2N− and H4N+ ion pair on the ground state. MSDFT is used with the PBE0 functional and the aug-cc-pVTZ basis set.

system from MSDFT are within 0.2 eV of experimental values or best estimates of computation.57−59 By symmetry, the branching space that lifts the energy degeneracy at the conical intersection corresponds to the largeamplitude interfragment N−N stretch mode and H2N out-ofplane bending. The PESs of the adiabatic ground and first excited states for the ammonia dimer near the intersection point are shown in Figure 2, whereas the diabatic states and coupling are depicted in Figure 3, using the ionic and excitedstate configurations as the diabatic reference states. Figure 3a shows the smooth crossing between the diabatic surfaces from GDAC decomposition, and their coupling depicted in Figure 3b fully splits the two surfaces to produce avoided crossings at geometries in the branching space except for one point, giving rise to a conical intersection in Figure 2. Similar to the findings in the monomer dissociation,40,61,62 the ground state of the NH2 product (2B1) correlates diabatically with the excited state of a locally excited NH3 in the dimer, whereas the excited state of the NH2 product (2A1) correlates diabatically with the ground state of the dimer complex, which has dominantly ionic character, close to the intersection point. Near the conical interaction, the excited-state surface (1A′) correlates diabatically with interfragment electron transfer in the ionic state of H2N−···NH+4 , and a proton transfer process occurs on the ground-state surface to yield two ammonia molecules as the (non-productive) photochemical reaction products. LiH. Finally, we present a GDAC representation of the four lowest 1Σ+ states of LiH as a function of the interatomic distance, R. Although LiH constitutes one of the simplest covalent bonds, the description of its potential energy curves is exceptionally challenging because multiple curve crossings are involved, spanning the entire interatomic separation. In the present discussion, it provides an interesting example to model multiple adiabatic states by a diabatic representation. The dissociation curves of LiH have been intensively inves-

Figure 3. Computed (a) diabatic states in ammonia dimer photodissociation as a function of the N−N separation and the H2N out-of-plane bending coordinates and (b) the associated diabatic coupling surface. MSDFT is used with the PBE0 functional and the aug-cc-pVTZ basis set. Notice that for a two-state system, the adiabatic state energies depend on only the squared diabatic coupling.

tigated,63−69 and we have shown that it is sufficient to use seven VB basis states in MSDFT to describe quantitatively the energies of the four lowest singlet states in comparison with experiments and MS-CASPT2 calculations.38 The seven VB configurations, in order of increasing energy (Supporting Information), except the first, include the |Ψion[(Li+)(H−)]⟩ ionic structure, the dissociation limits of the Li atom in the 2s, 2p, 3s, and 3p states, along with an excited configuration for the H(2s) atom and the reverse ionic configuration |Ψ−[(Li−)(H+)]⟩. The aug-cc-pVTZ basis set70 was used with the exception that the standard 3s-diffuse function for Li was replaced by the diffuse basis function in the Dunning−Hay basis set,71 which gives a slightly better energy for the 3s state of the atom. We first define the diabatic reference states to model the valence characters of the ionic |Ψion[(Li+)(H−)]⟩ and the three lowest covalent configurations |Ψ2s[(Li•)(H•)]⟩, Li(2p), and Li(3s). First, the ionic configuration exhibits a prolonged driving character in LiH bonding interactions in both the ground and excited states, and it crosses with all other covalent states toward the dissociation asymptotes. It has strong mixings at different regions with different covalent configurations, 6043

DOI: 10.1021/acs.jpclett.8b02472 J. Phys. Chem. Lett. 2018, 9, 6038−6046

Letter

The Journal of Physical Chemistry Letters

such as molecular dipoles and quadrupoles or the computation of nonadiabatic coupling terms in the diabatization process. It is also important to note that, unlike some ATD transformations, there is no issue with phase consistency when more than two states are involved because they are explicitly treated in the MSDFT Hamiltonian along with the required overlap integrals. GDAC is presented with the intention of providing a straightforward procedure for representing the relevant adiabatic PESs by the same number of diabatic states with valence characters corresponding to the dissociation products of the photochemical reaction. The GDAC diabatization method may be useful for defining chargelocalized states in studies of electron transfer and protoncoupled electron transfer reactions in proteins.

responsible for the interactions and avoided curve crossings of the four lowest states in consideration. To reflect its dominant role in the LiH bond in the ground state, the first diabatic reference state is chosen as |ΨD1 (R) = d11|Ψion[Li+)(H−)]⟩ + d21|Ψ2s[(Li•)(H•)]⟩, where the coefficients are optimized at a given geometry (R) with those for other states fixed at zero. Next, considering the ground-state dissociation products, we used a mixture of the covalent configuration |Ψ2s[(Li•)(H•)]⟩ and the reverse ionic configuration to define the second reference vector. The third reference vector was defined as an admixture of the 2p and the excited 2s of hydrogen configurations; |Ψ 3D (R) = d 3 3 |Ψ 2 p [(Li • )(H • )]⟩ + d73|Ψ2s/2s[(Li•)(H•)]⟩. Similarly, the coefficients d33 and d73were determined by solving a 2 × 2 secular equation at a given interatomic distance. The final diabatic reference state is a single 3s state, |ΨD4 (R) = |Ψ3s[(Li•)(H•)]⟩. The result of the GDAC decomposition to transform the four lowest adiabatic states is given in Figure 4. The key ionic−



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.8b02472.



Detailed discussion and solution of the minimization problem of eq 5, summary of the valence bond state definition for MSDFT calculations of the adiabatic states, definition of reference diabatic states, and optimized structures in the bond dissociation of methoxychloromethane (PDF)

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (A.G.). *E-mail: [email protected] (J.G.).

Figure 4. Computed diabatic states (U1−U4) from the generalized singular value (GDAC), reproducing exactly the adiabatic ground state (V1) and the first three excited states (V2−V4) of LiH as a function of interatomic distance.

ORCID

Adam Grofe: 0000-0002-8531-4396 Jiali Gao: 0000-0003-0106-7154 Notes

The authors declare no competing financial interest.

covalent curve crossings to produce the avoided curves between S0 and S1 and S1 and S2 are correctly generated, extending as far as 10 Å. Of course, the four diabatic states (U1−U4) reproduce exactly the four lowest adiabatic energies. In Figure 4, the “ionic” and 2s states cross at about 2.85 Å, switching bonding character from 2s to 2p VB states for Li. The curve crossing between the first and second excited states occurs at about 5 Å. For comparison, the two curve crossings were found to take place at 2.9 and 3.7 Å from the four-fold way approach9,35 and 3.2 and 5.3 Å from a diabatization based on minimization of nonadiabatic coupling terms.16 The second and third excited states have two avoided crossings at 2.5 and 10 Å, respectively, due to 3s and 3p hybridization mixed with lower states and 3s−ionic configuration interactions (Figure S1). Both crossings due to switching different valence characters are correctly represented by the GDAC diabatic state (U3 and U4). Overall, the GDAC approach based on GSVD can provide a good representation of the diabatic states in comparison with other methods. This procedure, which projects the reference valence-state functions onto the adiabatic space, differs from traditional approaches through adiabatic to diabatic transformations based on characteristic orbitals or restrained by some desired properties. In GDAC, the valence characters of diabatic states are defined by state functions rather than specific orbitals, and it does not require additional restraints



ACKNOWLEDGMENTS The work was supported in part by grants from the National Natural Science Foundation of China (Grant Number 91541124) and the Ministry of Science and Technology (Grant Number 2017YFB0203400). A.G. is supported by a National Postdoctoral Fellowship of China.



REFERENCES

(1) Baer, M. Beyond Born−Oppenheimer: Electronic Nonadiabatic Coupling Terms and Conical Intersections; Wiley: New York, 2006. (2) Yarkony, D. R. Diabolical conical intersections. Rev. Mod. Phys. 1996, 68, 985−1013. (3) Domcke, W.; Yarkony, D. R. Role of Conical Intersections in Molecular Spectroscopy and Photoinduced Chemical Dynamics. Annu. Rev. Phys. Chem. 2012, 63, 325−352. (4) Nakamura, H.; Truhlar, D. G. The direct calculation of diabatic states based on configurational uniformity. J. Chem. Phys. 2001, 115, 10353−10372. (5) Japser, A. W.; Kendrick, C. A.; Mead, C. A.; Truhlar, D. G. In Modern Trends in Chemical Reaction Dynamics: Eperiment and Theory (Part I); Yang, X., Liu, K., Eds.; World Scientific: Singapore, 2004; pp 329−391. (6) Koppel, H. Diabatic Representation: Methods for the Constrution of Diabatic Electronic States. In Conical Intersections: Electronic Structure, Dynamics and Spectroscopy; Domcke, W., Yarkony,

6044

DOI: 10.1021/acs.jpclett.8b02472 J. Phys. Chem. Lett. 2018, 9, 6038−6046

Letter

The Journal of Physical Chemistry Letters D. R., Koppel, H., Eds.; World Scientific Publishing Co.: Hackensack, NJ, 2004; pp 175−204. (7) Song, L.; Gao, J. On the Construction of Diabatic and Adiabatic Potential Energy Surfaces Based on Ab Initio Valence Bond Theory. J. Phys. Chem. A 2008, 112, 12925−12935. (8) Zhu, X. L.; Yarkony, D. R. Toward eliminating the electronic structure bottleneck in nonadiabatic dynamics on the fly: An algorithm to fit nonlocal, quasidiabatic, coupled electronic state Hamiltonians based on ab initio electronic structure data. J. Chem. Phys. 2010, 132, 104101. (9) Li, S. H. L.; Truhlar, D. G.; Schmidt, M. W.; Gordon, M. S. Model space diabatization for quantum photochemistry. J. Chem. Phys. 2015, 142, 064106. (10) Halvick, P.; Truhlar, D. G. A new diabatic representation of the coupled potential energy surfaces for sodium(3p 2P) + molecular hydrogen.fwdarw. sodium(3s 2S) + molecular hydrogen or sodium hydride + atomic hydrogen. J. Chem. Phys. 1992, 96, 2895−2909. (11) Jasper, A. W.; Zhu, C. Y.; Nangia, S.; Truhlar, D. G. Introductory lecture: Nonadiabatic effects in chemical dynamics. Faraday Discuss. 2004, 127, 1−22. (12) Cederbaum, L. S.; Koppel, H.; Domcke, W. Multimode Vibronic Coupling Effects in Molecules. Int. J. Quantum Chem. 1981, 20, 251−267. (13) Pacher, T.; Cederbaum, L. S.; Koppel, H. Approximately diabatic states from block diagonalization of the electronic Hamiltonian. J. Chem. Phys. 1988, 89, 7367. (14) Domcke, W.; Woywod, C. Direct Construction of Diabatic States in the Casscf Approach - Application to the Conical Intersection of the (1)a(2) and B-1(1) Excited-States of Ozone. Chem. Phys. Lett. 1993, 216, 362−368. (15) Domcke, W.; Woywod, C.; Stengle, M. Diabatic Casscf Orbitals and Wave-Functions. Chem. Phys. Lett. 1994, 226, 257−262. (16) Boutalib, A.; Gadea, F. X. Abinitio Adiabatic and Diabatic Potential-Energy Curves of the Lih Molecule. J. Chem. Phys. 1992, 97, 1144−1156. (17) Berriche, H.; Gadea, F. X. Ab-Initio Adiabatic and Diabatic Permanent Dipoles for the Low-Lying States of the Lih Molecule - a Direct Illustration of the Ionic Character. Chem. Phys. Lett. 1995, 247, 85−88. (18) Zhang, P.; Morokuma, K.; Wodtke, A. M. High-level ab initio studies of unimolecular dissociation of the ground-state N-3 radical. J. Chem. Phys. 2005, 122, 014106. (19) Viel, A.; Eisfeld, W.; Evenhuis, C. R.; Manthe, U. Photoionization-induced dynamics of the ammonia cation studied by wave packet calculations using curvilinear coordinates. Chem. Phys. 2008, 347, 331−339. (20) Mota, V. C.; Varandas, A. J. C. HN2((2)A ’) electronic manifold. II. Ab initio based double-sheeted DMBE potential energy surface via a global diabatization angle. J. Phys. Chem. A 2008, 112, 3768−3786. (21) Subotnik, J. E.; Yeganeh, S.; Cave, R. J.; Ratner, M. A. Constructing diabatic states from adiabatic states: Extending generalized Mulliken-Hush to multiple charge centers with Boys localization. J. Chem. Phys. 2008, 129, 244101. (22) Zhang, A. J.; Zhang, P. Y.; Chu, T. S.; Han, K. L.; He, G. Z. Quantum dynamical study of the electronic nonadiabaticity in the D +DBr → Br(Br*) + D2 reaction on new diabatic potential energy surfaces. J. Chem. Phys. 2012, 137, 194305. (23) Godsi, O.; Evenhuis, C. R.; Collins, M. A. Interpolation of multidimensional diabatic potential energy matrices. J. Chem. Phys. 2006, 125, 104105. (24) Ou, Q.; Bellchambers, G. D.; Furche, F.; Subotnik, J. E. Firstorder derivative couplings between excited states from adiabatic TDDFT response theory. J. Chem. Phys. 2015, 142, 064114. (25) Zhu, X. L.; Malbon, C. L.; Yarkony, D. R. An improved quasidiabatic representation of the 1, 2, 3(1)A coupled adiabatic potential energy surfaces of phenol in the full 33 internal coordinates. J. Chem. Phys. 2016, 144, 124312.

(26) Zhu, X. L.; Yarkony, D. R. Constructing diabatic representations using adiabatic and approximate diabatic data Coping with diabolical singularities. J. Chem. Phys. 2016, 144, 044104. (27) Atchity, G. J.; Ruedenberg, K. Determination of diabatic states through enforcement of configurational uniformity. Theor. Chem. Acc. 1997, 97, 47−58. (28) Nakamura, H.; Truhlar, D. G. Direct diabatization of electronic states by the fourfold way. II. Dynamical correlation and rearrangement processes. J. Chem. Phys. 2002, 117, 5576−5593. (29) Nakamura, H.; Truhlar, D. G. Extension of the fourfold way for calculation of global diabatic potential energy surfaces of complex, multiarrangement, non-Born-Oppenheimer systems: Application to HNCO(S-0,S-1). J. Chem. Phys. 2003, 118, 6816−6829. (30) Subotnik, J. E.; Vura-Weis, J.; Sodt, A. J.; Ratner, M. A. Predicting Accurate Electronic Excitation Transfer Rates via Marcus Theory with Boys or Edmiston-Ruedenberg Localized Diabatization. J. Phys. Chem. A 2010, 114, 8665−8675. (31) Alguire, E.; Subotnik, J. E. Diabatic couplings for charge recombination via Boys localization and spin-flip configuration interaction singles. J. Chem. Phys. 2011, 135, 044114. (32) Jin, Z. X.; Subotnik, J. E. Localized diabatization applied to excitons in molecular crystals. J. Chem. Phys. 2017, 146, 244110. (33) Simah, D.; Hartke, B.; Werner, H. J. Photodissociation dynamics of H2S on new coupled ab initio potential energy surfaces. J. Chem. Phys. 1999, 111, 4523−4534. (34) Fulscher, M. P.; Serrano-Andres, L. Quasi diabatic CASSCF state functions. Mol. Phys. 2002, 100, 903−909. (35) Hoyer, C. E.; Parker, K.; Gagliardi, L.; Truhlar, D. G. The DQ and DQ Phi electronic structure diabatization methods: Validation for general applications. J. Chem. Phys. 2016, 144, 194101. (36) Hoyer, C. E.; Xu, X. F.; Ma, D. X.; Gagliardi, L.; Truhlar, D. G. Diabatization based on the dipole and quadrupole: The DQ method. J. Chem. Phys. 2014, 141, 114104. (37) Mead, C. A.; Truhlar, D. G. Conditions for the definition of a strictly diabatic electronic basis for molecular systems. J. Chem. Phys. 1982, 77, 6090−6098. (38) Grofe, A.; Qu, Z. X.; Truhlar, D. G.; Li, H.; Gao, J. L. DiabaticAt-Construction Method for Diabatic and Adiabatic Ground and Excited States Based on Multistate Density Functional Theory. J. Chem. Theory Comput. 2017, 13, 1176−1187. (39) Mo, Y.; Gao, J. An Ab Initio Molecular Orbital-Valence Bond (MOVB) Method for Simulating Chemical Reactions in Solution. J. Phys. Chem. A 2000, 104, 3012−3020. (40) Gao, J.; Grofe, A.; Ren, H.; Bao, P. Beyond Kohn−Sham Approximation: Hybrid Multistate Wave Function and Density Functional Theory. J. Phys. Chem. Lett. 2016, 7, 5143−5149. (41) Cembran, A.; Song, L.; Mo, Y.; Gao, J. Block-localized density functional theory (BLDFT), diabatic coupling, and its use in valence bond theory for representing reactive potential energy surfaces. J. Chem. Theory Comput. 2009, 5, 2702−2716. (42) Grofe, A.; Chen, X.; Liu, W. J.; Gao, J. L. Spin-Multiplet Components and Energy Splittings by Multistate Density Functional Theory. J. Phys. Chem. Lett. 2017, 8, 4838−4845. (43) Ren, H. S.; Provorse, M. R.; Bao, P.; Qu, Z. X.; Gao, J. L. Multistate Density Functional Theory for Effective Diabatic Electronic Coupling. J. Phys. Chem. Lett. 2016, 7, 2286−2293. (44) Alter, O.; Brown, P. O.; Botstein, D. Generalized singular value decomposition for comparative analysis of genome-scale expression data sets of two different organisms. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 3351−3356. (45) Bohen, S. P.; Troyanskaya, O. G.; Alter, O.; Warnke, R.; Botstein, D.; Brown, P. O.; Levy, R. Variation in gene expression patterns in follicular lymphoma and the response to rituximab. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 1926−1930. (46) Abdi, H.: Generalized Singular Value Decomposition (GSVD). In Encyclopedia of Measurement and Statistics; Salkind, N., Ed.; Sage: Thousand Oaks, CA, 2007. (47) Petersen, K. B.; Pedersen, M. S. The Matrix Cookbook. http:// matrixcookbook.com (Nov. 15, 2012). 6045

DOI: 10.1021/acs.jpclett.8b02472 J. Phys. Chem. Lett. 2018, 9, 6038−6046

Letter

The Journal of Physical Chemistry Letters

properties for the (CaLi)(+) ionic molecule. Mol. Phys. 2016, 114, 1568−1582. (70) Dunning, T. H., Jr. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (71) Dunning, T. H.; Hay, P. J. In Methods of Electronic Structure Theory; Schaefer, H. F., III, Ed.; Plenum Press, 1977; Vol. 3; pp 1−28.

(48) Venghaus, F.; Eisfeld, W. Block-diagonalization as a tool for the robust diabatization of high-dimensional potential energy surfaces. J. Chem. Phys. 2016, 144, 114110. (49) King, H. F.; Stanton, R. E.; Kim, H.; Wyatt, R. E.; Parr, R. G. Corresponding orbitals and the nonorthogonality problems in molecular quantum mechanics. J. Chem. Phys. 1967, 47, 1936. (50) Mo, Y.; Gao, J. Ab initio QM/MM simulations with a molecular orbital-valence bond (MOVB) method: application to an SN2 reaction in water. J. Comput. Chem. 2000, 21, 1458−1469. (51) Mo, Y. R.; Bao, P.; Gao, J. L. Energy decomposition analysis based on a block-localized wavefunction and multistate density functional theory. Phys. Chem. Chem. Phys. 2011, 13, 6760−6775. (52) Cohen, A. J.; Mori-Sanchez, P.; Yang, W. T. Challenges for Density Functional Theory. Chem. Rev. 2012, 112, 289−320. (53) Farazdel, A.; Dupuis, M. On the Determination of the Minimum on the Crossing Seam of 2 Potential-Energy Surfaces. J. Comput. Chem. 1991, 12, 276−282. (54) Smith, D. M. A.; Rosso, K. M.; Dupuis, M.; Valiev, M.; Straatsma, T. P. Electronic coupling between heme electron-transfer centers and its decay with distance depends strongly on relative orientation. J. Phys. Chem. B 2006, 110, 15582−15588. (55) Tully, J. C. Molecular dynamics with electronic transitions. J. Chem. Phys. 1990, 93, 1061−1071. (56) Zhu, C. Y. Restoring electronic coherence/decoherence for a trajectory-based nonadiabatic molecular dynamics. Sci. Rep. 2016, 6, 24198. (57) Horke, D. A.; Watts, H. M.; Smith, A. D.; Jager, E.; Springate, E.; Alexander, O.; Cacho, C.; Chapman, R. T.; Minns, R. S. Hydrogen Bonds in Excited State Proton Transfer. Phys. Rev. Lett. 2016, 117, 163002. (58) Farmanara, P.; Radloff, W.; Stert, V.; Ritze, H. H.; Hertel, I. V. Real time observation of hydrogen transfer: Femtosecond timeresolved photoelectron spectroscopy in the excited ammonia dimer. J. Chem. Phys. 1999, 111, 633−642. (59) Farmanara, R.; Ritze, H. H.; Stert, V.; Radloff, W.; Hertel, I. V. Ultrafast photodissociation dynamics and energetics of the electronically excited H atom transfer state of the ammonia dimer and trimer. J. Chem. Phys. 2002, 116, 1443−1456. (60) Park, J. K.; Iwata, S. Ab initio study of photochemical reactions of ammonia dimer systems. J. Phys. Chem. A 1997, 101, 3613−3618. (61) Nangia, S.; Truhlar, D. G. Direct calculation of coupled diabatic potential-energy surfaces for ammonia and mapping of a fourdimensional conical intersection seam. J. Chem. Phys. 2006, 124, 124309. (62) Smith, A. D.; Watts, H. M.; Jager, E.; Horke, D. A.; Springate, E.; Alexander, O.; Cacho, C.; Chapman, R. T.; Minns, R. S. Resonant multiphoton ionisation probe of the photodissociation dynamics of ammonia. Phys. Chem. Chem. Phys. 2016, 18, 28150−28156. (63) Stwalley, W. C.; Zemke, W. T. Spectroscopy and Structure of the Lithium Hydride Diatomic-Molecules and Ions. J. Phys. Chem. Ref. Data 1993, 22, 87−112. (64) Gadea, F. X.; Boutalib, A. Computation and Assignment of Radial Couplings Using Accurate Diabatic Data for the Lih Molecule. J. Phys. B: At., Mol. Opt. Phys. 1993, 26, 61−74. (65) Gadea, F. X.; Leininger, T. Accurate ab initio calculations for LiH and its ions, LiH+ and LiH-. Theor. Chem. Acc. 2006, 116, 566− 575. (66) Bande, A.; Nakashima, H.; Nakatsuji, H. LiH potential energy curves for ground and excited states with the free complement local Schrodinger equation method. Chem. Phys. Lett. 2010, 496, 347−350. (67) Lamoudi, N.; Bouledroua, M.; Alioua, K.; Allouche, A. R.; Aubert-Frecon, M. Theoretical investigation of the lithium 2p < - 2s photoabsorption spectra perturbed by atomic hydrogen. Phys. Rev. A: At., Mol., Opt. Phys. 2013, 87, 052713. (68) Gim, Y.; Lee, C. W. Studies of singlet Rydberg series of LiH derived from Li(nl) + H(1s), with n ≤ 6 and l ≤ 4. J. Chem. Phys. 2014, 141, 144313. (69) Habli, H.; Mejrissi, L.; Ghalla, H.; Yaghmour, S. J.; Oujia, B.; Gadea, F. X. Ab initio investigation of the electronic and vibrational 6046

DOI: 10.1021/acs.jpclett.8b02472 J. Phys. Chem. Lett. 2018, 9, 6038−6046