We have recently introduced a diabatization scheme, which simultaneously fits and diabatizes adiabatic *ab initio* electronic wave functions, Zhu and Yarkony J. Chem. Phys. **140**, 024112 (2014). The algorithm uses derivative couplings in the defining equations for the diabatic Hamiltonian, **H**^{d}, and fits all its matrix elements simultaneously to adiabatic state data. This procedure ultimately provides an accurate, quantifiably diabatic, representation of the adiabatic electronic structure data. However, optimizing the large number of nonlinear parameters in the basis functions and adjusting the number and kind of basis functions from which the fit is built, which provide the essential flexibility, has proved challenging. In this work, we introduce a procedure that combines adiabatic state and diabatic state data to efficiently optimize the nonlinear parameters and basis function expansion. Further, we consider using direct properties based diabatizations to initialize the fitting procedure. To address this issue, we introduce a systematic method for eliminating the debilitating (diabolical) singularities in the defining equations of properties based diabatizations. We exploit the observation that if approximate diabatic data are available, the commonly used approach of fitting each matrix element of **H**^{d} individually provides a starting point (seed) from which convergence of the full **H**^{d} construction algorithm is rapid. The optimization of nonlinear parameters and basis functions and the elimination of debilitating singularities are, respectively, illustrated using the 1,2,3,4^{1}A states of phenol and the 1,2^{1}A states of NH_{3}, states which are coupled by conical intersections.

## I. INTRODUCTION

Recent advances in full dimensional quantum mechanical descriptions of nonadiabatic processes,^{1} portend, at least for molecules comprised of 10-15 atoms, accuracy sufficient to make computational simulations equal partners with experiments. For the foreseeable future, these simulations will not employ electronic structure data provided on the fly or directly, but rather from representations fit to high quality *ab initio* electronic structure data. Given the discontinuous energy gradients and singular interstate couplings that plague the adiabatic representation near conical intersections, quasi-diabatic representations of the electronic structure data will be used. Here, the attribute quasi, which we will suppress below, is meant to indicate that for molecules larger than diatoms rigorous diabatic bases, that is, bases for which the residual derivative coupling, defined mathematically in Section II, vanishes globally, do not exist.^{2–5} Reflecting this result, definitions of diabaticity^{6} include (global) minimization of derivative couplings, smoothness of molecular properties, and configuration uniformity. Diabatizations using smoothness of molecular properties fall into two classes: (i) those that diagonalize properties matrices and (ii) those that are equivalent to systems of orthogonality equations, with the former being a special case of the latter. There has been a growing interest in methods for determining diabatic bases for adiabatic electronic wave functions,^{7–14} an area with a long history.^{3,4,15–20}

Previously, we have introduced a diabatization procedure, which simultaneously fits and diabatizes *ab initio* electronic structure data, producing a coupled diabatic state representation denoted as **H**^{d}.^{21} A unique feature of this approach is that it uses derivative coupling data to determine diabatic states, so that the residual derivative coupling can be determined. This procedure ultimately provides an accurate, quantifiably diabatic, representation of the adiabatic electronic structure data but optimizing the large number of nonlinear parameters in the basis functions and adjusting the number and kind of the basis functions from which the fit is built, which provide the essential flexibility, has proved challenging. This work introduces a procedure, which uses adiabatic state and diabatic state data to efficiently optimize nonlinear parameters and basis function expansions. Also considered is the feasibility of using direct, properties based, diabatizations to initiate the algorithm. Here the attribute direct indicates that the diabatization at point R uses only information at that point. Incorporation of properties based diabatizations necessitates consideration of a much more general issue, the limitations imposed by the recently described^{22} singularities in the defining equations of properties based diabatizations. Owing to their debilitating effects and unanticipated existence, these singularities were termed diabolical singularities in Ref. 22, which is referred to as ZY1 below. A procedure for systematically removing these diabolical singularities is introduced.

Our analyses exploit the observation that if diabatic data are available, the much simpler task of fitting each matrix element of **H**^{d} individually provides a starting point (seed) from which convergence of the full **H**^{d} construction algorithm is rapid. The restructured algorithm can be used to assess and ultimately improve alternative approaches to diabatization, which do not use derivative coupling information.

Section II explains the procedures for fitting matrix elements of **H**^{d} using adiabatic state and diabatic state data and describes the relation to the above described issues. Section III introduces the issue of singularities in properties based diabatizations and provides a method to overcome them. This section also provides numerical examples. Section IV summarizes and describes directions for future work.

## II. THEORETICAL APPROACH

### A. H^{d} defining equations using adiabatic or approximate diabatic data

