Fast Methods For Multisite Charge Transfer Processes II: Analytic Nuclear Gradients and Nonadiabatic Dynamics For cCASSCF(1,M) and cCASSCF(2M-1,M) Wavefunctions

TL;DR

This paper develops and implements analytic nuclear gradients and derivative couplings for a specialized CASSCF method, enabling efficient nonadiabatic dynamics simulations of charge transfer processes.

Contribution

It introduces a new analytic gradient and coupling approach for constrained CASSCF(1,M) and (2M-1,M) wavefunctions, facilitating charge transfer modeling.

Findings

Enabled initial surface-hopping simulations of proton-coupled electron transfer.

Provided a computational framework for smooth potential energy surfaces.

Demonstrated applicability to phenoxyl-phenol system.

Abstract

We derive and implement analytic nuclear gradients and derivative couplings for a constrained Complete Active Space Self-Consistent Field with a small active space designed to model electron or hole transfer. Using a Lagrangian formalism, we are able to differentiate both the CASSCF energy and the constraint (which is required for globally smooth surfaces), and the resulting efficient algorithm can be immediately applied to nonadiabatic dynamics simulations of charge transfer processes. Here, we run initial surface-hopping simulations of a proton coupled electron transfer event for a phenoxyl-phenol system.

Equations

$$g_i(\mathbf{x};\mathbf{R})=0.$$

$$\mathcal{L}(\mathbf{x},\mathbf{R},\mathbf{\lambda })=f(\mathbf{x},\mathbf{R})-{\mathbf{\lambda }}^{\top }\mathbf{g}(\mathbf{x},\mathbf{R}),$$

$$\frac{df}{d\mathbf{R}}={\nabla }_{\mathbf{R}}\mathcal{L}(\mathbf{x},\mathbf{R},\mathbf{\lambda })={\nabla }_{\mathbf{R}}f(\mathbf{x},\mathbf{R})-{\left({\nabla }_{\mathbf{R}}\mathbf{g}(\mathbf{x},\mathbf{R})\right)}^{\top }\mathbf{\lambda }.$$

$${\left({\nabla }_{\mathbf{R}}\mathbf{g}(\mathbf{x},\mathbf{R})\right)}_{ij}={\nabla }_{R_j}{g}_i(\mathbf{x},\mathbf{R}),$$

$$\mathbf{\lambda }={\left({\nabla }_{\mathbf{x}}\mathbf{g}(\mathbf{x},\mathbf{R})\right)}^{-1\top }{\nabla }_{\mathbf{x}}f(\mathbf{x},\mathbf{R}).$$

$$\frac{df}{d\mathbf{R}}={\nabla }_{\mathbf{R}}f+{\left(\frac{d\mathbf{x}}{d\mathbf{R}}\right)}^{\top }{\nabla }_{\mathbf{x}}f.$$

$$\mathbf{g}(\mathbf{x},\mathbf{R})$$

$$\Rightarrow \frac{d\mathbf{g}(\mathbf{x},\mathbf{R})}{d\mathbf{R}}$$

$$\Rightarrow \frac{d\mathbf{x}}{d\mathbf{R}}$$

$${\mathcal{L}}_1(\mathbf{x},\mathbf{R},\mathbf{\lambda })=f(\mathbf{x},\mathbf{R})-{\mathbf{\lambda }}^{\top }\mathbf{g}(\mathbf{x},\mathbf{R}),$$

$$\mathbf{g}(\mathbf{x},\mathbf{R})$$

$${\nabla }_{\mathbf{x}}{\mathcal{L}}_1(\mathbf{x},\mathbf{R},\mathbf{\lambda })$$

$${\mathcal{L}}_2(\mathbf{x},\mathbf{R},\mathbf{\lambda },\mathbf{\mu },\mathbf{\nu })=h(\mathbf{x},\mathbf{R})-{\mathbf{\mu }}^{\top }\mathbf{g}(\mathbf{x},\mathbf{R})-{\mathbf{\nu }}^{\top }{\nabla }_{\mathbf{x}}{\mathcal{L}}_1(\mathbf{x},\mathbf{R},\mathbf{\lambda }),$$

$$\frac{dh}{d\mathbf{R}}={\nabla }_{\mathbf{R}}{\mathcal{L}}_2={\nabla }_{\mathbf{R}}h-({\nabla }_{\mathbf{R}}\mathbf{g}{)}^{\top }\mathbf{\mu }-({\nabla }_{\mathbf{R}}{\nabla }_{\mathbf{x}}{\mathcal{L}}_1{)}^{\top }\mathbf{\nu },$$

$$({\nabla }_{\mathbf{R}}{\nabla }_{\mathbf{x}}{\mathcal{L}}_1{)}_{ij}=({\nabla }_{R_j}{\nabla }_{x_i}{\mathcal{L}}_1).$$

$${\nabla }_{\mathbf{x}}{\mathcal{L}}_2(\mathbf{x},\mathbf{R},\mathbf{\lambda },\mathbf{\mu },\mathbf{\nu })$$

$${\nabla }_{\mathbf{\lambda }}{\mathcal{L}}_2(\mathbf{x},\mathbf{R},\mathbf{\lambda },\mathbf{\mu },\mathbf{\nu })$$

$$\left(\begin{array}{cc}{\nabla }_{\mathbf{xx}}^2{\mathcal{L}}_1& ({\nabla }_{\mathbf{x}}\mathbf{g}{)}^{\top }\\ {\nabla }_{\mathbf{x}}\mathbf{g}& 0\end{array}\right)\left(\begin{array}{c}\mathbf{\nu }\\ \mathbf{\mu }\end{array}\right)=\left(\begin{array}{c}{\nabla }_{\mathbf{x}}h\\ 0\end{array}\right).$$

