Stability conditions and local minima in multicomponent Hartree-Fock and density functional theory

Yang, Yang; Culpitt, Tanner; Tao, Zhen; Hammes-Schiffer, Sharon

doi:10.1063/1.5040353

Multicomponent quantum chemistry allows the quantum mechanical treatment of electrons and specified protons on the same level. Typically the goal is to identify a self-consistent-field (SCF) solution that is the global minimum associated with the molecular orbital coefficients of the underlying Hartree-Fock (HF) or density functional theory (DFT) calculation. To determine whether the solution is a minimum or a saddle point, herein we derive the stability conditions for multicomponent HF and DFT in the nuclear-electronic orbital (NEO) framework. The gradient is always zero for an SCF solution, whereas the Hessian must be positive semi-definite for the solution to be a minimum rather than a saddle point. The stability matrices for NEO-HF and NEO-DFT have the same matrix structures, which are identical to the working matrices of their corresponding linear response time-dependent theories (NEO-TDHF and NEO-TDDFT) but with a different metric. A negative eigenvalue of the stability matrix is a necessary but not sufficient condition for the corresponding NEO-TDHF or NEO-TDDFT working equation to have an imaginary eigenvalue solution. Electron-proton systems could potentially exhibit three types of instabilities: electronic, protonic, and electron-proton vibronic instabilities. The internal and external stabilities for theories with different constraints on the spin and spatial orbitals can be analyzed. This stability analysis is a useful tool for characterizing SCF solutions and is helpful when searching for lower-energy solutions. Initial applications to HCN, HNC, and 2-cyanomalonaldehyde, in conjunction with NEO ∆SCF calculations, highlight possible connections between stationary points in nuclear coordinate space for conventional electronic structure calculations and stationary points in orbital space for NEO calculations.

I. INTRODUCTION

The theoretical development of multicomponent quantum chemistry methods has garnered increased attention in recent years.^1–13 In contrast to the conventional electronic structure theory, which treats only electrons quantum mechanically, multicomponent quantum chemistry enables the quantum mechanical treatment of more than one type of particle, such as electrons and protons. In this approach, nuclear quantum effects such as delocalization¹⁰ and zero-point energy¹¹ associated with the quantum nuclei are included directly in a single quantum chemistry calculation. Moreover, the widely used Born-Oppenheimer approximation, which treats electrons and nuclei separately, can be avoided for these specified nuclei. Thus, nonadiabatic processes such as proton-coupled electron transfer (PCET)^14,15 could potentially be described with this approach.

Two branches of multicomponent quantum chemistry have been developed and actively investigated: wave function theory and density functional theory (DFT). The wave function theory branch¹ starts with multicomponent Hartree-Fock (HF) theory, which is based on the product of an electronic and a nuclear determinant, and can be systematically improved up to the exact multicomponent full configuration interaction (CI) method. The DFT branch² adopts the multicomponent one-particle densities as basic variables and is in principle exact, although the (exchange-)correlation functionals are unknown and need to be approximated.

The nuclear-electronic orbital (NEO) method⁵ is one of the promising multicomponent quantum chemistry approaches. It treats electrons and specified key protons quantum mechanically, while treating two or more nuclei classically to avoid difficulties with translations and rotations. Wave function-based NEO methods have been extensively studied in the past two decades, including NEO-HF,⁵ multiconfigurational self-consistent-field (NEO-MCSCF) theory,⁵ second-order Moller-Plesset perturbation (NEO-MP2) theory,¹⁶ NEO-CI,⁵ explicitly correlated HF (NEO-XCHF) theory,⁶ and reduced XCHF (NEO-RXCHF).⁹ Unfortunately, these methods are unable to accurately describe delocalized proton densities or relative energies in a computationally tractable manner. By contrast, Kohn-Sham NEO-DFT^{7,10,11,17,18} is more computationally practical, with the computational cost and scaling similar to those of conventional electronic DFT. Recent work has shown that delocalized proton densities and proton affinities can be described accurately with the epc17-1¹⁰ and epc17-2¹¹ electron-proton correlation functionals, respectively. Furthermore, recently the multicomponent time-dependent DFT method within the NEO framework, denoted NEO-TDDFT, has been developed.¹² Because both electrons and protons can be excited in this formalism, NEO-TDDFT could potentially describe electron-proton vibronic excitations, which are important in nonadiabatic excited state processes such as photoinduced PCET.^19,20

As in conventional electronic structure theory, the converged SCF procedure identifies a stationary point, but this stationary point may not correspond to the global minimum with respect to the molecular orbital coefficients. The stationary point is a minimum only if the corresponding Hessian is positive semi-definite; otherwise the stationary point corresponds to a saddle point that is considered to be an instability. Moreover, as in conventional electronic TDDFT, imaginary eigenvalues may occur in NEO-TDDFT calculations, indicating the instability of the underlying stationary SCF solution. The characterization of the nature of the stationary states found in NEO-HF and NEO-DFT calculations requires a stability analysis,^21–28 which involves the calculation of the Hessian within the multicomponent orbital space. This stability analysis also facilitates the search for a lower-energy solution and can be extended to higher levels of correlated theory.^29,30 Analogous issues related to SCF stability, convergence problems, and circumventing SCF solutions corresponding to saddle points or local minima have been explored in the context of conventional electronic structure theory.^23,24,31

In this paper, we investigate the stability conditions for multicomponent quantum chemistry in the NEO framework. The stability conditions for NEO-HF and NEO-DFT are derived in Sec. II. The relationship between negative eigenvalues of the stability matrix and imaginary eigenvalues of the corresponding linear response time-dependent theory is discussed in Sec. III. The instabilities with respect to electron orbital rotation, proton orbital rotation, and coupled electron-proton orbital rotation are discussed in Sec. IV. The internal and external stabilities for theories with different spatial and spin constraints are presented in Sec. V. Examples of stability analyses for the FHF⁻ molecule, the conversion between HCN and HNC, and hydrogen transfer in 2-cyanomalonaldehyde are presented in Sec. VI. These studies utilize NEO ∆SCF calculations³² to identify higher-energy stationary solutions corresponding to local minima. Concluding remarks are provided in Sec. VII.

II. BASIC THEORY

Both NEO-HF and Kohn-Sham NEO-DFT are based on a product of two determinants corresponding to electrons and protons, respectively,

\begin{matrix} | Ψ_{0} ⟩) = | Ψ_{0}^{e} ⟩) | Ψ_{0}^{p} ⟩) . \end{matrix}

(1)

The determinants are composed of N^e and N^p occupied electron orbitals $| i ⟩)$ and proton orbitals $| I ⟩)$ ⁠, respectively. Throughout this paper, we use lowercase (uppercase) letters to indicate electron (proton) orbitals and superscripts e (p) to indicate quantities related to electrons (protons). We also adopt the convention that i, j, k, l denote occupied orbitals, a, b, c, d denote virtual orbitals, and p, q, r, s denote general orbitals.

Following the procedure in Ref. 21, a small rotation of the occupied orbitals can be expressed as