The diabatic Hamiltonian **H**^{d}, the electronic Hamiltonian in the diabatic state basis, $ \Psi \alpha d ( q ; R ) $, has the form

which can be reindex as

where **V** are the linear parameters determined by the fitting procedure, the *g*^{(l)} (** R**) are the monomial functions of an overcomplete set of internal coordinates,

*P*

^{u(l),v(l)}is the projection operator for the appropriate irreducible representation of the complete nuclear permutation inversion (CNPI) group,

^{23,24}and

**B**

^{u,v}are the

*N*×

^{state}*N*basis matrices with a 1 in the (

^{state}*u*,

*v*) and (

*v*,

*u*) matrix elements and the rest 0. Accompanying

**H**

^{d}is the electronic Schrödinger equation,

#### 1. Adiabatic state data

From Eq. (2), adiabatic energies, energy gradients, and derivative couplings ($ f j I , J , a , ( m ) $) are given, respectively, by

The diabatic states are given in terms of the adiabatic states, $ \Psi I a ( q ; R ) $, by

where **q** (**R**) denote the electronic (nuclear) coordinates. The residual derivative coupling is

where (*ab*) implies the quantity in question is determined from the *ab initio* wave functions. The residual coupling is related to the diabatic state coupling $ \Psi \alpha d ( q ; R ) | \u2207 k \Psi \beta d ( q ; R ) q $ by the following system of linear equations:

where

While $\delta f k I , J , a ( R ) $ is a natural measure of the error in the derivative coupling, near a conical intersection where a small change in geometry can produce a big change in $ f k I , J , a , ( x ) ( R ) $, it does not reflect the true accuracy of the representation.

The **V** are determined by the requirement,

which is satisfied in a least square sense for *j* = 0, 1, …, 3*N*^{at} − 6; *n* = 1, …, *N ^{point}*; 1 ≤

*K*≤

*L*≤

*N*. Here,

^{state}*j*= 0, which only applies for

*K*=

*L*, implies do nothing. The diabaticity of

**H**

^{d}arises from Eq. (5) for

*K*≠

*L*. The determination of

**V**by the iterative solution of Eq. (5) has been discussed previously.

^{25}Since the (

*I*,

*J*) in Eq. (3) represent adiabatic data, Eq. (5) treats all matrix elements of

**H**

^{d}simultaneously.

#### 2. Diabatic state data

When approximate diabatic data $ H \u0303 \alpha , \beta d ( R ( n ) ) $, $ \u2207 j H \u0303 \alpha , \beta d ( R ( n ) ) $, *j* = 1, …, 3*N*^{at} − 6, *n* = 1, …, *N ^{point}* are available, we can fit a single element of

**H**

^{d}denoted $ H \alpha , \beta D $ directly in a least squares sense. Undoing the transformation to adiabatic states in Eqs. (3a) and (3b) leads us to require

for each (α, β) where, *j* = 0, …, 3*N*^{at} − 6,

and *j* = 0 means do nothing. For the least squares solution to Eq. (6), the expansion coefficients $ V l \alpha , \beta $, see Eq. (1b), are given by the normal equations,^{26}

This procedure divides the least squares problem in Eq. (5) into much smaller blocks so that the procedure is much faster than 1 iteration of the full fitting procedure. Further, the absence of the **d**^{J} in Eqs. (6)–(8) makes them much simpler to solve than Eq. (5). Here and in Section II B, we will use **H**^{d} [**H**^{D}] to indicate a quasi diabatic Hamiltonian determined using Eq. (5) [Eq. (8)]. Otherwise, **H**^{d} and **H**^{D} are identical in form. We will use **H**^{d} to denote a generic diabatic Hamiltonian when no confusion will result. Seeded with the approximate **H**^{D} constructed by the solution to Eq. (8), convergence of the full algorithm to produce an **H**^{d} based on Eq. (5) is rapid.

It should be emphasized that the matrix element at a time (**H**^{D}) diabatization is not new, in fact its use in fitting diabatic representations is well established.^{27,28} However, since the **H**^{D} type fit treats each matrix element approximately, that is, in a least squares sense, a global **H**^{d} type optimization always lowers the RMS energy error compared to an exclusively **H**^{D} based procedure. Since the RMS adiabatic energy error is a defining quantity in a diabatic representation, the **H**^{d} method, which includes diabaticity enforcing equations, is preferred albeit much more challenging to implement.

### B. Overview of adiabatic data based determination of H^{d}

In order to explain the procedures to be introduced in this work, a brief overview of the method used to determine **H**^{d} is provided.

#### 1. Initialization