$$(\text{underaccent}\sim \mathbf{V}{)}_{pq}=\{\begin{array}{ll}{V}_{pq},& p<q,\\ 0,& p=q,\\ -V_{qp},& p>q.\end{array}$$

$${\mathcal{L}}_1=\sum _{i=1}^M{w}_i{E}_i-\sum _{j=1}^{M-1}{\lambda }_j{\mathbf{Q}}^j:{\mathbf{P}}^{\mathrm{active}}$$

$$w_j$$

$$a_j$$

$$\Delta E_j$$

$$E_i=\frac{1}{2}\left(2{\mathbf{H}}^0:{\mathbf{P}}^{i,t}+{\mathbf{P}}^{i,t}:{\mathbf{\Pi }}^J:{\mathbf{P}}^{i,t}-{\mathbf{P}}^{i,\alpha }:{\mathbf{\Pi }}^K:{\mathbf{P}}^{i,\alpha }-{\mathbf{P}}^{i,\beta }:{\mathbf{\Pi }}^K:{\mathbf{P}}^{i,\beta }\right),$$

$${\Pi }_{\rho \sigma \zeta \eta }^J=\int d{\mathbf{r}}^{1}d{\mathbf{r}}^2{\chi }_{\rho }({\mathbf{r}}^1){\chi }_{\sigma }({\mathbf{r}}^1)\frac{1}{|{\mathbf{r}}^1-{\mathbf{r}}^2|}{\chi }_{\zeta }({\mathbf{r}}^2){\chi }_{\eta }({\mathbf{r}}^2),$$

$${\Pi }_{\rho \sigma \zeta \eta }^K=\int d{\mathbf{r}}^{1}d{\mathbf{r}}^2{\chi }_{\rho }({\mathbf{r}}^1){\chi }_{\zeta }({\mathbf{r}}^1)\frac{1}{|{\mathbf{r}}^1-{\mathbf{r}}^2|}{\chi }_{\sigma }({\mathbf{r}}^2){\chi }_{\eta }({\mathbf{r}}^2),$$

$${\mathbf{Q}}^j={\mathbf{\Sigma }}^j-{\mathbf{\Sigma }}^{j+1},$$

$${\mathbf{\Sigma }}^j$$

$$\equiv \mathbf{S}\left(\begin{array}{ccc}0& \cdots & 0\\ ⋮& {\left({\mathbf{S}}^{jj,\mathrm{frag}}\right)}^{-1}& ⋮\\ 0& \cdots & 0\end{array}\right)\mathbf{S}$$

$${\mathbf{P}}^{i,\alpha }$$

$${\mathbf{P}}^{i,\beta }$$

$${\mathbf{P}}^{\mathrm{active}}$$

$$\mathbf{C}={\mathbf{C}}^{\mathrm{old}}{e}^{\mathbf{A}},$$

$$\left(\begin{array}{cc}{\nabla }_{\mathbf{AA}}^2{E}_{\mathrm{tot}}-{\nabla }_{\mathbf{AA}}^2{G}_0& ({\nabla }_{\mathbf{A}}\mathbf{G}{)}^{\top }\\ {\nabla }_{\mathbf{A}}\mathbf{G}& 0\end{array}\right)\left(\begin{array}{c}\mathbf{\nu }\\ \mathbf{\mu }\end{array}\right)=\left(\begin{array}{c}{\nabla }_{\mathbf{A}}{E}_k\\ 0\end{array}\right)$$

$$E_{\mathrm{tot}}$$

Full text

1 Introduction

Charge transfer (CT) processes play an important role in numerous chemical phenomena including redox reactions 1, 2, 3, photochemical processes 4, electrochemistry 5, and catalysis 6, 7. Fundamental theoretical understanding of electronic coupling in these processes has been widely developed 8, 9, 10. Their theoretical investigation is essential for advancing technologies in energy storage, molecular electronics, and designing of novel catalyst, requiring sophisticated theoretical frameworks to properly describe the electronic structure and dynamics involved.11, 12, 13 Despite their importance, accurate and efficient modeling of CT processes remains a challenge due to the complex interplay between electronic and nuclear motion, particularly near conical intersections where nonadiabatic effects become crucial14. On the one hand, single reference, mean-field theories such as Hartree Fock (HF) or density functional theory (DFT) often lack the accuracy to describe CT processes properly, as they fail to capture the multi-reference character of the wavefunction; Time-Dependent DFT (TD-DFT) has been widely applied to study CT processes,15, 10 but the latter treats excited states and the ground state in an unbalanced fashion (and cannot recover the correct topology of a crossing16). On the other hand, multi-reference methods such as Complete Active Space Self-Consistent Field (CASSCF) and Multi-Reference Configuration Interaction (MRCI) provide more accurate descriptions of electron correlation for CT processes.17, 18 However, while recent advances have expanded CASSCF capabilities through density matrix renormalization group (DMRG) theory and provided analytic derivative couplings for selected multi-state CAS implementations, these approaches can be extremely computationally demanding, a fact which prohibits their practical application for most large scale dynamics simulations. This fundamental trade-off between accuracy and efficiency underscores the need for alternative methods that can be both reliable and computational efficient.

In a recent manuscript 19, 20, we have presented a dynamically weighted state-averaged constrained CASSCF(1,M)/CASSCF(2M-1,M) approach, or eDSC/hDSC for short, to accurately generate PES and wavefunctions for studying CT processes in radical systems. For a two-state problem, our current implementation requires only twice the cost of HF calculations, and for a $M$-state problem, the algorithm requires at most $M$ times the cost of HF calculations. Given its low cost and balanced good accuracy, the development of analytic nuclear gradients and derivative couplings for this method has now become essential to enable its immediate application in nonadiabatic dynamics simulations. Gradients and derivative couplings are needed also for efficient transition state searches as well as to accurately capture electronic transitions at avoided crossings, which is critical for studying CT mechanisms with proper treatment of quantum effects.

With this background in mind, in this paper, we present the required mathematical formalism for computing such analytic nuclear gradients and derivative couplings through a Lagrange multiplier approach.21, 22, 23 We have implemented the formalisms in QChem24 and verified these expressions by comparing them with finite-difference results. As a proof-of-concept, we have used our analytic nuclear gradients and derivative couplings to run surface hopping (SH) simulations of a proton-coupled electron transfer (PCET) process, demonstrating the practical utility of our approach. With the current implementation, computing the analytic nuclear gradient and derivative coupling requires approximately the same computational cost as converging the electronic structure calculation within eDSC/hDSC for a system with 25 atoms and 232 basis functions (the phenoxyl-phenol system with 6-31g* basis set), which is only twice the cost of HF calculations. Thus, this work demonstrates the immediate feasibility of nonadiabatic dynamics simulation for studying CT processes of large systems, a capability which should open the door to numerous research fields including photocatalysis, electrochemistry, transition metal catalysis, proton-coupled electron transfer, and the design of molecular electronics.

An outline of this paper is as follows. In Sec. 2, we briefly review (in general) the Lagrange multiplier approach for computing the analytic nuclear gradient. In Sec. 3, we demonstrate how to apply such a theory to our constrained optimization problem. In Sec. 4, we construct the relevant equations for the $M$-state eDSC/hDSC problem. Specifically, in Secs. 5-7, we derive the necessary terms to compute the analytic nuclear gradient of the relevant adiabatic energies and in Sec. 8, we derive the necessary terms to for the derivative coupling. In Sec. 9, we apply our formalism to SH simulations of a proton-coupled electron transfer process. Finally, in Sec. 10 we summarize and point out new directions of this work.

2 Background: Review of Analytic Nuclear Gradient through Lagrange Multiplier Approach

Let $\{x_{i}\}$ be a set of independent variables and $\{R_{i}\}$ be a set of parameters. Let $g_{i}(\bm{x};\bm{R})$ be a set of constraints on $\bm{x}$ that depends parametrically on $\bm{R}$, where $\bm{x}$ is a shorthand for $(x_{1},x_{2},\cdots)^{\top}$ and $\bm{R}$ is shorthand for $(R_{1},R_{2},\cdots)^{\top}$. For clarity, one can think of $\bm{x}$ as some parameterization of the electronic wavefunction (e.g., the molecular orbital coefficients in the atomic orbital basis) while $\bm{R}$ represents the nuclear coordinates in a Hartree-Fock framework. A Lagrange multiplier approach provides a method to compute the total derivative of any target function w.r.t. $\bm{R}$ when working with fully determined systems 25, 26. To that end, suppose we have enough constraints $g_{i}(\bm{x};\bm{R})$ such that $\bm{x}$ can be determined locally by equations

[TABLE]

Then, for any arbitrary target function $f(\bm{x},\bm{R})$, one can build a Lagrangian

[TABLE]

where $\bm{\lambda}=(\lambda_{1},\lambda_{2},\cdots)^{\top}$ is a collection of Lagrange multipliers and $\bm{g}(\bm{x},\bm{R})=[g_{1}(\bm{x},\bm{R}),g_{2}(\bm{x},\bm{R}),\cdots]^{\top}$. According to a Lagrange multiplier approach, we can compute the total derivative of $f$ w.r.t. $\bm{R}$ on the constrained manifold (solution manifold of Eq. 1) by evaluating

[TABLE]

Some words about notation are now crucial. Throughout the paper, we will use $\nabla_{\bm{}}$ to represent the gradient taken on the free manifold (i.e., all quantities $\bm{x}$, $\bm{\lambda}$, and $\bm{R}$ are independent of each other) and $\frac{d{}}{d{\bm{R}}}$ to represent the gradient taken on the constrained manifold (i.e., where $\bm{x}$ and $\bm{\lambda}$ depend on $\bm{R}$). With this convention, $\nabla_{\bm{R}}$ in Eq. 3 is taken with $\bm{x}$ and $\bm{\lambda}$ fixed, $\nabla_{\bm{R}}\bm{g}(\bm{x},\bm{R})$ is a matrix indexed by

[TABLE]

and $\bm{\lambda}$ is the solution to the equations $\nabla_{\bm{x}}\mathcal{L}(\bm{x},\bm{R},\bm{\lambda})=\bm{0}$, i.e.,

[TABLE]

Note that since $\bm{x}$ can be fully determined by equations $\bm{g}(\bm{x},\bm{R})=\bm{0}$, $\nabla_{\bm{x}}\bm{g}(\bm{x},\bm{R})$ is a full-rank square matrix. With this approach, the total derivative of any arbitrary function $f(\bm{x},\bm{R})$ w.r.t. $\bm{R}$ can be computed (through Eqs. 3 and 5) without computing $\frac{d{\bm{x}}}{d{\bm{R}}}$. Note that there is no mystery as far as deriving Eq. 3. The simplest means to derive this result is to write

[TABLE]

And since $\frac{d{}}{d{\bm{R}}}$ is evaluated on the constrained manifold, one has

[TABLE]

Eq. 3 follows from Eqs. 6 and 9 provided that $\nabla_{\bm{x}}\bm{g}(\bm{x},\bm{R})$ is invertible, i.e., the constraints determine the system.

3 Theory: Analytic Nuclear Gradient for Constrained Optimization Problem

The theory above is not general enough to handle most problems in quantum chemistry. For electronic structure problems, it is often the case that, in order to find variables $\bm{x}$, we must first optimize a function $f(\bm{x},\bm{R})$ satisfying a constraint $g(\bm{x},\bm{R})=0$; thereafter, second, we need to differentiate a function $h(\bm{x},\bm{R})$. In other words, rather than investigating a fully determined by a constraint, instead we must first solve a constrained optimization problem and then further differentiate. This case is clearly more complicated than the optimization treated above.

To study such a case, let the Lagrangian for the initial constrained optimization problem be

[TABLE]

Note that the solution for $\bm{x}$ can be completely determined by the constraints that

[TABLE]

Therefore, for a target function $h(\bm{x},\bm{R})$ that we wish to differentiate, we can build a second Lagrangian

[TABLE]

such that $\bm{x}$ is now fully determined by the constraints that are contracted with the Lagrange multipliers $\bm{\mu}$ and $\bm{\nu}$. Following the approach discussed in Sec. 2, the gradient of $h(\bm{x},\bm{R})$ can be evaluated by

[TABLE]

where $(\nabla_{\bm{R}}\nabla_{\bm{x}}\mathcal{L}_{1})$ is a matrix indexed by

[TABLE]

As always, here a partial derivative $\nabla_{\bm{}}$ is taken by fixing the remaining four quantities $\{\bm{x},\bm{R},\bm{\lambda},\bm{\mu},\bm{\nu}\}$.

Now, all that remains is to determine Lagrange multipliers $\bm{\mu}$ and $\bm{\nu}$. The crucial observation here is that, if we follow Sec. 2 above, there are actually two sets of equations for determining these Lagrange multipliers, namely:

[TABLE]

While Eq. 16 is naturally derived given that $\bm{x}$ is our independent variable in Eqs. 11-12, the key insight here is to notice that Eqs. 11-12 determine not only $\bm{x}$ but also $\bm{\lambda}$; to some extent, $\bm{\lambda}$ plays the same role as $\bm{x}$ in a constrained optimization problem which leads to Eq. 17. In any event, combining Eqs. 16 and 17, the equations for $\bm{\mu}$ and $\bm{\nu}$ are

[TABLE]

Interestingly, Eq. 18 is very similar to the Newton-KKT equations 27. Once $\bm{\mu}$ and $\bm{\nu}$ are solved from this equation, the gradient of the target function $h(\bm{x},\bm{R})$ can be evaluated by Eq. 14 without computing $\frac{d{\bm{x}}}{d{\bm{R}}}$ explicitly.

4 Brief Review of eDSC/hDSC and Notation Used Below

Let us now show how to apply the standard Lagrange multiplier approach above to compute the analytic nuclear gradient of the adiabatic energies and derivative couplings for constrained dynamically weighted state-averaged self-consistent field calculations with either $(i)$ one electron in $M$ orbitals (corresponding to “electron transfer”) or $2M-1$ electrons in $2M$ orbitals (corresponding to “hole transfer”). While one might refer to these methods as cDW-SA-CASSCF$(1,M)$ or cDW-SA-CASSCF$(2M,2M-1)$ calculations, henceforward we will use the even simpler acronyms eDSC and hDSC (for electron/hole dynamically weighted, state averaged, CASSCF) 19, 20. Note that eDSC/hDSC use restricted orbitals, so that the spin up and spin down orbitals are equivalent spatially. Furthermore, to avoid the complexity of notation from indices of vector/matrix/tensor elements, here we introduce several shortcuts for notation used below.

• •

Einstein summation convention is applied when the same indices corresponding to vector/matrix/tensor elements appear twice.

• •

A contraction notation (${:}$) is introduced such that for any matrix $\bm{A}$ and $\bm{B}$, we define $\bm{A}{:}\bm{B}=A_{ij}B_{ij}={\rm Tr}(\bm{A}^{\top}\bm{B})$.

• •

The contraction notation is also applied between any matrix $\bm{A}$ and rank-4 tensor $\bm{\Pi}$ such that $\bm{\Pi}{:}\bm{A}$ is a matrix indexed by $(\bm{\Pi}{:}\bm{A})_{ij}=\Pi_{ijkl}A_{kl}$, and $\bm{A}{:}\bm{\Pi}$ is a matrix indexed by $(\bm{A}{:}\bm{\Pi})_{kl}=A_{ij}\Pi_{ijkl}$.

• •

The contraction between any two rank-4 tensors $\bm{\Pi}$ and $\bm{\Theta}$ is also a rank-4 tensor, which is indexed by $(\bm{\Pi}{:}\bm{\Theta})_{ijkl}=\Pi_{ijmn}\Theta_{mnkl}$.

• •

For given vector $\bm{V}$ indexed by $pq$ (with elements $V_{pq}$ where $p<q$), it will be helpful to construct the corresponding antisymmetric matrix. To that end, we define $\underaccent{\ \sim}{\bm{V}}$ as the matrix with elements:

[TABLE]

As far as matrices and basis sets for quantum chemistry are concerned, the one-electron Hamiltonian $\bm{H}^{0}$, the two-electron interaction tensors $\bm{\Pi}^{J}$ (Coulomb) and $\bm{\Pi}^{K}$ (exchange), and all density matrices $\bm{P}$s are expressed in an atomic orbital (AO) basis, $\{\ket{\chi_{\rho}}\}$. All quantities that are explicitly contracted with these quantities are also expressed in an AO basis. There will be occasions below where matrices are expressed in a molecular orbital (MO) basis; for the most part, such matrices are written with a tilde, $\tilde{\bm{o}}$. Finally, we will use a bar, $\bar{\bm{o}}$, to indicate a matrix expressed in an orthonormalized AO basis (and we will of course define the basis transformation explicitly when introducing such matrices). There is no distinction between upper and lower case variables. Lastly, as far as matrix/tensor indices are concerned, we use $\rho,\sigma,\zeta,\eta$ for indices in an AO basis, $p,q,r,s$ for indices in a MO basis, and $i,j,k,l,m,n$ for general indices.

At this point, let us briefly review the eDSC/hDSC formalism. The readers may refer to Refs. 19 and 20 for more details. For an $M$ state problem, the Lagrangian is

[TABLE]

where $w_{i}$ is the dynamical weighting factor, $E_{i}$ is the energy of the $i{\rm th}$ electron configuration, $\lambda_{j}$ is the Lagrange multiplier for the $j{\rm th}$ constraint, $\bm{Q}^{j}$ is the kernel of the $j{\rm th}$ constraint ($\bm{Q}_{j}$ is a matrix that depends only on the nuclear geometry $\bm{R}$), and $\bm{P}^{\rm active}$ is the projector onto the vector space spanned by all active orbitals (or equivalently, the density matrix from all active orbitals). For concreteness, here we will review each term in detail.

• •

The dynamical weighting factors are defined by

[TABLE]

where $T$ is a fixed “temperature” parameter and configuration 1 is assumed to have the lowest energy among all electron configurations.

• •

The energy for state $i$ can be written as

[TABLE]

where $\bm{H}^{0}$ is the one-electron Hamiltonian, $\bm{\Pi}^{J}$ is the Coulomb interaction tensor,

[TABLE]

$\bm{\Pi}^{K}$ is the exchange interaction tensor,

[TABLE]

• •

$\bm{P}^{i,\alpha}$ is the spin-up density matrix of the $i{\rm th}$ electron configuration, $\bm{P}^{i,\beta}$ is the spin-down density matrix of the $i{\rm th}$ electron configuration, and $\bm{P}^{i,t}$ is the total density matrix summed over both spins, i.e., $\bm{P}^{i,t}=\bm{P}^{i,\alpha}+\bm{P}^{i,\beta}$.

The idea behinds eDSC/hDSC is to constrain the vector space spanned by all active orbitals projects equally onto each atomic fragment (charge center) so as to treat electron and hole transfer. To achieve this feat, we write

[TABLE]

where $\bm{\Sigma}^{j}$ defines the projection onto the vector space spanned by the $j$th atomic fragment. In the atomic orbital (AO) basis, we find that it is convenient to express

[TABLE]

where $\bm{S}$ is the overlap matrix between AOs and $\bm{S}^{jj,{\rm frag}}$ is the diagonal block of $\bm{S}$ that belongs to the $j$th atomic fragment (charge center). Note that if the $j$th fragment includes all atoms, $\bm{\Sigma}^{j}=\bm{S}$ and ${\rm Tr}(\bm{\Sigma}^{j}\bm{P}^{\rm active})=M$.

Now, when choosing the property independent variable for optimization, note that the density matrices are suboptimal because of complications from the orthonormality constraints. Instead, a cleaner strategy is to parameterize the Lagrangian in Eq. 20 by the generator of the rotation matrix for the molecular orbitals (MO). Thus, we write

[TABLE]

and

[TABLE]

Here, the $\bm{N}$ matrices are occupation matrices, which are diagonal with unity diagonal elements if the orbital is occupied and zero otherwise. $\bm{C}$ is the MO coefficient and $\bm{A}$ is an anti-symmetric generator matrix. With this parameterization, the $\nabla_{\bm{x}}$ operator in Eq. 14 can be understood as $\nabla_{\bm{A}}$.

Let us now derive the working equations for the analytic nuclear gradient of adiabatic energies (Secs. 5-7) and the derivative couplings (Sec. 8). Obviously, the key obstacle is to extract the necessary lagrange multipliers, which we do now.

5 Extracting Lagrange Multipliers (Using Eq. 18)

In the context of eDSC/hDSC calculations, if we seek the analytic nuclear gradient of the $k$th adiabatic energy ($E_{k}$), Eq. 18 reads as follows

[TABLE]

where

[TABLE]

Let us begin by evaluating the terms on the L.H.S of Eq. 34.

5.1 Term $\nabla_{\bm{AA}}^{2}E_{\rm tot}$

We show in Appendix 12.A that, given the definition of $w_{i}$ in Eqs. 21-23, one can write

[TABLE]

where

[TABLE]

and $E_{i}$ is defined in Eq. 24 above. Again, we assume that $E_{1}$ is the lowest state energy. According to Eq. 39, we can evaluate $\nabla_{\bm{AA}}^{2}E_{\rm tot}$ once we compute $\nabla_{\bm{AA}}^{2}E_{i}$ and $\nabla_{\bm{A}}E_{i}$. In order to calculate the latter two terms, recalling Eq. 24, let us use the notation

[TABLE]

such that

[TABLE]

It is then straightforward to differentiate and we find (for $q>p$, $s>r$)

[TABLE]

where $[\bm{A},\bm{B}]=\bm{AB}-\bm{BA}$ is the commutator between two matrices $\bm{A}$ and $\bm{B}$, and

[TABLE]

In order to further simplify the notation that will be used later, one may also define a rank-4 tensor $\bm{\Theta}^{1}$

[TABLE]

(using Eq. 52 for $\tilde{\bm{W}}$) such that

[TABLE]

Similarly, for $\epsilon_{2}(\bm{\Pi},\bm{P})$, one finds

[TABLE]

where no summation is applied on $p,q,r,s$ and we have defined:

[TABLE]

Lastly, one can define another rank-4 tensor $\bm{\Theta}^{2}$

[TABLE]

which simplifies Eq. 56 to

[TABLE]

At this point, we can evaluate $\nabla_{\bm{A}}E_{i}$ and $\nabla_{\bm{AA}}^{2}E_{i}$ for each adiabatic state. To make our final state-averaged expression most easy to express, let us further define several analogues of Eq. 57-58 above for each state $i$ (where we contract the two electron tensors and each state’s density matrix, either in the AO or MO basis):

[TABLE]

and

[TABLE]

The spin-up (spin-down) Fock matrices for the $i$th electron configuration in the MO basis are naturally defined as follows:

[TABLE]

5.1.1 Final Form for $\frac{\partial{E_{i}}}{\partial{A_{pq}}}$

With the notation and derivation above, the first derivative of the energy for the $i$th electronic configuration is

[TABLE]

5.1.2 Final Form for $\frac{\partial{{}^{2}E_{i}}}{\partial{A_{rs}A_{pq}}}$

Finally, given Eqs. 54 and 61, we can write

[TABLE]

5.1.3 A Scheme for Iterative Inversion of Eq. 34

In practice, one cannot afford to invert the L.H.S. of Eq. 34. Rather, one must compute its contraction with a vector and solve iteratively. To that end, let us now derive an expression for the contraction of Eq. 71 with an arbitrary vector $\bm{V}$ (with antisymmetric analogue $\underaccent{\ \sim}{\bm{V}}$ as in Eq. 19).

To begin with, note that, for an arbitary vector/matrix $\bm{V}$/$\underaccent{\ \sim}{\bm{V}}$, one can easily verify that:

[TABLE]

Here, $\tilde{\bm{H}}$ and $\bm{N}$ are obviously symmetric matrices, and

[TABLE]

is also a symmetric matrix. One can further verify that

[TABLE]

where

[TABLE]

Given Eqs. 72-76, the path to computing $\left(\partial^{2}E_{i}/\partial A_{rs}\partial A_{pq}\right)V_{rs}$ is clear. Namely, one calculates the following quantities in sequential order:

[TABLE]

so that, in the end, we can write

[TABLE]

and

[TABLE]

5.1.4 Final Form for a Contraction of $\frac{\partial{{}^{2}E_{\rm tot}}}{\partial{A_{rs}A_{pq}}}$ with $\bm{V}$

If we define

[TABLE]

recall Eq. 39, we can finally write

[TABLE]

5.2 Term $\nabla_{\bm{AA}}^{2}G_{0}$

The $G_{0}$ term can be treated in the same way as the one-electron energy term in Eq. 46. If we define

[TABLE]

we can write

[TABLE]

Now, by using Eqs. 54 and 72-73, we have

[TABLE]

where

[TABLE]

5.3 Term $\nabla_{\bm{A}}\bm{G}$

In Eq. 34, we must make a distinction between the term $\nabla_{\bm{A}}\bm{G}$ and the termn $(\nabla_{\bm{A}}\bm{G})^{\top}$, as the two quantities contract with vectors in different ways. The term $(\nabla_{\bm{A}}\bm{G})^{\top}$ will be discussed in the next section. For a vector $\bm{V}$ that has the same dimension as in Eqs. 88 and 91, we can write

[TABLE]

where $G_{j}$ is defined in Eq. 37 and one has

[TABLE]

Each row on the R.H.S. of Eq. 94 can be written as

[TABLE]

5.4 Term $(\nabla_{\bm{A}}\bm{G})^{\top}$

For a vector $\bm{\mu}$ that has the length of $M-1$, the contraction between $(\nabla_{\bm{A}}\bm{G})^{\top}$ and $\bm{\mu}$ is given by

[TABLE]

This completes our analysis of 34; we have now a practical way to solve the necessary Lagrange multipliers $\bm{\mu}$ and $\bm{\nu}$.

6 Explicit Gradient w.r.t. Nuclear Geometry (Using Eq. 14

At this point, we are in a position to evaluate Eq. 14 for the nuclear gradient of the energy state $k$. Let us write out Eq. 14 as:

[TABLE]

Given $\bm{\nu}$ and $\bm{\mu}$, all that remains is to evaluate the four explicit gradient terms, $\nabla_{\bm{R}}E_{k}$, $\nabla_{\bm{R}}\bm{G}$, $\nabla_{\bm{R}}\nabla_{\bm{A}}E_{\rm tot}$, and $\nabla_{\bm{R}}\nabla_{\bm{A}}G_{0}$.

Now, when evaluating the explicit gradient terms above, e.g. $\nabla_{\bm{R}}E_{k}$, care must be still be taken when work in an atomic orbital basis – even though we have already derived the necessary Lagrange multipliers. To understand the subtlety, note the target function is of the form $E_{k}(\bm{R},\bm{C}(\bm{R}))$ (where $\bm{C}$ is the matrix of MOs in an atomic orbital basis), but we have chosen above (in Sec. 4) to substitute $E_{k}(\bm{R},\bm{C}(\bm{R}))$ with $E_{k}(\bm{R},\bm{A}(\bm{R}))$ using Eq. 33. One must ask: can we safely define (in the sense of Lagrange multiplier approach) as in Eq. 100?

[TABLE]

Unfortunately, this definition cannot be correct because $\bm{C}$ at different geometry $\bm{R}$ satisfies different orthonormality constraints:

[TABLE]

where $\bm{R}^{1}$ and $\bm{R}^{2}$ are two nuclear geometries, $\bm{S}(\bm{R}^{1})$ and $\bm{S}(\bm{R}^{2})$ are the AO overlap matrices at these two geometries, and $\bm{C}^{1}$ and $\bm{C}^{2}$ are MO coefficients in AO basis at these two geometries. In other words, in practice, we should actually consider the target function to be $E_{k}(\bm{R},\bm{C}(\bm{R}),\bm{S}(\bm{R}))$.

The distinction above indicates that, for a proper gradient, one must take into acccount the gradient of the nuclear overlap matrix, which leads to well-known Pulay forces in gradient theory. For our purposes, we can derive the relevant corrections simply by orthonormalize the AO basis and write $\bm{C}$ in the orthonormal AO basis as the independent variable. Specifically, one can define

[TABLE]

where $\bm{T}=\bm{S}^{1/2}$ is the square root of the AO overlap matrix $\bm{S}$, and treat $\bar{\bm{C}}$ as our independent variable (still parameterized by $\bm{A}$). The advantage of this choice is that $\bar{\bm{C}}$ is orthonormalized in any geometry, i.e., $\bar{\bm{C}}^{\top}\bar{\bm{C}}=\bm{I}$, and there is no need to introduce an extra constraint term for the orthonormality of orbitals. Moreover, since $\nabla_{\bm{A}}$ evaluates the response of MOs that does not depend on the basis in which MOs are expressed, changing $\bm{C}$ to $\bar{\bm{C}}$ does not influence any results derived above (one can verify that the AO basis or orthonormal AO basis are both contracted out in expressions for the $\nabla_{\bm{A}}$ terms). Given this clarification of independent variables, it is also clear what quantities are considered “explicitly geometry dependent” in Eq. 99. One should express $E_{k},\bm{G},E_{\rm tot}$, and $G_{0}$ in terms of $\bar{\bm{C}}$, and all quantities except $\bar{\bm{C}}$ (and the Lagrange multipliers $\bm{\lambda}$) are considered geometry dependent. We demonstrate this treatment using $E_{k}$ as an example below, the same approach applies to all terms in Eq. 99.

For $E_{k}$, we first write

[TABLE]

where

[TABLE]

given that $\bm{X}=\bm{S}^{-1/2}$ and $\bm{X}{\otimes}\bm{X}$ is a rank-4 tensor indexed by

[TABLE]

Now, the explicit geometry gradient should be applied to $\bar{\bm{H}}^{0}$, $\bar{\bm{\Pi}}^{J}$, and $\bar{\bm{\Pi}}^{K}$. For instance,

[TABLE]

where we have introduced a superscript ${[\bm{R}]}$ indicating that the geometry gradient is applied to each element of the corresponding matrix. Note that the $\bm{H}^{0{[\bm{R}]}}$ term itself involves two terms, that the change of Coulomb potential (and other possible one-electron potentials) and the change of the matrix element due to the change of a geometry dependent basis. With this setup, we will evaluate each term in Eq. 99.

6.1 Term $\nabla_{\bm{R}}E_{k}$

One can verify that the explicit nuclear gradient for $E_{k}$ is

[TABLE]

where

[TABLE]

and $\bm{T^{[\bm{R}]}}$ can be evaluated easily from $\bm{S}^{[\bm{R}]}$, see Appendix 12.B.

6.1.1 Term $(\nabla_{\bm{R}}\bm{G})^{\top}\bm{\mu}$

In a similar fashion, we can write

[TABLE]

6.1.2 Term $(\nabla_{\bm{R}}\nabla_{\bm{A}}G_{0})^{\top}\bm{\nu}$

The $(\nabla_{\bm{R}}\nabla_{\bm{A}}G_{0})^{\top}\bm{\nu}$ term and the term $(\nabla_{\bm{R}}\nabla_{\bm{A}}E_{\rm tot})^{\top}\bm{\nu}$ in Eq. 99 are slightly different from the former two terms but can be treated in the similar way if the contraction with the Lagrange multiplier $\bm{\nu}$ is performed first. By introducing a matrix $\underaccent{\ \sim}{\bm{\nu}}$ from $\bm{\nu}$ in the same way as expressed in Eq. 19, we can write

[TABLE]

where

[TABLE]

for any two matrices $\bm{A}$ and $\bm{B}$, and

[TABLE]

6.1.3 Term $(\nabla_{\bm{R}}\nabla_{\bm{A}}E_{\rm tot})^{\top}\bm{\nu}$

Finally,

[TABLE]

where $\tilde{\bm{D}}^{j}$ is defined in Eqs. 77-79 by replacing $\underaccent{\ \sim}{\bm{V}}$ with $\underaccent{\ \sim}{\bm{\nu}}$ and

[TABLE]

where $\tilde{\bm{Z}}^{i}$ is defined in Eqs. 80-84 by replacing $\underaccent{\ \sim}{\bm{V}}$ with $\underaccent{\ \sim}{\bm{\nu}}$.

7 Summary of Working Equations for the Nuclear Gradient

Here we summarize the working equations for the analytic nuclear gradient of the adiabatic energy $E_{k}$. The expression is

[TABLE]

where $\bm{\mu}$ and $\bm{\nu}$ are solved by Eq. 34 through iterative methods, given that $\nabla_{\bm{A}}E_{k}$ is evaluated by Eq. 70; for any vector $\bm{V}$ of the same size as $\bm{\nu}$ (note that $\bm{V}$ is usually indexed by two indices $V_{pq},p<q$), $\left(\nabla_{\bm{AA}}^{2}E_{\rm tot}\right){:}\bm{V}$ is evaluated by Eq. 88, $\left(\nabla_{\bm{AA}}^{2}G_{0}\right){:}\bm{V}$ is evaluated by Eq. 91, $(\nabla_{\bm{A}}G_{j>0}){:}\bm{V}$ is evaluated by Eq. 97; and for any vector $\bm{\mu}$ of length $M-1$, $(\nabla_{\bm{A}}\bm{G})^{\top}\bm{\mu}$ is evaluated by Eq. 98.

After solving for $\bm{\mu}$ and $\bm{\nu}$, one can evaluate Eq. 123 term by term, where $\nabla_{\bm{R}}E_{k}$ is evaluated by Eq. 114, $(\nabla_{\bm{R}}\bm{G})^{\top}\bm{\mu}$ is evaluated by Eq. 117, $(\nabla_{\bm{R}}\nabla_{\bm{A}}G_{0})^{\top}\bm{\nu}$ is evaluated by Eq. 118, and finally, $\left(\nabla_{\bm{R}}\nabla_{\bm{A}}E_{\rm tot}\right)^{\top}\bm{\nu}$ is evaluated by Eqs. 121 and 122.

8 Derivative Coupling

Having derive the relevant equations for the energy gradient, let us turn our attention to the the derivative coupling between two adiabatic states. A result can be found with only minor revisions on the procedure proposed in the previous section. Let us define

[TABLE]

as the derivative coupling between adiabatic states $\left|\Psi_{k}\right>$ and $\left|\Psi_{l}\right>$, where both $\left|\Psi_{k}\right>$ and $\left|\Psi_{l}\right>$ can be expressed as a single Slater determinant of all the occupied orbitals. Since different adiabatic states differ only by their active orbitals, Eq. 124 is equivalent to

[TABLE]

where $\{\left|\phi\right>\}$ are one-electron orbitals and $a_{l}$ is the index of the $l{\rm th}$ active orbital.

Now, in an AO basis $\{\left|\chi_{\rho}\right>\}$, any one-electron orbital can be expressed as

[TABLE]

where $C_{\rho i}$ defines a matrix of molecular orbitals, $\bm{C}$. Eq. 125 can then be expanded as

[TABLE]

The first term on the R.H.S. of Eq. 128 is straightforward to compute; our goal is to compute the second term, i.e., $\sum_{\rho\sigma}C_{\rho a_{k}}S_{\rho\sigma}\frac{d{}}{d{\bm{R}}}C_{\sigma a_{l}}$.

Next, consider an analogous target function defined in a basis of orthonormal atomic orbitals,

[TABLE]

where both $\bar{\bm{C}}^{0}$ and $\bar{\bm{C}}$ are defined in Eq. 103 – except that henceforward, $\bar{\bm{C}}^{0}$ is considered a constant matrix when taking any derivatives. One can easily show that

[TABLE]

At this point, we notice that the first term of Eq. 131 is again trivial to compute, and the second term is identical to $\sum_{\rho\sigma}C_{\rho a_{k}}S_{\rho\sigma}\frac{d{}}{d{\bm{R}}}C_{\sigma a_{l}}$ found above. Therefore, we may conclude that:

[TABLE]

Eq. 132 offers us an approach to compute the derivative coupling using the formalism derived above in Sec. 7. After all, the Lagrangian approach is general for computing the analytic nuclear gradient of any target function so that, in principle, computing $\frac{d{D^{kl}}}{d{\bm{R}}}$ should be no more difficult than computing $\frac{d{E_{k}}}{d{\bm{R}}}$. In fact, only two changes are required.

(1) For the first term on the right hand side of Eq. 123, one must replace $\nabla_{\bm{R}}E_{k}$ with

[TABLE]

The vanishing of the matrix element above arises from the argument that $\bar{\bm{C}}$ is the independent variable and $\bar{\bm{C}}^{0}$ is a constant matrix. Therefore, we can rewrite Eq. 123 explicitly for the purpose of computing derivative couplings as:

[TABLE]

(2) In Eq. 34, one must replace $\nabla_{\bm{A}}E_{k}$ with $\nabla_{\bm{A}}D^{kl}$ where

[TABLE]

After solving for the Lagrange multipliers $\bm{\mu},\bm{\nu}$ from the revised Eq. 34, $\frac{d{D^{kl}}}{d{\bm{R}}}$ is computed through Eq. 134, and finally the derivative coupling $\bm{d}^{kl}$ is computed through Eq. 132.

9 Numerical Results

The derivations in Secs. 5 and 8 have been verified by comparing our computations with finite different results. As a proof-of-concept, we will now demonstrate the practical utility of our approach by using analytic nuclear gradients and derivative couplings to study a proton-coupled electron transfer (PCET) process in a phenoxy-phenol system through a surface hopping (SH) simulation. 28 Our target is the phenoxy-phenol molecule, which as been well studied by the Hammes-Schiffer group29. To begin our analysis, using the analytic hDSC gradient, we optimized the nuclear geometry of the complex and relaxed the system to its ground-state minimum.19 The geometry is shown in Fig. 1A. Next, we excited the system to the first electronic excited state and ran SH dynamics as initialized with atomic velocities from a Maxwell-Boltzmann distribution corresponding to 300 K. The dynamics employed 100 independent trajectories, each propagated with a classical time step of 0.242 fs for a total simulation time of 60.5 fs. A two-quantum time step 30 scheme was applied to interpolate the Hamiltonian and the wavefunction propagation.

In Fig. 2, we plot the excited-state population $(A)$ and the ensemble average of three different bond lengths $(B)$. From our simulation, we find that proton transfer occurs in approximately 10 fs (see the red and blue curves in Fig. 2B). After the transfer, the bridging hydrogen atom never transfers back, but dissipates the energy into the phenoxy ring and oscillates around the oxygen atom. A schematic plot of the hydrogenic motion is shown in Fig. 1B. Note that all trajectory hops from the excited state to the ground state occur in a short time duration that aligns with the proton transfer, in agreement with the notion that the proton transfers through a diabatic process.

As far as the efficiency of our implementation is concerned, each SCF calculation takes 25 seconds while each gradient calculation requires 34 seconds on a 48-core node. 14 seconds out of the 34 seconds needed for gradient calculations are spent on preparing basis function integrals and tensor contractions engine (such as an engine to calculate the contraction between $\bm{\Pi}^{[\bm{R}]}$ and the density matrix), which are prime candidates for further optimization in the near future. Even with these inefficiencies, however, we find the complete surface hopping simulation for the phenoxy-phenol system requires approximately 1 minute per classical time step and 4 hours per trajectory. Thus, the present simulations demonstrates that our method provides a computationally feasible approach for studying CT processes and nonadiabatic dynamics in systems of practical interest going forward.

10 Discussion and Conclusions

In this paper, we have derived and implemented the analytic nuclear gradient and derivative coupling for the recently developed $(1,M)$ and $(2M-1,M)$ dynamically-weighted constrained CASSCF method (so-called eDSC/hDSC) through a Lagrangian formalism. Our implementation provides an efficient way to study CT processes with globally smooth surfaces, making nonadiabatic dynamics simulations of large molecular systems feasible for short and long times. This work can be immediately applied to various research fields, including the study of excited-state dynamics in photocatalytic systems with explicit treatment of quantum effects at avoided crossings, modeling electron transfer reactions at electrode interfaces, describing the multi-reference character of transition metal complexes during catalytic cycles, simulating the charge transport mechanism in molecular junctions, as well as many other physical processes.

Future directions for this work include further optimization of the integral computation routines, which currently account for a significant portion of the gradient calculation time. We are also extending eDSC/hDSC to incorporate the generalized Hartree-Fock (GHF) formalism with spin-orbit coupling (SOC), enabling the study of spin dynamics in radical systems where both charge transfer and spin evolution play critical roles. There is today a great deal of research pointing out that electron transfer may well be spin-polarized – even in organic systems without much spin polarization.31 To that end, a GHF framework for eDSC/hDSC calculations will serve as a zeroth order algorithm for implementing a phase-space Hamiltonian approach that couples nuclear motion with electronic and spin degrees of freedom and propagates electron transfer and nonadiabatic dynamics in a manner that conserves the total linear and angular momentum32, 33. Thus, at the end of day, our hope is that the present algorithm should be an initial step in a exciting search for efficient and accurate large-scale nonadiabatic simulations of electron transfer in complex molecular systems.

11 Acknowledgments

12 Appendix

12.A Second Differential of $E_{\rm tot}$

In this section we derive the second differential of $E_{\rm tot}$ given the dynamical weighting functions in Eqs. 21-23. We start with $a_{j>1}=\frac{1-e^{-\beta\Delta E_{j}}}{\beta\Delta E_{j}}$ (note that $a_{1}=1$), which gives

[TABLE]

We can define

[TABLE]

such that

[TABLE]

Next,

[TABLE]

where we have defined

[TABLE]

Now, instead of working with $E_{\rm tot}$, it is beneficial to define

[TABLE]

such that one has

[TABLE]

where we have defined

[TABLE]

Finally,

[TABLE]

If we define

[TABLE]

then we can write

[TABLE]

12.B Nuclear Gradient of $\bm{T}=\bm{S}^{1/2}$

To compute $\bm{T}^{[\bm{R}]}$ from $\bm{S}^{[\bm{R}]}$, note that

[TABLE]

Now, let us suppose we diagonalize $\bm{T}$:

[TABLE]

Therefore, using the above two equations, it follows that:

[TABLE]

where

[TABLE]

Since $\bm{\Lambda}$ is diagonal, the solution to $\bm{V}$ is straightforward,

[TABLE]

and finally,

[TABLE]