| i^{'} ⟩) = | i ⟩) + \sum_{a} c_{i}^{a} | a ⟩) = (1 + \sum_{a} c_{i}^{a} a_{a}^{†} a_{i}) | i ⟩),

(2)

| I^{'} ⟩) = | I ⟩) + \sum_{A} d_{I}^{A} | A ⟩) = (1 + \sum_{A} d_{I}^{A} b_{A}^{†} b_{I}) | I ⟩),

(3)

where c and d with associated subscripts and superscripts denote the amplitudes of mixing for electron and proton orbitals, respectively, $a^{†}$ (a) denote creation (annihilation) operators for electrons, and $b^{†}$ (b) denote the corresponding operators for protons. Note that the primed orbitals expressed in Eqs. (2) and (3) are not normalized. Without loss of generality, only rotations mixing virtual orbitals and occupied orbitals are considered because the mixing between occupied orbitals can always be rotated back, and the energy is invariant to the rotations within the occupied orbital space.

The multicomponent ansatz under the orbital mixing can be approximated by keeping terms to second order with respect to the mixing amplitudes c and d,

\begin{matrix} | Ψ_{0}^{'} ⟩) & \approx & () 1 + \sum_{i a} c_{i}^{a} a_{a}^{†} a_{i} + \sum_{I A} d_{I}^{A} b_{A}^{†} b_{I} \\ + \frac{1}{2} \sum_{i j a b} c_{i}^{a} c_{j}^{b} a_{b}^{†} a_{j} a_{a}^{†} a_{i} \\ + \frac{1}{2} \sum_{I J A B} d_{I}^{A} d_{J}^{B} b_{B}^{†} b_{J} b_{A}^{†} b_{I} \\ + \sum_{i a I A} c_{i}^{a} d_{I}^{A} b_{A}^{†} b_{I} a_{a}^{†} a_{i})) | Ψ_{0} ⟩) . \end{matrix}

(4)

Note that this choice of orbital rotation is not unique and can be either unitary or non-unitary. Herein we follow the orbital mixing procedure in Ref. 21.

A. NEO-HF stability conditions

NEO-HF theory is based on the Hamiltonian

Ĥ = {\hat{T}}^{e} + {\hat{T}}^{p} + {\hat{V}}^{e e} + {\hat{V}}^{p p} + {\hat{V}}^{e p} + {\hat{v}}_{e x t}^{e} + {\hat{v}}_{e x t}^{p},

(5)

where the terms represent the kinetic energy of the electrons, the kinetic energy of the protons, the electron-electron Coulomb repulsion, the proton-proton Coulomb repulsion, the electron-proton Coulomb attraction, the interaction of the electrons with the external potential, and the interaction of the protons with the external potential, respectively. The energy change with respect to the orbital rotation is

Δ E = \frac{⟨ Ψ_{0}^{'} | Ĥ | Ψ_{0}^{'} ⟩}{⟨ Ψ_{0}^{'} | Ψ_{0}^{'} ⟩} - \frac{⟨ Ψ_{0} | Ĥ | Ψ_{0} ⟩}{⟨ Ψ_{0} | Ψ_{0} ⟩} .

(6)

Evaluating the above expression, the first-order energy change with respect to c and d is

\begin{matrix} Δ E^{(1)} & = & \sum_{i a} c_{i}^{a} ⟨ Ψ_{0} | (Ĥ - E) a_{a}^{†} a_{i} | Ψ_{0} ⟩ + c c \\ + \sum_{I A} d_{I}^{A} ⟨ Ψ_{0} | (Ĥ - E) b_{A}^{†} b_{I} | Ψ_{0} ⟩ + c c, \end{matrix}

(7)

where cc denotes the complex conjugate and E is the energy of the system before rotation [i.e., the second term of Eq. (6)]. When an SCF solution is obtained, this first-order energy change is zero due to Brillouin’s theorem for NEO-HF, which is briefly proven in the supplementary material. The second-order energy change is

\begin{matrix} Δ E^{(2)} & = & \frac{1}{2} \sum_{i j a b} c_{i}^{a} c_{i}^{b} ⟨Ψ_{0}| (\hat{H} - E) a_{b}^{†} a_{j} a_{a}^{†} a_{i} |Ψ_{0}⟩ + c c \\ + \frac{1}{2} \sum_{I J A B} d_{I}^{A} d_{J}^{B} ⟨Ψ_{0}| (\hat{H} - E) b_{B}^{†} b_{J} b_{A}^{†} b_{I} |Ψ_{0}⟩ + c c \\ + \sum_{i a I A} c_{i}^{a} d_{I}^{A} ⟨Ψ_{0}| (\hat{H} - E) b_{A}^{†} b_{I} a_{a}^{†} a_{i} |Ψ_{0}⟩ + c c \\ + \sum_{i j a b} c_{i}^{a} {(c_{j}^{b})}^{*} ⟨Ψ_{0}| a_{j}^{†} a_{b} (\hat{H} - E) a_{a}^{†} a_{i} |Ψ_{0}⟩ \\ + \sum_{i J A B} d_{I}^{A} {(d_{J}^{B})}^{*} ⟨Ψ_{0}| b_{J}^{†} b_{B} (\hat{H} - E) b_{A}^{†} b_{I} |Ψ_{0}⟩ \\ + \sum_{i a J B} c_{i}^{a} {(d_{J}^{B})}^{*} ⟨Ψ_{0}| b_{J}^{†} b_{B} (\hat{H} - E) a_{a}^{†} a_{i} |Ψ_{0}⟩ \\ + \sum_{b j I A} {(c_{j}^{b})}^{*} d_{I}^{A} ⟨Ψ_{0}| a_{j}^{†} a_{b} (\hat{H} - E) b_{A}^{†} b_{I} |Ψ_{0}⟩ . \end{matrix}

(8)

It can be written in the matrix form as

\begin{matrix} Δ E^{(2)} & = & \frac{1}{2} [\begin{matrix} C^{†} & C^{T} & D^{†} & D^{T} \end{matrix}] [\begin{matrix} A^{e} & B^{e} & T^{*} & R^{*} \\ B^{e *} & A^{e *} & R & T \\ T^{T} & R^{†} & A^{p} & B^{p} \\ R^{T} & T^{†} & B^{p *} & A^{p *} \end{matrix}] [\begin{matrix} C \\ C^{*} \\ D \\ D^{*} \end{matrix}] \\ \equiv & \frac{1}{2} u^{†} M u, \end{matrix}

(9)

with the matrix elements defined as

C_{i a} = c_{i}^{a},

(10)

D_{I A} = d_{I}^{A},

(11)

A_{i a, j b}^{e} = ⟨ Ψ_{0} | a_{i}^{†} a_{a} (Ĥ - E_{0}) a_{b}^{†} a_{j} | Ψ_{0} ⟩,

(12)

B_{i a, j b}^{e} = ⟨ Ψ_{0} | a_{i}^{†} a_{a} a_{j}^{†} a_{b} (Ĥ - E_{0}) | Ψ_{0} ⟩,

(13)

A_{I A, J B}^{p} = ⟨ Ψ_{0} | b_{I}^{†} b_{A} (Ĥ - E_{0}) b_{B}^{†} b_{J} | Ψ_{0} ⟩,

(14)

B_{I A, J B}^{p} = ⟨ Ψ_{0} | b_{I}^{†} b_{A} b_{J}^{†} b_{B} (Ĥ - E_{0}) | Ψ_{0} ⟩,

(15)

R_{i a, J B} = ⟨ Ψ_{0} | (Ĥ - E_{0}) b_{B}^{†} b_{J} a_{a}^{†} a_{i} | Ψ_{0} ⟩,

(16)

T_{i a, J B} = ⟨ Ψ_{0} | b_{J}^{†} b_{B} (Ĥ - E_{0}) a_{a}^{†} a_{i} | Ψ_{0} ⟩ .

(17)

The matrix M is the stability matrix, or the Hessian. It is a Hermitian matrix and therefore has only real eigenvalues. Evaluation of the matrix elements with the Hamiltonian in Eq. (5) using canonical SCF orbitals gives

A_{i a, j b}^{e} = δ_{i j} δ_{a b} (ε_{a} - ε_{i}) + ⟨ a j | | i b ⟩,

(18)

B_{i a, j b}^{e} = ⟨ a b | | i j ⟩,

(19)

A_{I A, J B}^{p} = δ_{I J} δ_{A B} (ε_{A} - ε_{I}) + ⟨ A J | | I B ⟩,

(20)

B_{I A, J B}^{p} = ⟨ A B | | I J ⟩,

(21)

R_{i a, J B} = - ⟨ i J | a B ⟩,

(22)

T_{i a, J B} = - ⟨ i B | a J ⟩,

(23)

where δ is the Kronecker delta and the ε’s with associated subscripts are the NEO-HF orbital energies. For both electron and proton orbitals,

⟨ p q | | r s ⟩ = ⟨ p q | r s ⟩ - ⟨ p q | s r ⟩

(24)

and

⟨ p q | r s ⟩ = \int d x d x^{'} \frac{ϕ_{p}^{*} (x) ϕ_{q}^{*} (x^{'}) ϕ_{r} (x) ϕ_{s} (x^{'})}{| r - r^{'} |},

(25)

where x is the combined variable of spatial and spin coordinates (r,σ) and the ϕ’s with associated subscripts are the NEO-HF spin orbitals. Note that the stability matrix M in Eq. (9) is the same as the NEO-TDHF matrix, although NEO-TDHF has a non-identity metric in its working equation.¹² Also note that if the matrix elements are evaluated with non-canonical orbitals, the elements in the two A matrices must be modified to be $A_{i a, j b}^{e} = ⟨ j | i ⟩ f_{a b}^{e} - ⟨ a | b ⟩ f_{j i}^{e} + ⟨ a j | | i b ⟩$ and $A_{I A, J B}^{p} = ⟨ J | I ⟩ f_{A B}^{p} - ⟨ A | B ⟩ f_{J I}^{p} + ⟨ A J | | I B ⟩,$ where f^e and f^p are the Fock matrices for the electrons and protons, respectively. However, the form of the elements in the other submatrices remains the same for non-canonical orbitals.

An SCF solution is stable if ΔE⁽²⁾ is never negative for any orbital rotation [i.e., any trial vector u in Eq. (9)]. This requirement is equivalent to the condition that the stability matrix M always has non-negative eigenvalues, or in other words, M is positive semi-definite. The detailed derivations of these equations are provided in the supplementary material.

B. NEO-DFT stability conditions

In NEO-DFT, the energy of a system can be expressed as

\begin{matrix} E [ρ^{e}, ρ^{p}] & = & E^{e} [ρ^{e}; ρ^{p}] + E^{p} [ρ^{p}; ρ^{e}] \\ - J^{e p} [ρ^{e}, ρ^{p}] - E_{c}^{e p} [ρ^{e}, ρ^{p}], \end{matrix}

(26)

\begin{matrix} E^{e} [ρ^{e}; ρ^{p}] & = & T_{s}^{e} [ρ^{e}] + J^{e e} [ρ^{e}] + E_{x c}^{e} [ρ^{e}] + V_{e x t}^{e} [ρ^{e}] \\ + J^{e p} [ρ^{e}, ρ^{p}] + E_{c}^{e p} [ρ^{e}, ρ^{p}], \end{matrix}

(27)

\begin{matrix} E^{p} [ρ^{p}; ρ^{e}] & = & T_{s}^{p} [ρ^{p}] + J^{p p} [ρ^{p}] + E_{x c}^{p} [ρ^{p}] + V_{e x t}^{p} [ρ^{p}] \\ + J^{e p} [ρ^{e}, ρ^{p}] + E_{c}^{e p} [ρ^{e}, ρ^{p}] . \end{matrix}

(28)

Here the electronic energy E^e is the same as that of regular DFT except for the addition of the mean-field electron-proton Coulomb energy J^ep, where both the electrons and the quantum protons are represented by densities, and the electron-proton correlation energy $E_{c}^{e p}$ ⁠. The electron-proton Coulomb energy J^ep and the electron-proton correlation energy $E_{c}^{e p}$ are subtracted in the total energy expression to eliminate double counting.

Under the orbital rotations in Eqs. (2) and (3), the changes in the non-interacting kinetic energy can be evaluated in a similar manner as for NEO-HF. The changes in the rest of the terms explicitly depend on the changes of the densities, which can be evaluated with

Δ ρ^{e} (r) = \frac{⟨ Ψ_{0}^{'} | {\hat{ρ}}^{e} (r) | Ψ_{0}^{'} ⟩}{⟨ Ψ_{0}^{'} | Ψ_{0}^{'} ⟩} - \frac{⟨ Ψ_{0} | {\hat{ρ}}^{e} (r) | Ψ_{0} ⟩}{⟨ Ψ_{0} | Ψ_{0} ⟩},

(29)

Δ ρ^{p} (r) = \frac{⟨ Ψ_{0}^{'} | {\hat{ρ}}^{p} (r) | Ψ_{0}^{'} ⟩}{⟨ Ψ_{0}^{'} | Ψ_{0}^{'} ⟩} - \frac{⟨ Ψ_{0} | {\hat{ρ}}^{p} (r) | Ψ_{0} ⟩}{⟨ Ψ_{0} | Ψ_{0} ⟩},

(30)

where ${\hat{ρ}}^{e}$ and ${\hat{ρ}}^{p}$ are the density operators for electrons and protons, respectively. For a specific term V[ρ^e, ρ^p] in the energy expression given by Eqs. (26)–(28), the change in V, up to second order, is

\begin{matrix} Δ V & = & \int d r \{) V [ρ^{e} (r) + Δ ρ^{e} (r), ρ^{p} (r) + Δ ρ^{p} (r)] \\ - V [ρ^{e} (r), ρ^{p} (r)] \}) \\ \approx & \int d r [\frac{δ V}{δ ρ^{e} (r)} Δ ρ^{e} (r) + \frac{δ V}{δ ρ^{p} (r)} Δ ρ^{p} (r)] \\ + \int d r d r^{'} [) \frac{1}{2} \frac{δ^{2} V}{δ ρ^{e} (r) δ ρ^{e} (r^{'})} Δ ρ^{e} (r) Δ ρ^{e} (r^{'}) \\ + \frac{1}{2} \frac{δ^{2} V}{δ ρ^{p} (r) δ ρ^{p} (r^{'})} Δ ρ^{p} (r) Δ ρ^{p} (r^{'}) \\ + \frac{δ^{2} V}{δ ρ^{e} (r) δ ρ^{p} (r^{'})} Δ ρ^{e} (r) Δ ρ^{p} (r^{'})]) . \end{matrix}

(31)

The first-order change in the total energy is

Δ E^{(1)} = \sum_{i a} c_{i}^{a} f_{a i}^{e} + c c + \sum_{I A} d_{I}^{A} f_{A I}^{p} + c c,

(32)

where f^e (f^p) is the DFT Fock operator for electrons (protons) related to the energy expression E^e (E^p). These first-order terms vanish when an SCF solution is achieved because the solutions of the Hermitian Fock operator are orthogonal. The second-order changes can be cast into the same matrix form as Eq. (9), but with slightly modified definitions of the matrix elements,

\begin{matrix} A_{i a, j b}^{e} & = & δ_{i j} δ_{a b} (ε_{a} - ε_{i}) + ⟨ a j | i b ⟩ + ⟨ a j | \frac{δ^{2} (E_{x c}^{e} + E_{c}^{e p})}{δ ρ^{e} δ ρ^{e}} | i b ⟩, \end{matrix}

(33)

B_{i a, j b}^{e} = ⟨ a b | i j ⟩ + ⟨ a b | \frac{δ^{2} (E_{x c}^{e} + E_{c}^{e p})}{δ ρ^{e} δ ρ^{e}} | i j ⟩,

(34)

\begin{matrix} A_{I A, J B}^{p} & = & δ_{I J} δ_{A B} (ε_{A} - ε_{I}) + ⟨ A J | I B ⟩ + ⟨ A J | \frac{δ^{2} (E_{x c}^{p} + E_{c}^{e p})}{δ ρ^{p} δ ρ^{p}} | I B ⟩, \end{matrix}

(35)

B_{I A, J B}^{p} = ⟨ A B | I J ⟩ + ⟨ A B | \frac{δ^{2} (E_{x c}^{p} + E_{c}^{e p})}{δ ρ^{p} δ ρ^{p}} | I J ⟩,

(36)

R_{i a, J B} = - ⟨ i J | a B ⟩ + ⟨ i J | \frac{δ^{2} E_{c}^{e p}}{δ ρ^{e} δ ρ^{p}} | a B ⟩,

(37)

T_{i a, J B} = - ⟨ i B | a J ⟩ + ⟨ i B | \frac{δ^{2} E_{c}^{e p}}{δ ρ^{e} δ ρ^{p}} | a J ⟩,

(38)

where the integral associated with a given (exchange-)correlation kernel K is defined as

\begin{matrix} ⟨ p q | K | r s ⟩ = \int d x d x^{'} ϕ_{p}^{*} (x) ϕ_{q}^{*} (x^{'}) K (x, x^{'}) ϕ_{r} (x) ϕ_{s} (x^{'}) . \end{matrix}

(39)

Therefore, the NEO-DFT stability matrix M is the same as the previously derived NEO-TDDFT matrix, but NEO-TDDFT has a non-identity metric in its working equation.¹² A NEO-DFT solution is stable when the stability matrix M is positive semi-definite. The detailed derivations of these equations are provided in the supplementary material.

III. CONNECTION OF STABILITY ANALYSIS TO NEO-TDHF AND NEO-TDDFT IMAGINARY EIGENVALUES

In this section, we discuss the connection between the NEO-HF stability analysis and NEO-TDHF, with the understanding that an analogous connection between the NEO-DFT stability analysis and NEO-TDDFT also exists. An imaginary eigenvalue within NEO-TDHF theory is direct evidence of the instability of an underlying SCF solution. Following the same logic as Thouless,³³ this statement can be proven by contradiction. Assume NEO-TDHF has at least one imaginary eigenvalue, but the underlying NEO-HF solution is stable. Suppose one NEO-TDHF root with the imaginary eigenvalue E_n satisfies

\begin{matrix} [\begin{matrix} A^{e} & B^{e} & T^{*} & R^{*} \\ B^{e *} & A^{e *} & R & T \\ T^{T} & R^{†} & A^{p} & B^{p} \\ R^{T} & T^{†} & B^{p *} & A^{p *} \end{matrix}] [\begin{matrix} X_{n}^{e} \\ Y_{n}^{e} \\ X_{n}^{p} \\ Y_{n}^{p} \end{matrix}] = E_{n} [\begin{matrix} I & 0 & 0 & 0 \\ 0 & - I & 0 & 0 \\ 0 & 0 & I & 0 \\ 0 & 0 & 0 & - I \end{matrix}] [\begin{matrix} X_{n}^{e} \\ Y_{n}^{e} \\ X_{n}^{p} \\ Y_{n}^{p} \end{matrix}], \end{matrix}

(40)

where X and Y represent changes in the density matrices.¹² Multiplying by $[X_{n}^{e^{†}} Y_{n}^{e^{†}} X_{n}^{p^{†}} Y_{n}^{p^{†}}]$ on the left of both sides gives

\begin{matrix} [X_{n}^{e^{†}} Y_{n}^{e^{†}} X_{n}^{p^{†}} Y_{n}^{p^{†}}] [\begin{matrix} A^{e} & B^{e} & T^{*} & R^{*} \\ B^{e *} & A^{e *} & R & T \\ T^{T} & R^{†} & A^{p} & B^{p} \\ R^{T} & T^{†} & B^{p *} & A^{p *} \end{matrix}] [\begin{matrix} X_{n}^{e} \\ Y_{n}^{e} \\ X_{n}^{p} \\ Y_{n}^{p} \end{matrix}] \\ = E_{n} (X_{n}^{e^{†}} X_{n}^{e} - Y_{n}^{e^{†}} Y_{n}^{e} + X_{n}^{p^{†}} X_{n}^{p} - Y_{n}^{p^{†}} Y_{n}^{p}) . \end{matrix}

(41)

Because the stability matrix is Hermitian, the left-hand side of the equation is real. Because E_n is imaginary, the real number $(X_{n}^{e^{†}} X_{n}^{e} - Y_{n}^{e^{†}} Y_{n}^{e} + X_{n}^{p^{†}} X_{n}^{p} - Y_{n}^{p^{†}} Y_{n}^{p})$ must be zero and both sides of the equation are zero. Since we have assumed that the NEO-HF solution is stable and therefore has non-negative eigenvalues, according to the variation principle, the vector that yields expectation value zero is an eigenvector with eigenvalue zero,

\begin{matrix} [\begin{matrix} A^{e} & B^{e} & T^{*} & R^{*} \\ B^{e *} & A^{e *} & R & T \\ T^{T} & R^{†} & A^{p} & B^{p} \\ R^{T} & T^{†} & B^{p *} & A^{p *} \end{matrix}] [\begin{matrix} X_{n}^{e} \\ Y_{n}^{e} \\ X_{n}^{p} \\ Y_{n}^{p} \end{matrix}] = 0 [\begin{matrix} X_{n}^{e} \\ Y_{n}^{e} \\ X_{n}^{p} \\ Y_{n}^{p} \end{matrix}] = 0 [\begin{matrix} I & 0 & 0 & 0 \\ 0 & - I & 0 & 0 \\ 0 & 0 & I & 0 \\ 0 & 0 & 0 & - I \end{matrix}] [\begin{matrix} X_{n}^{e} \\ Y_{n}^{e} \\ X_{n}^{p} \\ Y_{n}^{p} \end{matrix}] . \end{matrix}

(42)

As a consequence, E_n must be zero and real, which contradicts the assumption that it is imaginary. Therefore, an imaginary eigenvalue in NEO-TDHF is a sufficient condition for the instability of an underlying NEO-HF solution.

However, an unstable SCF solution with a negative eigenvalue for the Hessian does not necessarily lead to an imaginary eigenvalue in NEO-TDHF theory. A simple example is to assume that $A^{e} = A^{e *} = [\begin{matrix} 1 & 2 \\ 2 & 1 \end{matrix}]$ and all other submatrices are 0. In this case, the stability matrix has a negative eigenvalue −1, but the NEO-TDHF working equation has only real eigenvalue solutions of ±1 and ±3. We will also present a practical example of this situation for a molecular system in Sec. VI.

IV. ELECTRONIC, PROTONIC, AND VIBRONIC INSTABILITIES

If we set the proton orbital rotation to be zero (i.e., D = D^* = 0) and only allow the electron orbitals to vary, the second-order energy change reduces to

\begin{matrix} Δ E^{(2)} = \frac{1}{2} [\begin{matrix} C^{†} & C^{T} \end{matrix}] [\begin{matrix} A^{e} & B^{e} \\ B^{e *} & A^{e *} \end{matrix}] [\begin{matrix} C \\ C^{*} \end{matrix}] \equiv \frac{1}{2} u^{e †} M^{e} u^{e} . \end{matrix}

(43)

This form is the same as that of the second-order energy change in pure electronic HF or DFT,^21,25 although the proton quantization impacts the electronic orbitals in the NEO framework. Moreover, in NEO-DFT, as shown in Eqs. (33) and (34), the second derivative of the electron-proton correlation functional will contribute to the A^e and B^e blocks. If the electronic matrix M^e is unstable with a negative eigenvalue, the full NEO matrix M must be unstable, which we define to be an electronic instability. Similarly, we define a protonic instability when the submatrix $[\begin{matrix} A^{p} & B^{p} \\ B^{p *} & A^{p *} \end{matrix}]$ has a negative eigenvalue.

When there is no electronic or protonic instability, the NEO solution is still not guaranteed to be stable because the matrices R and T that couple the electron and proton parts can lead to instabilities. In those cases, the system will be stable with respect to both pure electron orbital rotations and pure proton orbital rotations, but it is unstable with respect to coupled electron-proton orbital rotations, which is defined to be an electron-proton vibronic instability. Although mathematically well-defined, we have not yet observed this type of vibronic instability in practice.

V. STABILITY THEORY WITH SPIN AND SPATIAL CONSTRAINTS

Analogous to conventional electronic HF theory, NEO-HF can be classified into generalized NEO-HF (NEO-GHF) with the allowance of mixed-spin orbitals, unrestricted NEO-HF (NEO-UHF) with pure-spin orbitals but not necessarily the same set of spatial orbitals for different spins, and restricted NEO-HF (NEO-RHF) with the same set of spatial orbitals shared by different spins. Note that while NEO-GHF has never been implemented, the derivation of the stability analysis in Sec. II is general and also applies to NEO-GHF. Additionally, because NEO treats more than one type of particle quantum mechanically, it is possible to treat one type of particle, such as electrons, with RHF, while treating the other type of particle, such as protons, with UHF. Practical NEO calculations often utilize this mixed treatment, which is denoted as NEO-MHF in this section. Furthermore, for each type of theory, the orbitals may or may not be constrained to be real. The relationships among these methods are depicted schematically in Fig. 1. Following the work by Seeger and Pople,²³ we will discuss the stability conditions for real and complex NEO-GHF, NEO-UHF, NEO-RHF, and NEO-MHF. The results are summarized in the main paper, and the detailed derivations are given in the supplementary material. The stability conditions for the various types of NEO-DFT are similar, and we will not discuss them separately.

FIG. 1.

View large Download slide

Schematic depiction of spatial and spin constraints in NEO-HF theory, starting from the maximum constraints at the top and the most general case at the bottom. The solid black lines correspond to the same relationships as found in conventional electronic structure theory,²³ while the dashed lines indicate new relationships arising from the possibility of MHF. Note that arrows indicate only relationships between adjacent levels.

A. Complex NEO-GHF

For the most general complex NEO-GHF case, the stability can be evaluated by directly diagonalizing the stability matrix M in Eq. (9).

B. Real NEO-GHF

For the real NEO-GHF case, the stability matrix M is real and symmetric, and the second-order energy change can be rewritten as

\begin{matrix} Δ E^{(2)} & = & \frac{1}{2} [\begin{matrix} C^{†} & C^{T} & D^{†} & D^{T} \end{matrix}] [\begin{matrix} A^{e} & B^{e} & T & R \\ B^{e} & A^{e} & R & T \\ T^{T} & R^{T} & A^{e} & B^{e} \\ R^{T} & T^{T} & B^{p} & A^{p} \end{matrix}] [\begin{matrix} C \\ C^{*} \\ D \\ D^{*} \end{matrix}] \\ = & \frac{1}{2} [\begin{matrix} \frac{1}{\sqrt{2}} (C^{†} + C^{T}) & \frac{1}{\sqrt{2}} (D^{†} + D^{T}) \end{matrix}] [\begin{matrix} A^{e} + B^{e} & T + R \\ T^{T} + R^{T} & A^{p} + B^{p} \end{matrix}] [\begin{matrix} \frac{1}{\sqrt{2}} (C + C^{*}) \\ \frac{1}{\sqrt{2}} (D + D^{*}) \end{matrix}] \\ + \frac{1}{2} [\begin{matrix} \frac{1}{\sqrt{2}} (C^{†} - C^{T}) & \frac{1}{\sqrt{2}} (D^{†} - D^{T}) \end{matrix}] [\begin{matrix} A^{e} - B^{e} & T - R \\ T^{T} - R^{T} & A^{p} - B^{p} \end{matrix}] [\begin{matrix} \frac{1}{\sqrt{2}} (C - C^{*}) \\ \frac{1}{\sqrt{2}} (D - D^{*}) \end{matrix}] \\ = & \frac{1}{2} v_{1}^{†} N_{1} v_{1} + \frac{1}{2} v_{2}^{†} N_{2} v_{2} . \end{matrix}

(44)

The vector v₁ is always real, and therefore the corresponding matrix N₁ characterizes the internal stability within real NEO-GHF. When N₁ is not positive semi-definite, there is an instability within the real NEO-GHF space. The vector v₂ is always imaginary, and the corresponding matrix N₂ characterizes the external stability with respect to complex NEO-GHF. When N₂ is not positive semi-definite, a lower energy can be achieved by allowing the orbitals to be complex. Furthermore, T = R and T^T = R^T when the orbitals are real, and therefore the electron and proton sub-blocks in N₂ are decoupled, with

\frac{1}{2} v_{2}^{†} N_{2} v_{2} = \frac{1}{2} v_{2}^{e †} N_{2}^{e} v_{2}^{e} + \frac{1}{2} v_{2}^{p †} N_{2}^{p} v_{2}^{p} .

(45)

The matrices $N_{2}^{e}$ and $N_{2}^{p}$ characterize external electronic and protonic stabilities, respectively, with respect to complex NEO-GHF.

C. Complex NEO-UHF

For NEO-UHF, the orbital rotations can be classified as spin-conserved (SC) rotations, which mix orbitals with the same spins, and spin-flip (SF) rotations, which mix orbitals with different spins. For the complex NEO-UHF case, the SC and SF parts can be decoupled for each sub-block of the stability matrix M, and the second-order energy change can be written as

\begin{matrix} Δ E^{(2)} = \frac{1}{2} u_{S C}^{†} M_{S C} u_{S C} + \frac{1}{2} u_{S F}^{e †} M_{S F}^{e} u_{S F}^{e} + \frac{1}{2} u_{S F}^{p †} M_{S F}^{p} u_{S F}^{p} . \end{matrix}

(46)

The SC matrix couples electron and protons, whereas the SF matrices only contain either pure electron or pure proton sub-blocks. This distinction arises because the SF coupling is non-zero only in the Fock exchange kernel or non-collinear density functional kernels,³⁴ but there are no interactions of these types between electrons and protons. The SC matrix M_SC characterizes the internal stability within complex NEO-UHF, whereas the SF matrices $M_{S F}^{e}$ and $M_{S F}^{p}$ characterize the external stability with respect to complex NEO-GHF, which we refer to as SF stability.

D. Real NEO-UHF

For the real NEO-UHF case, the stability matrices can be obtained from those of complex NEO-UHF with the same technique as used in Eq. (44). The SC matrix N_1,SC characterizes the internal stability within real NEO-UHF. The SF matrices $N_{1, S F}^{e}$ and $N_{1, S F}^{p}$ characterize the external stability with respect to real NEO-GHF. The SC matrices $N_{2, S C}^{e}$ and $N_{2, S C}^{p}$ characterize the external stability with respect to complex NEO-UHF. The SF matrices $N_{2, S F}^{e}$ and $N_{2, S F}^{p}$ characterize the external stability with respect to complex NEO-GHF.

E. Complex NEO-RHF

The complex NEO-RHF case, in which both electrons and protons are treated as closed-shell systems, is rarely used for molecular systems because proton orbitals are highly localized compared to electronic orbitals, and the Coulomb repulsion energy would be extremely high if two protons formed a pair and occupied the same spatial orbital. Nevertheless, for completeness, we still analyze its stability conditions from a purely mathematical perspective. For the complex NEO-RHF case, the spatial variations for α and β electrons and protons are required to be the same, and therefore we define the new quantities

{}^{1}C = \frac{1}{\sqrt{2}} (C_{α α} + C_{β β}),

(47)

{}^{3}C = \frac{1}{\sqrt{2}} (C_{α α} - C_{β β}),

(48)

{}^{1}D = \frac{1}{\sqrt{2}} (D_{α α} + D_{β β}),

(49)

{}^{3}D = \frac{1}{\sqrt{2}} (D_{α α} - D_{β β}) .

(50)

In this case, ¹C (¹D) corresponds to a rotation of the electron (proton) orbitals within the complex NEO-RHF space, while ³C (³D) corresponds to a rotation of the electron (proton) orbitals outside the complex NEO-RHF space. In other words, ¹C (¹D) denotes singlet rotations, while ³C (³D) denotes triplet rotations. As a consequence, the SC stability matrix M_SC factorizes into a singlet one (¹M) and a triplet one (³M). Furthermore, the electron part (³M^e) and proton part (³M^p) do not couple in the triplet matrix. The singlet matrix ¹M characterizes the internal stability within complex NEO-RHF, while the triplet matrices ³M^e and ³M^p characterize the external stability with respect to complex NEO-UHF. Furthermore, the same SF stability as defined for complex NEO-UHF with respect to complex NEO-GHF also applies here.

F. Real NEO-RHF

For the real NEO-RHF case, the same technique as used in Eq. (44) can be applied to the complex NEO-RHF matrices to obtain the stability matrices. Starting from ¹M, the internal stability within real NEO-RHF is characterized by the first singlet matrix ¹N₁, and the external stability with respect to complex NEO-RHF is characterized by the second singlet matrix ¹N₂. The same technique can be applied to the triplet matrix ³M, and the external stability with respect to real NEO-UHF (or complex NEO-UHF) is characterized by the triplet matrix ³N₁ (or ³N₂) Again, the electron and proton parts decouple for all matrices except the matrix ¹N₁. Analogous to real NEO-UHF, SF stability matrices can be defined to characterize external stabilities with respect to real and complex NEO-GHF, as discussed in the supplementary material.

G. Complex NEO-MHF

For complex NEO-MHF, in which electrons are constrained to be closed-shell singlet while protons are assumed to be all high-spin (i.e., each proton occupies a different spatial orbital), the singlet and triplet quantities are defined for electrons only. The singlet electronic part couples to the proton rotation and gives rise to a matrix ^aM, while the triplet electronic part does not couple to the proton rotation and is denoted ^bM. Therefore, the internal stability within NEO-MHF is characterized by the coupled matrix ^aM, whereas the external stability with respect to complex NEO-UHF is characterized by the pure electronic ^bM. Similar to complex NEO-RHF, the SF stability with respect to complex NEO-GHF also applies here.

H. Real NEO-MHF

For the real NEO-MHF case, the same technique as given by Eq. (44) can be applied to the complex NEO-MHF stability matrices, as described in the supplementary material. For the remainder of this paper, we will assume that NEO-HF corresponds to real NEO-MHF for notational simplicity.

VI. STABILITY ANALYSIS FOR MOLECULAR SYSTEMS

In this section, we present some examples of instabilities from NEO-DFT calculations on molecular systems. Our implementation for these molecules, which have one quantum proton, includes only the SC block for the proton because the SF block is not of interest. First we investigated the FHF⁻ molecule, where the proton and all electrons are treated quantum mechanically. The geometry was determined from a geometry optimization at the conventional electronic DFT/B3LYP level with the def2-QZVP electronic basis set. For the NEO calculations, the cc-pVDZ electronic basis set and the 8s8p8d nuclear basis set¹⁰ were used. Although a larger electronic basis set with diffuse functions may improve the quantitative accuracy, this basis set is sufficient for illustrative purposes in the context of the NEO-DFT stability analysis. The electronic and nuclear basis function centers associated with the quantum hydrogen were placed at the midpoint between the two fluorine atoms, which were separated by 2.30 Å.

For the epc17-1 and epc17-2 functionals, the FHF⁻ molecule exhibits an instability in the SCF solutions with energies −199.791 a.u. and −199.757 a.u., respectively. By contrast, the lowest-energy SCF solutions found for epc17-1 and epc17-2 correspond to energies −200.260 a.u. and −200.235 a.u., respectively. For each of the higher-energy solutions, the NEO-DFT SCF equations are satisfied, and therefore the gradient in the orbital space is zero (i.e., they are stationary points). The instability indicates that the higher-energy solutions correspond to saddle points rather than local minima in the orbital space. Further investigation suggests that they are electronically stable but not protonically stable because diagonalization of the submatrix $[\begin{matrix} A^{e} & B^{e} \\ B^{e *} & A^{e *} \end{matrix}]$ yields all positive eigenvalues, while diagonalization of the submatrix $[\begin{matrix} A^{p} & B^{p} \\ B^{p *} & A^{p *} \end{matrix}]$ yields negative eigenvalues. Moreover, both higher-energy solutions correspond to internal NEO-HF instabilities (i.e., within the set of solutions with closed-shell RHF electrons and a single proton). Analysis of the proton density, in conjunction with the extremely high energy, suggests that these solutions do not correspond to excited states of interest.

Despite the instabilities identified in the SCF solutions, the corresponding NEO-TDDFT working equations were found to have all real eigenvalue solutions. This example illustrates that an instability of the SCF solution does not necessarily imply imaginary eigenvalues for the NEO-TDDFT working equation. Additionally, we have also been able to find examples of the more common case, in which negative eigenvalues of the stability matrix and imaginary eigenvalues of the NEO-TDDFT working equation both exist for the HCN molecule (see the supplementary material for details). Note that the higher-energy SCF solutions discussed above are not necessarily physically meaningful and only serve an illustrative purpose. However, in some cases, such as the conversion between HCN and HNC and the hydrogen transfer in 2-cyanomalonaldehyde presented below, the higher-energy solutions are physically meaningful.

We also investigated the two isomers HCN and HNC, where again the proton and all electrons are treated quantum mechanically in the NEO calculations. The schematic picture of the two isomers is shown in Fig. 2. Geometry optimizations were performed at the conventional electronic DFT/B3LYP level with the def2-QZVP electronic basis set for both HCN and HNC, and the resulting C–N bond lengths were averaged for a final C–N distance of 1.155 Å. Then geometry optimizations were performed for HCN and HNC to determine the location of the H atom with this fixed C–N distance. Electronic and nuclear basis function centers associated with the quantum hydrogen were placed at these two positions: one set of basis functions was placed at a distance of 0.996 Å from the nitrogen atom and another set was placed at a distance of 1.066 Å from the carbon atom. NEO calculations using the cc-pVDZ electronic basis set and the 8s8p8d nuclear basis set were converged to both the HCN and HNC solutions by starting with an initial guess localized on one side or the other. This type of calculation may be described as a NEO ∆SCF calculation.

FIG. 2.

View large Download slide

Schematic picture of the HCN and HNC molecules. One set of electron and proton basis functions for hydrogen is positioned at each end of the molecule, as depicted by the solid and dashed oblong shapes.

For the epc17-1 and epc17-2 functionals, the stability matrices for both the HCN and HNC solutions have only positive eigenvalues, indicating that they are minima in the orbital space. The total energy of HNC is 0.025 a.u. higher than the total energy of HCN for both the epc17-1 and epc17-2 functionals. This energy difference is consistent with the energy difference between the isomers HNC and HCN in conventional electronic DFT. It is also similar to the energy difference of 0.026 a.u. obtained when the hydrogen nucleus is treated quantum mechanically with a three-dimensional Fourier grid Hamiltonian method^35,36 at the B3LYP/cc-pVDZ level of theory. Thus, the HNC solution is a local minimum in the orbital space, while the HCN solution is presumed to be the global minimum because a lower-energy solution was not found. As expected, for both SCF solutions, the NEO-TDDFT working equations yield real eigenvalue solutions. This example of HNC and HCN isomers can be viewed as a special case of an asymmetric double well potential for the hydrogen with two well-separated minima connected by a high barrier.

A more typical example of an asymmetric double-well proton transfer system is 2-cyanomalonaldehyde, where the proton transfers between the two oxygen atoms (Fig. 3). In this case, the transferring proton and all electrons are treated quantum mechanically in the NEO calculations. For illustrative purposes, we utilized a geometry for malonaldehyde obtained previously³⁷ by averaging the reactant and product geometries that were optimized at the MP2/6-31G(d,p) level. Starting with this geometry for malonaldehyde, the molecular symmetry was broken by replacing a hydrogen with a CN group, followed by optimizing only the coordinates of the CN group at the B3LYP/6-31+G(d,p) level. The NEO calculations used two electronic and nuclear basis function centers for the quantum hydrogen, corresponding to the minima of the left and right wells, respectively, of the asymmetric double well proton transfer potential. The positions of these basis function centers were determined by optimizing the hydrogen on the left or right oxygen at the conventional DFT/B3LYP/6-31+G(d,p) level. The cc-pVQZ electronic basis set and 8s8p8d nuclear basis set were used for each basis function center associated with the quantum proton, and the cc-pVDZ electronic basis set was used for the other nuclei. The stationary solutions localized in the left well and the right well were obtained by starting with an initial guess localized on one side or the other, corresponding to a NEO ΔSCF calculation for the higher-energy solution. As a reference, the hydrogen vibrational wave functions and energy levels were obtained using a three-dimensional Fourier grid Hamiltonian method^35,36 at the DFT/B3LYP/6-31+G(d,p) level of theory for the fixed coordinates of the other nuclei.

FIG. 3.

View large Download slide

The upper part shows the structures for 2-cyanomalonaldehyde obtained as described in the text, where all nuclei except the transferring hydrogen are fixed to an average reactant/product geometry. The lower part shows the one-dimensional slice of the proton potential energy surface along the proton transfer axis, which is defined by the line connecting the optimized positions of the transferring hydrogen on each oxygen. The lowest two proton vibrational energy levels and the corresponding one-dimensional slices of the proton vibrational wave functions are also depicted. The energies corresponding to the barrier height, difference in minima, and splitting between the lowest two vibrational energy levels computed with the Fourier grid Hamiltonian method are given, and the NEO-DFT epc17-2 energy level splitting is given in parentheses. The coordinates of the two structures shown are provided in the supplementary material.

For the epc17-1 and epc17-2 functionals, the stability matrices for both localized solutions have only positive eigenvalues, indicating that they are minima in the NEO orbital space. The epc17-2 solution localized in the right well is 575 cm⁻¹ (1.64 kcal/mol) higher in energy than the solution localized in the left well. This energy splitting is almost identical to the grid-based reference value, with a difference of only 48 cm⁻¹ (0.14 kcal/mol). This splitting is also similar to the energy difference of the two minimum energy geometries obtained by conventional DFT, suggesting that the zero-point energies associated with the proton in each well are very similar. Furthermore, NEO-TDDFT calculations starting from each of the two localized NEO-DFT solutions give only real eigenvalues that correspond to transitions to higher proton vibrational excited states. These excited vibrational states are predominantly localized within the same well as the initial reference and agree reasonably well with the grid-based reference values for the lowest five vibrational states, as given in the supplementary material. Delocalized solutions with even higher energies are more challenging to obtain and will be investigated in the future.

In conventional electronic DFT, HCN and HNC, as well as the two proton positions in 2-cyanomalonaldehyde, correspond to two different isomers that are both minima in the nuclear coordinate space. In NEO-DFT, where specified protons are treated quantum mechanically on the same level as the electrons, all of these isomers correspond to minima in the electron/proton orbital space. The mathematical relationship between stationary points in nuclear coordinate space for conventional electronic DFT and orbital space for NEO-DFT is not clear. However, these two examples illustrate that in some cases, a stationary point in orbital space of a NEO-DFT calculation corresponds to a stationary point in nuclear coordinate space of a conventional electronic DFT calculation. Moreover, in these cases, the stationary points from NEO-DFT correspond to a physically meaningful form of the molecule.

VII. CONCLUSIONS

We have derived the stability conditions for the NEO-HF and NEO-DFT methods. If the stability matrix is positive semi-definite, the solution is stable and corresponds to either the global minimum or a local minimum, whereas otherwise it is a saddle point within the orbital space. Internal and external stability matrices for complex and real NEO-GHF, NEO-UHF, NEO-RHF, and NEO-MHF (i.e., RHF for closed-shell electrons and UHF for high-spin protons) have been defined. Moreover, the multicomponent SCF framework allows the potential for electronic, protonic, and electron-proton vibronic instabilities. The stability matrices for NEO-HF and NEO-DFT are the same as those found in NEO-TDHF and NEO-TDDFT, but NEO-TDHF and NEO-TDDFT have non-identity metrics in their working equations. An imaginary eigenvalue solution for NEO-TDHF or NEO-TDDFT indicates that there must be at least one negative eigenvalue in the stability matrix. However, a negative eigenvalue in the stability matrix does not necessarily indicate an imaginary eigenvalue solution of the NEO-TDHF or NEO-TDDFT working equation. Thus, a negative eigenvalue of the stability matrix is a necessary but not sufficient condition for the corresponding NEO-TDHF or NEO-TDDFT working equation to have an imaginary eigenvalue solution.

We have found and characterized different types of SCF solutions that are stationary points in orbital space but are higher in energy than the presumed global minimum. For example, NEO-DFT calculations for the FHF⁻ and HCN molecules revealed instabilities corresponding to saddle points in orbital space. In particular, we found SCF solutions for which the eigenvalue of the stability matrix was negative but the eigenvalue solution of the NEO-TDDFT working equation was real, as well as the more typical unstable SCF solutions for which the eigenvalue solution of the NEO-TDDFT working equation was imaginary. In addition, NEO-DFT calculations using a ∆SCF procedure for the HCN molecule illustrated the possibility of a local minimum in orbital space, corresponding to the HNC molecule in this case. NEO-DFT calculations using this ∆SCF procedure for 2-cyanomalonaldehyde illustrated an analogous situation for two different proton transfer isomers.

This type of stability analysis is a useful tool for characterizing the nature of an SCF solution in orbital space. When instabilities are found, this analysis will also be helpful in the search for lower-energy solutions, following the same procedure as in electronic structure theory.^23,24 In the NEO framework, some of the stationary points corresponding to minima or saddle points in nuclear coordinate space in conventional electronic structure calculations may be converted to stationary points in orbital space. This analysis will be particularly interesting when the stationary points in nuclear coordinate space are distinguished by differences in the positions of the hydrogen nuclei, as illustrated by the HCN and HNC calculations described herein. Another example is a proton transfer system that exhibits two minima in nuclear coordinate space associated with the proton bound to the proton donor or acceptor, as illustrated by the 2-cyanomalonaldehyde system studied herein. These types of stationary points may also be stationary points in orbital space for NEO calculations and could correspond to either saddle points or minima in orbital space. The theory described herein will enable a detailed analysis of the topology of orbital space in NEO calculations and a comparison to the topology of nuclear coordinate space in conventional electronic structure calculations.

SUPPLEMENTARY MATERIAL

See supplementary material for the complete derivation of the stability analysis for NEO-HF and NEO-DFT, as well as the detailed formulation of internal and external instabilities for theories with different spatial and spin constraints. Coordinates for 2-cyanomalonaldehyde and excitation energies calculated with NEO-TDDFT for each localized solution with comparison to proton vibrational energy level splittings obtained from the grid-based approach are also included.

ACKNOWLEDGMENTS

This work was supported by the National Science Foundation Grant Nos. CHE-1361293 and CHE-1762018, as well as the National Science Foundation Partnerships for International Research and Education Grant No. 1545907. Some of this material (the cyanomalonaldehyde calculations) is based upon work supported by the U.S. Department of Energy, Office of Science, Offices of Basic Energy Sciences and Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program. We are grateful to Kurt Brorsen, Josh Goings, Fabijan Pavosevic, and Patrick Schneider for helpful discussions.

REFERENCES

1.

I. L.

Thomas

,

Phys. Rev.

185

,

90

(

1969

).

https://doi.org/10.1103/physrev.185.90

Crossref

2.

J. F.

Capitani

,

R. F.

Nalewajski

, and

R. G.

Parr

,

J. Chem. Phys.

76

,

568

(

1982

).

https://doi.org/10.1063/1.442703

Crossref

3.

M.

Tachikawa

,

K.

Mori

,

H.

Nakai

, and

K.

Iguchi

,

Chem. Phys. Lett.

290

,

437

(

1998

).

https://doi.org/10.1016/s0009-2614(98)00519-3

Crossref

4.

T.

Kreibich

and

E. K. U.

Gross

,

Phys. Rev. Lett.

86

,

2984

(

2001

).

https://doi.org/10.1103/physrevlett.86.2984

Crossref

PubMed

5.

S. P.

Webb

,

T.

Iordanov

, and

S.

Hammes-Schiffer

,

J. Chem. Phys.

117

,

4106

(

2002

).

https://doi.org/10.1063/1.1494980

Crossref

6.

C.

Swalina

,

M. V.

Pak

,

A.

Chakraborty

, and

S.

Hammes-Schiffer

,

J. Phys. Chem. A

110

,

9983

(

2006

).

https://doi.org/10.1021/jp0634297

Crossref

PubMed

7.

M. V.

Pak

,

A.

Chakraborty

, and

S.

Hammes-Schiffer

,

J. Phys. Chem. A

111

,

4522

(

2007

).

https://doi.org/10.1021/jp0704463

Crossref

PubMed

8.

H.

Nakai

,

Int. J. Quantum Chem.

107

,

2849

(

2007

).

https://doi.org/10.1002/qua.21379

Crossref

9.

A.

Sirjoosingh

,

M. V.

Pak

,

C.

Swalina

, and

S.

Hammes-Schiffer

,

J. Chem. Phys.

139

,

034102

(

2013

).

https://doi.org/10.1063/1.4812257

Crossref

PubMed

10.

Y.

Yang

,

K. R.

Brorsen

,

T.

Culpitt

,

M. V.

Pak

, and

S.

Hammes-Schiffer

,

J. Chem. Phys.

147

,

114113

(

2017

).

https://doi.org/10.1063/1.4996038

Crossref

PubMed

11.

K. R.

Brorsen

,

Y.

Yang

, and

S.

Hammes-Schiffer

,

J. Phys. Chem. Lett.

8

,

3488

(

2017

).

https://doi.org/10.1021/acs.jpclett.7b01442

Crossref

PubMed

12.

Y.

Yang

,

T.

Culpitt

, and

S.

Hammes-Schiffer

,

J. Phys. Chem. Lett.

9

,

1765

(

2018

).

https://doi.org/10.1021/acs.jpclett.8b00547

Crossref

PubMed

13.

S.

Bubin

,

M.

Pavanello

,

W.-C.

Tung

,

K. L.

Sharkey

, and

L.

Adamowicz

,

Chem. Rev.

113

,

36

(

2013

).

https://doi.org/10.1021/cr200419d

Crossref

PubMed

14.

S.

Hammes-Schiffer

and

A. A.

Stuchebrukhov

,

Chem. Rev.

110

,

6939

(

2010

).

https://doi.org/10.1021/cr1001436

Crossref

PubMed

15.

S.

Hammes-Schiffer

,

J. Am. Chem. Soc.

137

,

8860

(

2015

).

https://doi.org/10.1021/jacs.5b04087

Crossref

PubMed

16.

C.

Swalina

,

M. V.

Pak

, and

S.

Hammes-Schiffer

,

Chem. Phys. Lett.

404

,

394

(

2005

).

https://doi.org/10.1016/j.cplett.2005.01.115

Crossref

17.

A.

Chakraborty

,

M. V.

Pak

, and

S.

Hammes-Schiffer

,

Phys. Rev. Lett.

101

,

153001

(

2008

).

https://doi.org/10.1103/physrevlett.101.153001

Crossref

PubMed

18.

A.

Chakraborty

,

M. V.

Pak

, and

S.

Hammes-Schiffer

,

J. Chem. Phys.

131

,

124115

(

2009

).

https://doi.org/10.1063/1.3236844

Crossref

PubMed

19.

Y.

Liao

,

Acc. Chem. Res.

50

,

1956

(

2017

).

https://doi.org/10.1021/acs.accounts.7b00190

Crossref

PubMed

20.

P.

Goyal

and

S.

Hammes-Schiffer

,

ACS Energy Lett.

2

,

512

(

2017

).

https://doi.org/10.1021/acsenergylett.6b00723

Crossref

21.

D. J.

Thouless

,

Nucl. Phys.

21

,

225

(

1960

).

https://doi.org/10.1016/0029-5582(60)90048-1

Crossref

22.

J.

Čížek

and

J.

Paldus

,

J. Chem. Phys.

47

,

3976

(

1967

).

https://doi.org/10.1063/1.1701562

Crossref

23.

R.

Seeger

and

J. A.

Pople

,

J. Chem. Phys.

66

,

3045

(

1977

).

https://doi.org/10.1063/1.434318

Crossref

24.

H. B.

Schlegel

and

J. J. W.

McDouall

,

Computational Advanced Organic Chemistry Molecular Structure Reaction

(

Springer

,

Dordrecht

,

1991

), pp.

167

–

185

.

Google Scholar

25.

R.

Bauernschmitt

and

R.

Ahlrichs

,

J. Chem. Phys.

104

,

9047

(

1996

).

https://doi.org/10.1063/1.471637

Crossref

26.

Y.

Yamaguchi

,

I. L.

Alberts

,

J. D.

Goddard

, and

H. F.

Schaefer

,

Chem. Phys.

147

,

309

(

1990

).

https://doi.org/10.1016/0301-0104(90)85046-y

Crossref

27.

J. J.

Goings

,

F.

Ding

,

M. J.

Frisch

, and

X.

Li

,

J. Chem. Phys.

142

,

154109

(

2015

).

https://doi.org/10.1063/1.4918561

Crossref

PubMed

28.

T.

Yamada

and

S.

Hirata

,

J. Chem. Phys.

143

,

114112

(

2015

).

https://doi.org/10.1063/1.4929354

Crossref

PubMed

29.

U.

Bozkaya

,

J. M.

Turney

,

Y.

Yamaguchi

,

H. F.

Schaefer

, and

C. D.

Sherrill

,

J. Chem. Phys.

135

,

104103

(

2011

).

https://doi.org/10.1063/1.3631129

Crossref

PubMed

30.

U.

Bozkaya

,

J. Chem. Phys.

135

,

224103

(

2011

).

https://doi.org/10.1063/1.3665134

Crossref

PubMed

31.

A. C.

Vaucher

and

M.

Reiher

,

J. Chem. Theory Comput.

13

,

1219

(

2017

).

https://doi.org/10.1021/acs.jctc.7b00011

Crossref

PubMed

32.

T.

Ziegler

,

A.

Rauk

, and

E. J.

Baerends

,

Theor. Chim. Acta

43

,

261

(

1977

).

https://doi.org/10.1007/bf00551551

Crossref

33.

D. J.

Thouless

,

Nucl. Phys.

22

,

78

(

1961

).

https://doi.org/10.1016/0029-5582(61)90364-9

Crossref

34.

F.

Wang

and

T.

Ziegler

,

Int. J. Quantum Chem.

106

,

2545

(

2006

).

https://doi.org/10.1002/qua.21050

Crossref

35.

C. C.

Marston

and

G. G.

Balint-Kurti

,

J. Chem. Phys.

91

,

3571

(

1989

).

https://doi.org/10.1063/1.456888

Crossref

36.

S. P.

Webb

and

S.

Hammes-Schiffer

,

J. Chem. Phys.

113

,

5214

(

2000

).

https://doi.org/10.1063/1.1289528

Crossref

37.

A.

Hazra

,

J. H.

Skone

, and

S.

Hammes-Schiffer

,

J. Chem. Phys.

130

,

054108

(

2009

).

https://doi.org/10.1063/1.3068526

Crossref

PubMed

2018

Author(s)

Stability conditions and local minima in multicomponent Hartree-Fock and density functional theory

I. INTRODUCTION

II. BASIC THEORY

A. NEO-HF stability conditions

B. NEO-DFT stability conditions

III. CONNECTION OF STABILITY ANALYSIS TO NEO-TDHF AND NEO-TDDFT IMAGINARY EIGENVALUES

IV. ELECTRONIC, PROTONIC, AND VIBRONIC INSTABILITIES

V. STABILITY THEORY WITH SPIN AND SPATIAL CONSTRAINTS

A. Complex NEO-GHF

B. Real NEO-GHF

C. Complex NEO-UHF

D. Real NEO-UHF

E. Complex NEO-RHF

F. Real NEO-RHF

G. Complex NEO-MHF

H. Real NEO-MHF

VI. STABILITY ANALYSIS FOR MOLECULAR SYSTEMS

VII. CONCLUSIONS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

REFERENCES

Supplementary Material

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Stability conditions and local minima in multicomponent Hartree-Fock and density functional theory

I. INTRODUCTION

II. BASIC THEORY

A. NEO-HF stability conditions

B. NEO-DFT stability conditions

III. CONNECTION OF STABILITY ANALYSIS TO NEO-TDHF AND NEO-TDDFT IMAGINARY EIGENVALUES

IV. ELECTRONIC, PROTONIC, AND VIBRONIC INSTABILITIES

V. STABILITY THEORY WITH SPIN AND SPATIAL CONSTRAINTS

A. Complex NEO-GHF

B. Real NEO-GHF

C. Complex NEO-UHF

D. Real NEO-UHF

E. Complex NEO-RHF

F. Real NEO-RHF

G. Complex NEO-MHF

H. Real NEO-MHF

VI. STABILITY ANALYSIS FOR MOLECULAR SYSTEMS

VII. CONCLUSIONS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

REFERENCES

Supplementary Material

Citing articles via

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only