The initial set of **R**^{(n)}, denoted as the skeletal data, is a limited set of points including minima, saddle points, minimum energy conical intersections, and points in the product channels sufficient to coarsely describe the processes under consideration. In addition, paths between these well spaced points are required to help in guiding the diabatization of the electronic wave functions at these points. Sufficient data are required so that the emerging **H**^{d} gets the state ordering qualitatively correct.

#### 2. Constructing the domain of definition. Fill. A quasi classical surface hopping trajectories based procedure

As the name indicates the skeletal data are far from sufficient to describe nonadiabatic dynamics. To generate a sufficient set of points, the domain of definition, a quasi classical surface hopping trajectories guided procedure, termed Fill, which is similar to the Grow procedure introduced by Collins,^{29,30} is used. This procedure accomplishes two things. It constructs **H**^{d} by *gradually adding* **R**^{(n)} to the fit data set. If too many points are added in one step early in the process, the procedure for solving Eq. (5) does not converge. Fill also determines a set of **R**^{(n)} sufficient to simulate nuclear dynamics.

### C. Sources of diabatic data

Key to this work are the sources of the diabatic data. Two sources of approximate diabatic data are considered, one suitable for improving **H**^{d} and one that can facilitate construction of the skeletal data.

#### 1. Improving the *g*^{(l)}

Assume that an **H**^{d} has been constructed using an initial set of *g*^{(l)}, but the fit is inadequate, so that an extended and/or optimized set of *g*^{(l)} must be found. If the full set of data, the **R**^{(n)} comprising the domain of definition, is used to construct **H**^{d} expanded in a new basis, the procedure based on Eq. (5) is not in general convergent. Instead a Fill approach, analogous to that used in the initial construction, had to be used. This is quite time consuming.

This Fill procedure can now be avoided. The previous **H**^{d} can be used to provide the diabatic data for all **R**^{(n)} in the domain of definition. Using these data, the intermediate and as noted above much easier step of constructing **H**^{D} using Eq. (8) is performed. Then, that **H**^{D} is used to seed the full **H**^{d} determination (Eq. (5)). Prior to each fit, the linear dependent *g*^{(l)} are removed using the rank reducing procedure described in the Appendix. As shown in Sec. III, with **H**^{D} as the seed, convergence of the full **H**^{d} is rapid.

#### 2. Facilitating the initial start

As suggested above, it can be challenging to properly diabatize the skeletal data set without considerable user intervention. This bottleneck can be circumvented by using electronic structure data that have been approximately, and directly diabatized by, for example, a dipole properties based diabatization, including the Dipole-Quadrupole (DQ)^{13} method or the Boys localization (BL)^{6} method. In this case, the diabatic data $ H \u0303 d $ and $\u2207 H \u0303 d $ are given at each **R**^{(n)} in the skeletal data set by

and

where **U** is the adiabatic-to-diabatic (AtD) states transformation obtained from the properties diabatization. The transformation in Eq. (9b) defining $\u2207 H \u0303 d $ works well except near a conical intersection.

### D. Analyzing and improving H^{D}

The pointwise data used to facilitate the initial start are provided by an alternative diabatization procedure. Consequently, it is possible to use that data to construct an **H**^{D} that reproduces the diabatic energies of the alternative diabatization. That **H**^{D} can be used to compute the corresponding residual derivative coupling, via Eq. (4). It can also be used to seed the solution to Eq. (5) showing the extent to which a diabatization, which does not reflect the derivative couplings can be improved by the solution to Eq. (5). This issue will addressed computationally in a future work.

## III. APPLICATIONS

In this section, two examples are considered which illustrate the procedures and issues introduced here. In the first example, the basis optimization procedure is used to significantly improve the quality of the expansion, the *g*^{(l)}, in a large system, phenol, in its full 33 internal coordinates where it is challenging to design an optimum expansion using the Fill method alone.^{31} For the second example, the BL diabatization,^{6} a direct dipole based diabatization, will be applied to the 1,2^{1}A states of ammonia, which are coupled by a seam of conical intersections and for which an **H**^{d} has been previously reported.^{32} The ultimate goal of this treatment is to use BL or other direct dipole based diabatization data as an approximate **H**^{D} to seed the determination of **H**^{d}. The algorithm will also enable us to quantify the changes in properties based diabatizations attributable to explicit consideration of the derivative coupling. This will be the subject of a future work. Here, we investigate the somewhat surprising and intriguing question, whether BL determined direct dipole diabatization data are suitable for the task. This question arises because of the existence of singularities in the equations defining properties based diabatizations, denoted diabolical singularities. In this work, we demonstrate the debilitating effects of these singularities and show how they can be circumvented. Other diabatization procedures including the standard dipole moment properties diabatization,^{13,20,33} known as the Werner-Meyer approach^{34} and the DQ^{13} approach, will be considered in future work. All electronic structure calculations reported in this work employ the COLUMBUS electronic structure codes.^{35–37}

### A. Improving a representation

The photodissociation of phenol,

has been the subject of much recent research.^{38–46} In the example used here, it is described using a multireference configuration interaction single and double excitation (MRCISD) expansion of approximately 88 × 10^{6} configuration state functions (CSFs), in the full 33 internal coordinates using four quasi-diabatic states. The initial representation considered here, obtained from the Fill procedure discussed in Sec. II B 2, was based on 5425 **R**^{(n)}, 14 940 projected *g*^{(l)}, resulting in 602 131 equations, which were solved in a least squares sense. The initially chosen *g*^{(l)} representation had a 413.35 cm^{−1} RMS energy fitting error and a 9.93% error in the energy gradients. Neglecting the time to compute the electronic structure data, we estimate the wall clock time to obtain this representation using the grow approach is 72 h.

Table I demonstrates how starting from that initial representation, the representation was significantly improved by repeatedly modifying the basis, determining **H**^{D} and then using **H**^{D} as the seed for the determination of **H**^{d} using Eq. (5). In 14 such passes, the RMS energy fitting error was reduced to 263 cm^{−1} and the error in the energy gradients was reduced to 7.16%, while increasing the representation size to 28 766 *g*^{(l)}. Each of the 14 lines in Table I represents a reconfiguration of the *g*^{(l)}, an extension (passes 1-4, 8-13), or modification (passes 5-7, 14) of the preceding expansion. Note modifications (passes 5-7, 14) are, in principle, particularly challenging since some functions are removed and replaced. However, the procedure performs flawlessly in this instance. The degree of improvement going from **H**^{D} to **H**^{d}, that is, the value of the **H**^{d} optimization, as measured by the RMS energy error, can be anywhere from 5% to 50% or even larger. It took less than 8 h to converge each of the full optimizations including the **H**^{D} seed, an enormous saving, so that the whole optimization in Table I took less than 5 days, considerably less than the ∼40 days the Fill approach alone would require. Note that each representation is built on the preceding expansion, so that the different optimizations could not be run in parallel.

Pass . | Grad error (%) . | RMS energy error (cm^{−1})
. | N^{c}
. |
---|---|---|---|

0 | 9.93 | 413.35 | 14 940 |

1 | 9.11 | 398.4 | 16 108 |

2 | 8.92 | 393.07 | 16 740 |

3 | 8.61 | 368.42 | 19 224 |

4 | 8.50 | 358.7 | 20 532 |

5 | 8.11 | 340.8 | 20 532 |

6 | 8.00 | 331.69 | 20 532 |

7 | 7.95 | 327.94 | 20 532 |

8 | 7.79 | 318.56 | 21 712 |

9 | 7.65 | 305.27 | 23 100 |

10 | 7.54 | 292.25 | 24 532 |

11 | 7.45 | 287.63 | 25 828 |

12 | 7.29 | 275.53 | 27 340 |

13 | 7.19 | 266.83 | 28 716 |

14 | 7.16 | 263.68 | 28 716 |

Pass . | Grad error (%) . | RMS energy error (cm^{−1})
. | N^{c}
. |
---|---|---|---|

0 | 9.93 | 413.35 | 14 940 |

1 | 9.11 | 398.4 | 16 108 |

2 | 8.92 | 393.07 | 16 740 |

3 | 8.61 | 368.42 | 19 224 |

4 | 8.50 | 358.7 | 20 532 |

5 | 8.11 | 340.8 | 20 532 |

6 | 8.00 | 331.69 | 20 532 |

7 | 7.95 | 327.94 | 20 532 |

8 | 7.79 | 318.56 | 21 712 |

9 | 7.65 | 305.27 | 23 100 |

10 | 7.54 | 292.25 | 24 532 |

11 | 7.45 | 287.63 | 25 828 |

12 | 7.29 | 275.53 | 27 340 |

13 | 7.19 | 266.83 | 28 716 |

14 | 7.16 | 263.68 | 28 716 |

### B. Removing diabolical singularities

Diabolical singularities will be considered using a generalized form of the orthogonality equations introduced in Ref. 6 to construct BL diabatizations. As discussed in ZY1, this approach has the advantage of being able to describe the commonly used generalized Mulliken-Hush,^{47} DQ,^{13} and Werner-Meyer^{34} diabatization methods. The BL approach achieves diabatization by maximizing the norm of the difference of the electric dipole moment for pairs of states, extending the Boys/Ruedenberg^{48} orbital localization method to electronic states. This extremum condition is equivalent to a system of orthogonalization equations. See Eq. (11), Ref. 6, and ZY1 for details.

To construct a generalized BL diabatic representation define the matrix elements of a set of Hermitian operators *O*^{(n)} by

and also define

where *d* (*a*) denote diabatic (adiabatic) states. Then, a generalized BL diabatic representation can be defined as maximizing

It was shown in Ref. 6 that this is equivalent to the orthogonality conditions,

In the Boys localization formulation of Subotnik *et al.*,^{6} we take *m* = 3, and **O** = μ, the electronic dipole operator, that is, O^{(1)} = μ^{(x)}, O^{(2)} = μ^{(y)}, and O^{(3)} = μ^{(z)}. For *N ^{state}* = 2, the solution to Eq. (11b) is

where *θ*_{i,j} is the rotation angle for the two state AtD transformation. The derivative coupling, assumed to be in the adiabatic representation, is given by

Note that *θ*_{i,j} is singular on the 3*N ^{at}* − 8 =

*N*− 2 dimensional space for which the numerator and denominator in Eq. (11c) vanish. The existence of this seam of singularities is a consequence of the defining equations, Eqs. (11a)–(11c). The precise locus of the seam is determined by the choice of operators

^{int}**O**and the choice of wave functions.

Low energy points on this subspace can be found by minimizing the following Lagrangian:

The minimization of *L*(**R**, α, β) requires the gradients of $ O i , j n , ( a ) ( R ) $, which are costly to evaluate even with modern analytic gradient techniques.

As we will demonstrate below (see also ZY1), the singularities in Eq. (11c) are in many respects very similar to those encountered in conical intersections of potential energy surfaces (PESs). Further, they cause the resulting diabatic surfaces to be significantly deteriorated, which is not useful as diabatic or seed data, even in regions away from the actual points of singularity. This we show is a consequence of a topological or geometric phase effect^{49} whose precise form depends on the diabatization procedure. Fortunately, unlike conical intersections of PESs, there is no practical need to determine the precise location of these singularities. These singularities, as well as the pathological artifacts they create on the diabatic surfaces, can be removed through a surprisingly simple procedure without the need to calculate any additional *ab initio* quantities.

The 1,2^{1}A states of ammonia, which will be considered, describe an archetypical example of nonadiabatic photodissociation,

The calculations reported here build on our previous work on this system.^{22,50,51} They were performed using molecular orbitals obtained from a state averaged multiconfiguration self consistent field treatment that averaged two states with equal weights and used wave functions obtained from a 8 electron in 9 orbital, complete active space expansion. State averaging and diabatic states have been discussed in Ref. 52. The atomic orbital basis consists of Dunning’s aug-cc-pVTZ basis plus an extra Rydberg s-function on nitrogen and Dunning’s aug-cc-pVTZ basis on the hydrogens. Dynamic correlation was incorporated at the second order configuration interaction levels (single and double excitations relative to the active space) with the interacting space restrictions. The resulting expansion, used for all the electronic structure data, consists of 33.9 × 10^{6} CSFs.^{50}

The potential problems with direct properties based diabatizations owing to diabolical singularities are described below. We use the BL diabatization with *N*^{state} = 2, and as noted above start with *m* = 3 and *O*^{(1)}, *O*^{(2)}, *O*^{(3)} = *μ ^{x}*,

*μ*,

^{y}*μ*. In Fig. 1(a), we report the norms of adiabatic state derivative couplings ||

^{z}**f**

^{i,j,a}|| determined from the

*ab initio*wave functions and the BL procedure, Eq. (11d), along a linear synchronous transit path connecting the planar D

_{3h}2

^{1}A equilibrium structure (R(N − H) = 1.9814 a

_{0}, E

_{1}= 0, E

_{2}= 43 884.72 cm

^{−1}) and the C

_{2v}minimum energy conical intersection (R(N − H

^{1}) = 3.7206 a

_{0}, R(N − H

^{2}) = R(N − H

^{3}) = 1.9360 a

_{0}, ∠H

^{1}NH

^{2}= ∠H

^{1}NH

^{3}= 54.60 ˚, E

_{1}= E

_{2}= 39 075.99 cm

^{−1}). The structures and energies are obtained from Ref. 50. The path was chosen because of its importance in the quenching dynamics and because in ZY1 there was shown to be a diabolical singularity on the path. The diabolical singularity is of the accidental symmetry allowed type since the numerator in Eq. (11c) vanishes by symmetry for planar geometries. Thus, only one condition, the vanishing of the denominator, is required and a 1-dimensional path is sufficient to find it. As explained in ZY1, the data in Fig. 1(a) indicate the complete breakdown of the BL determined diabatic representation in the vicinity of the diabolical singularity at R(N − H) = 2.3 a

_{0}.

^{22}The solid grey line in Fig. 1(a), the BL determined ||

**f**

^{i,j,a}||, exhibits a peak near R(N − H) = 3.72 a

_{0}, which precisely and correctly reproduces the singular

*ab initio*determined ||

**f**

^{i,j,a,(ab)}|| at the true conical intersection. However, the peak at R(N − H) ∼ 2.3 a

_{0}is entirely fallacious and, as shown in ZY1, is due to the singularity in Eq. (11c). As indicated in Fig. 1(a), the true (

*ab initio*) adiabatic derivative coupling is quite small in that region. The BL determined derivative coupling is seen to correspond to diabatic state derivative coupling.

Fig. 2 considers the consequences of the diabolical singularity in the present context reporting an analysis of the BL diabatization in a 2 dimensional subspace of nuclear coordinate space, illustrating an issue, not previously described, but essential to this work and attributable to diabolical singularities, a geometric phase effect. The two dimensional space used here has its origin at the C_{2v} minimum energy conical intersection and is spanned by ϕ the out of plane angle defined in the Fig. 2 caption, and R(N-H^{1}) which are, respectively, the **h** vector and approximately the **g** vector, comprising the branching plane^{53} of the C_{2v} minimum energy crossing point. It should be pointed out that this plane does not contain the diabolical singular point in Fig. 1(a) but does contain a nearby point on the seam.

Figs. 2(a) and 2(b) report the norm of the BL determined $ G 1 , 2 ( x ) $ and $ O 1 , 2 ( x ) $, before the diabatization, *x* = *a*, and after the diabatization, *x* = *d*, respectively. Figs. 2(c) and 2(d), respectively, report the diabatic PESs and the coupling surface. It is evident from Fig. 2(a) that before the BL procedure, the $|| G 1 , 2 ( x ) ||$ and $|| O 1 , 2 ( x ) ||$, *x* = *a* in the region near the minimum energy conical intersection, *R*(N − H) ∼ 3.75 a.u. are discontinuous, as a result of the strong coupling attributable to the true conical intersection. This should be contrasted with this region in the diabatic representation, Fig. 2(b), where $|| G 1 , 2 ( x ) ||$ and $|| O 1 , 2 ( x ) ||$, *x* = *d* are well behaved. On the other hand, in the shorter bond length region, near *R*(N − H) ∼ 2.3 a.u., where an accidental crossing of $|| G 1 , 2 ( x ) ||$ and $|| O 1 , 2 ( x ) ||$, *x* = *a*, and hence a diabolical singularity occur, the BL diabatization, as a result of the singularity, produces a double cone shaped intersection of the dipole norms that is very similar to a conical intersection of PESs. See Fig. 2(b). The resulting diabatic representation, Figs. 2(c) and 2(d), is smooth near the real (energetic, *R*(N − H) ∼ 3.75 a.u.) conical intersection but discontinuous near the diabolical singularity (*R*(N − H) = 2.3 a.u.). In fact, in the entire short bond length region, the diabatic representation is strongly coupled, causing extra features in the diabatic potential that are very difficult to fit, even in regions significantly away from the actual singularity. The most significant problem that prevents the diabatic representation from being truly diabatic and conveniently fit to an analytical representation is a geometric phase effect that causes the two diabatic potentials to be topographically connected. This can be seen in Fig. 2(c). Due to the singularity, much like the case of geometric phase in conical intersections, the AtD transformation is no longer single valued, but accumulates a relative phase of π/2 along a closed loop circling the singularity. As a result, which can be seen in Fig. 2(c), the two diabatic representations switch along such a loop. On the short R(N-H) side of the singularity, the yellow diabatic surface continuously joins into the blue surface. It is therefore not possible to uniquely identify each diabatic state. The consequence of this type of geometric phase effect is that even if one ignores region where the singularity occurs, the global topological property of the diabatic representation still prevents a satisfactory fit to any single valued analytical form.

Thus, the diabolical singularity destroys the utility of the associated diabatization. In order to fix this problem, we introduce a simple modification that can be applied to any operator based direct diabatization method, that can be rewritten in the form of Eq. (11b), that being adding as an extra component *wH ^{e}* to the operator vector in Eqs. (10a) and (10b). Here,

*w*is a scale factor, and

*H*is the electronic Hamiltonian. For the BL diabatization, this amounts to adding a fourth operator

^{e}*O*

^{(4)}=

*wH*. In this case, Eq. (11b) becomes

^{e} where we have made explicit the contribution of the *O*^{(4)} term in Eq. (11b). It should be noted that the extra term, *O*^{(4)} is not added as an ad hoc means of fixing the singularities. Rather, its utility can be derived from our previous analysis of the problem in ZY1. From Eq. (13), this modified BL diabatization is seen to add a positive term (coming from $ G i , j ( a ) $) to the denominator. The addition of *wH ^{e}* as an extra component in any direct operator diabatization method eliminates, for

*w*sufficiently large, the singularities but at the same time preserves the correct behavior of these direct diabatization methods at conical intersections of the PESs. We explain this below.

In the two state case, this problem manifests itself in regions where two vectors $ G 1 , 2 ( x ) $ and $ O 1 , 2 ( x ) $ *x* = *a*, have the same norm. A fix to this problem should therefore prevent the norm of $ G 1 , 2 ( x ) $ and $ O 1 , 2 ( x ) $ *x* = *a*, from being the same, while at the same time maintaining the correct behavior near the true conical intersections. It is evident that *wH ^{e}* precisely fulfills these requirements because: (i)

*H*does not contribute to either $ G 1 , 2 ( x ) $ or $ O 1 , 2 ( x ) $

^{e}*x*=

*a*near points on the conical intersection seams because of the degeneracy, and (ii) away from the conical intersection seam,

*H*contributes only to $ G 1 , 2 ( x ) $ vector, but not $ O 1 , 2 ( x ) $ vector,

^{e}*x*=

*a*, thus, creating a difference between the norms of the $ O 1 , 2 ( a ) $ and $ G 1 , 2 ( a ) $ vectors away from the conical intersection seam.

Unless the seam of singularities in Eq. (11b) intersects the conical intersection seam itself, which would be crossing between two *N*^{int} − 2 dimensional subspaces and thus quite rare, for any individual point a sufficient scaling factor for the *wH ^{e}* operator can ensure that the magnitude of $ G i , j ( x ) $ is always larger than $ O i , j ( x ) $ in the adiabatic representation (

*x*=

*a*). It is readily seen that if $ G i , j ( x ) $ and $ O i , j ( x ) $ have the same norm in the diabatic representation, then, the norms will be the same in any representation at this particular geometry. Therefore, when $ G i , j ( x ) $ and $ O i , j ( x ) $ norms are different in adiabatic representation, they cannot be the same in diabatic representation. Thus, using the largest value of scaling factor for all data points considered, one can ensure that the $ G i , j ( x ) $ and $ O i , j ( x ) $,

*x*=

*a*, surfaces never cross, thus, avoiding the singularity.

It is worth pointing out that for most operators such as dipole and angular momentum, neither $ G i , j ( x ) $ nor $ O i , j ( x ) $ surfaces will possess very large norms for any dynamically relevant geometries. Therefore, the scaling factors are in fact quite small, and the additional term does not strongly perturb the representation in regions that behave correctly.

Figs. 1(a)-1(d) compare the *ab initio* and BL determined norm of the derivative couplings ||**f**^{i,j,a}|| for *w* = 0, 4, 8, 10, respectively. This comparison quantifies the accuracy of the modified BL diabatization and the effect of the increasing weight *w* on the true and spurious singularities in ||**f**^{i,j,a}||. The diabolical singularity for R(N − H) ∼ 2.4 a.u. is evident for *w* = 0, 4. See Figs. 1(a) and 1(b). However for *w* ≥ 8, the $ G 1 , 2 ( a ) $ surface is always above the $ O 1 , 2 ( a ) $ surface and the resulting representation becomes free of the spurious spike in the BL determined derivative couplings near R(N − H) = 2.4 a.u. See Figs. 1(c) and 1(d). Figs. 1(a) and 1(b) show that despite the spurious singularity for short R(N-H), near the conical intersection, R(N − H) ∼ 3.75 a.u., the *ab initio* and BL determined ||**f**^{I,J,a}|| are in good agreement as expected. However, away from the conical intersection the BL determined ||**f**^{I,J,a}|| (solid lines) coupling is in poor agreement with the exact, small, *ab initio* determined derivative couplings (open circles). For *w* ≥ 8, the modified BL diabatization only degrades slightly for short R(N-H) and is certainly useful as a seed. The spurious spike in the BL determined **f**^{I,J,a} for small *w* is anathema for nonadiabatic dynamics. For a discussion of **f**^{d} in Fig. 1, see ZY1.

Returning to the AtD transformation, comparing Figs. 3(a)-3(d) with Figs. 2(a)-2(d) confirms that *w* = 4 is not sufficiently large to remove the singularity near R(N − H) = 2.4 a.u. However, that is not the case in Figs. 4(a)-4(d) and 5(a)-5(d) where *w* = 8 and *w* = 10, respectively. There it is evident that the modified operator no longer possess a spurious conical intersection like feature between the $ G i , j ( x ) $ and $ O 1 , 2 ( x ) $ *x* = *a* surfaces. The corresponding diabatic $ G i , j ( x ) $ and $ O 1 , 2 ( x ) $ *x* = *d*, are smooth and well behaved in Figs. 4(b) and 5(b). As result the diagonal and off-diagonal diabatic potentials Figs. 4(c), 4(d), 5(c), and 5(d) are smooth and well behaved in marked contrast to the results in Figs. 3(c) and 3(d) where *w* = 4. Thus, modified BL operators with *w* > 7 create a good approximate diabatic representation that can be used as an initial seed in the **H**^{d} fitting procedure.

The analysis in this section, which depends on Eq. (11), is restricted to *N*^{state} = 2 but otherwise is quite general. The extension of this analysis to the *N*^{state} > 2 case is a topic of current research.

## IV. SUMMARY AND CONCLUSIONS

A procedure is introduced to improve the performance and ease of use of an algorithm for the construction of a coupled quasi-diabatic state Hamiltonian, **H**^{d}, which represents adiabatic potential energy surfaces coupled by conical intersections. As a result of using derivative couplings in the defining equations, **H**^{d} is quantifiably quasi-diabatic. We found that the task of fitting each element of these representations individually to approximately diabatic data to produce a quasi diabatic representation **H**^{D} is much simpler than the (preferred) task of fitting all elements of **H**^{d} simultaneously to adiabatic state data. This is a consequence, at least in part, of the need to include the eigenvectors of the adiabatic states in the later procedure. Despite this, we found that **H**^{D} provides a starting point (seed) from which the iterative construction of **H**^{d} converges rapidly. This provides an efficient means of optimizing the large number of nonlinear parameters and adjusting the number and kind of basis functions, which provide the flexibility in the functional representation of **H**^{d}.

We also considered the effects of the recently discovered diabolical singularities on the diabatic representations produced by commonly used properties (principally dipole) based diabatizations. The effects of the diabolical singularities, which are a consequence of the equations defining the diabatizations, were shown to be quite deleterious compromising the utility of the constructed diabatizations. The deleterious effects were shown to be attributable to a type of geometric phase effect. A method to eliminate these singularities was introduced and tested.

In future work, the general utility of the **H**^{D} seed introduced here will be assessed. Of particular interest will be its use in assessing and improving diabatizations, which do not use derivative coupling information. The changes induced in those diabatizations by the explicit minimization of derivative couplings errors will be of particular interest. The determination of the prevalence and impact of diabolical singularities will provide a second direction of investigation.

## Acknowledgments

This work was supported by National Science Foundation Grant No. CHE-1361121 to D.R.Y.

### APPENDIX: RANK REDUCTION

For each block (*α*, *β*) of **H**^{d}, the null space satisfies for each (*n* = 1, …, *N ^{point}*,

*j*= 0, …,

*N*= 3

^{int}*N*− 6),

^{at}where

and ∇_{0} means do nothing. The pair (*n*, *j*) runs over the number of equations *N ^{eq}* and

*N*> >

^{eq}*N*. This procedure has the following desirable properties. It reflects the

^{c}**R**

^{n}used in the fit-initially removing large numbers of unneeded functions and then decreasing the number of dependent functions as

*N*grows; and it can be used to eliminate functions before the fit is carried out. Because each block is treated separately, the null-space analysis is a series of eigenvalue decompositions of matrices considerable smaller than

^{point}**H**

^{d}.

The disadvantage of this approach is that for each block of **H**^{d} the result is expressed in the eigenbasis, which complicates the elimination of individual basis functions. Consequently, we have implement the following rank-revealing QR decomposition:^{54}

where for clarity the dimensions are given parenthetically. Here, the columns of **Q** are orthonormal, P^{†} is a permutation matrix, and **R** (this is standard notation and it is unrelated to the R used to locate atoms) has the form

where the diagonal elements of R satisfy *ρ*_{l,l} > *ρ*_{l−1,l−1}, and ||**R**_{22}|| < tol. Thus, once a tolerance (tol) is picked, the last *N ^{c}* −

*k*columns are linearly dependent on the first

*k*columns. Then, discarding the last

*N*−

^{c}*k*columns of

*P*provides a reduced set of polynomials for the expansion in Eq. (1).