Fast Methods For Multisite Charge Transfer Processes I: Constrained, State Averaged CASSCF(1,M) and CASSCF(2M-1,M) Simulations

TL;DR

This paper introduces a new constrained, state-averaged CASSCF method tailored for simulating charge transfer across multiple molecular fragments, enabling efficient nonadiabatic dynamics simulations.

Contribution

It develops a dynamically-weighted, constrained CASSCF approach for multi-fragment charge transfer, improving computational efficiency and applicability to complex systems.

Findings

Successfully modeled exponential decay of diabatic couplings in SSH chain

Efficiently solved constrained optimization with DIIS-SQP algorithm

Enabled nonadiabatic dynamics simulations for multi-state charge transfer

Abstract

We design a dynamically-weighted state-averaged constrained CASSCF to treat \ul{e}lectrons or \ul{h}oles moving between $n$ molecular fragments (where $n$ can be larger than 2). Within such a so-called eDSCn/hDSCn approach, we consider configurations that are mutually single excitations of each other, and we apply a generalized set of constraints to tailor the method for studying charge transfer problems. The constrained optimization problem is efficiently solved using a DIIS-SQP algorithm, thus maintaining computational efficiency. We demonstrate the method for a finite Su-Schrieffer-Heeger (SSH) chain, successfully reproducing the expected exponential decay of diabatic couplings with distance. When combined with a gradient, the current extension immediately enables efficient nonadiabatic dynamics simulations of complex multi-state charge transfer processes.

Equations

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

$$|{\Psi }_1⟩$$

$$|{\Psi }_2⟩$$

$$|{\Psi }_1⟩$$

$$|{\Psi }_2⟩$$

$$w_1(\Delta E)$$

$$w_2(\Delta E)$$

$$d{E}_{\mathrm{tot}}$$

$$w_1^{\prime }$$

$$w_2^{\prime }$$

$$\mathrm{Tr}({\hat{P}}^D{\hat{P}}^{\mathrm{active}}-{\hat{P}}^A{\hat{P}}^{\mathrm{active}})=0.$$

$$\mathcal{L}=w_1{E}_1+w_2{E}_2-\lambda \mathrm{Tr}({\hat{P}}^D{\hat{P}}^{\mathrm{active}}-{\hat{P}}^A{\hat{P}}^{\mathrm{active}})$$

$$E_{\mathrm{tot}}=\sum _{j=1}^{M^c}{w}_j{E}_j,$$

$$|{\Psi }_j⟩=|1\bar{1},2\bar{2},...,N\bar{N},N+j⟩,$$

$$|{\Psi }_j⟩=|1\bar{1},2\bar{2},...N-j,\overline{N-j},N-j+1,N-j+2,\overline{N-j+2},...,N\bar{N}⟩,$$

$$\mathrm{Tr}\left({\hat{P}}^1{\hat{P}}^{\mathrm{active}}\right)=\mathrm{Tr}\left({\hat{P}}^2{\hat{P}}^{\mathrm{active}}\right)=\cdots =\mathrm{Tr}\left({\hat{P}}^{M^f}{\hat{P}}^{\mathrm{active}}\right).$$

$$\mathrm{Tr}\left(({\hat{P}}^j-{\hat{P}}^{j+1}){\hat{P}}^{\mathrm{active}}\right)=0,j=1\cdots M^f-1.$$

$$w_j=\frac{a_j}{{\sum }_j{a}_j},$$

$$a_j$$

$$\Delta E_j$$

$$d{E}_{\mathrm{tot}}=\sum _j{w}_j^{\prime }d{E}_j,$$

$$w_{j>1}^{\prime }$$

$$w_1^{\prime }$$

$$\mathcal{L}=\sum _{j=1}^{M^c}{w}_j{E}_j-\sum _{j=1}^{M^f-1}{\lambda }_j\mathrm{Tr}\left(({\hat{P}}^j-{\hat{P}}^{j+1}){\hat{P}}^{\mathrm{active}}\right).$$

$$E_j=\frac{1}{2}\mathrm{Tr}\left[({\mathbf{H}}^0+{\mathbf{F}}^{j\alpha }){\mathbf{P}}^0\right]+\frac{1}{2}\mathrm{Tr}\left[({\mathbf{H}}^0+{\mathbf{F}}^{j\beta })({\mathbf{P}}^0+{\mathbf{P}}^j)\right]$$

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

$$\mathrm{Tr}\left[{\mathbf{Q}}^j{\mathbf{P}}^{\mathrm{active}}\right]=0,$$

$${\mathbf{P}}^{\mathrm{active}}=\sum _{j=1}^{M^c}{\mathbf{P}}^j.$$

$$\mathcal{L}=$$

$$-$$

$$\left({\nabla }_{{\mathbf{P}}^0}\mathcal{L}\right){\mathbf{P}}^0$$

$$\left({\nabla }_{{\mathbf{P}}^{\mathbf{k}}}\mathcal{L}\right){\mathbf{P}}^k$$

$${\mathbf{M}}^0$$

$${\mathbf{M}}^{j>0}$$

$${\mathbf{Q}}^0$$

Full text

1 Introduction

Electron transfer through multiple charge centers play a fundamental role in complex physical and chemical processes. For example, multi-site electron transfer governs the rates of electrocatalytic reactions 1, 2, determines the efficiency of multi-junction solar cells 3, and controls charge transport in molecular electronics 4, 5, 6. Obviously, understanding the coupling between electronic states is crucial, as it is known that $(i)$ in the limit of weak coupling, there is incoherent sequential charge transfer between sites and (ii) in the limit of strong coupling, one finds very fast band transport.7, 6, 8, 9, 10 These rules of thumb hold in general, and can be ascertained easily for models of noninteracting electronic Hamiltonians. That being said, however, if we seek to go beyond model Hamiltonians and investigate multisite electron transfer (and ideally electron transport) within an *ab initio *context, we face several obstacles. First, the presence of multiple states introduces competing pathways for charge transfer (CT) and increased quantum effects, such that one cannot treat CT between each pair of charge centers separately. Second, to account for electron-phonon coupling, a meaningful theoretical description must switch smoothly between charge-localized state and delocalized states across multiple sites that may change dramatically with molecular geometry, external fields, and chemical environment. Third, it is critical that we recover balanced electron correlation energies for electrons in different environments while also maintaining efficiency and computational accessibility for dynamics simulations; in particular, it is essential that one treat multiple electronic configurations on an equal footing.

To that end, single-reference methods such as Hartree-Fock (HF) or Density Functional Theory (DFT) are clearly not the appropriate methods of choice, as they fail to capture the multi-reference nature of CT processes; Time-Dependent DFT (TD-DFT) has emerged as a widely-adopted approach for modeling charge transfer processes,11, 12 but the latter treats excited states and the ground state in an unbalanced fashion (and cannot recover the correct topology of a crossing13). Recent work in our group14 (and analogous work from other groups15, 16) to restore balance can indeed work for electron transfer between two sites but not has not yet been generalized to multiple site problems. Multi-reference methods such as Complete Active Space Self-Consistent Field (CASSCF), Multi-Reference Configuration Interaction (MRCI),17, 18 and density matrix renormalization group (DMRG) theory are natural candidates for future work as they can provide accurate descriptions – but these methods do suffer from rapid computational scaling with the size of the system and the active space. Moreover, it is worth noting that working with CASSCF becomes more and more difficult for systems with more and and more electronic states not only as far as memory costs, but also because, in such contexts, one often finds multiple solutions and discontinuous potential energy surfaces.

With this background in mind, in a recent article 19, 20, we took an initial step toward solving the problem above by introducing an electron/hole dynamically-weighted state-averaged constrained CASSCF (eDSC/hDSC) method for studying two-state CT processes in radical systems. By constraining the active orbital space to project equally onto donor and acceptor fragments, the eDSC/hDSC approach was able to generate nearly globally smooth, accurate potential energy surfaces at approximately twice the cost of restricted open-shell Hartree-Fock (ROHF) calculations. In addition, the method was able to correct describe not only weakly but also strongly coupled CT processes (where methods like constrained DFT (cDFT) fail). That being said, the method delineated in Ref. 19 was restricted to a two-site system; and given the discussion above, there would be an enormous benefit (with many applications) if one could extend the eDSC/hDSC approach so as to treat multi-state electron or hole transfer problems in an efficient manner.

To that end, in the present paper, we will show that one can indeed extend the eDSC/hDSC methodology to treat CT through $n$ multiple charge centers / molecular fragments, resulting in an algorithm that we will refer to as eDSCn/hDSCn. Just as for eDSC/hDSC20, the algorithm is able to efficiently treat multiple electronic states and constraints efficiently by using a DIIS-SQP solver for the constrained optimization problem. At the end fo the day, we believe the method opens the door to efficient nonadiabatic dynamics simulations of complex multi-state CT processes, which has broad applications in study of energy conversion, molecular electronics, and biological electron transfer.

An outline of the paper is as follows. In Sec. 2, we briefly review the eDSC/hDSC method. In Sec. 3, we extend the formalism to the case of more than two molecular fragments and configurations. In Sec. 4, we derive the working equations for solving the extended formalism with the DIIS-SQP algorithm. In Sec. 5, we introduce a generalized constraint formulation that provides flexibility in treating different problems. In Sec. 6, we present the results for the finite Su–Schrieffer–Heeger chains. In Sec. 7, we discuss the physical meaning of a convergence trick applied in the algorithm. In Sec. 8, we conclude this manuscript and provide outlook for future research directions.

2 Review of eDSC/hDSC

Before describing the eDSCn/hDSCn algorithm detail, let us review the two-state eDSC/hDSC approach from Ref. 19. For the two state problem, we define the target function for eDSC (exciting the electron)/hDSC (exciting the hole) as

[TABLE]

where $E_{1}$ and $E_{2}$ are energies of two restricted-open-shell electron configurations that both have one unpaired electron. For eDSC, the two configurations are

[TABLE]

and for hDSC, the two configurations are

[TABLE]

$w_{1}$ and $w_{2}$ in Eq. 1 are dynamical weights that depend on a “temperature” parameter $T$ and the energy difference $\Delta E=E_{2}-E_{1}$, where we assume $E_{2}>E_{1}$. For the two-state problem with only one active electron (and effectively zero correlation energy) it is natural to define the weights as:

[TABLE]

The differential of Eq. 1 is

[TABLE]

where

[TABLE]

Note that $w_{1}^{\prime}\neq\frac{\partial{w_{1}}}{\partial{\Delta E}}$.

At the heart of the eDSC/hDSC is the need to suppress any and all local excitations; in so doing, one can find nearly globally smooth surfaces as relevant for charge transfer. In order to achieve such a suppression, we impose a constraint that the vector space spanned by active orbitals (i.e., orbitals $N+1$ and $N+2$ for eDSC and orbitals $N$ and $N+1$ for hDSC) projects equally onto the vector space spanned by orbitals from the donor and acceptor molecular fragments, i.e.,

[TABLE]

where $\hat{P}^{D}$ and $\hat{P}^{A}$ are vector space spanned by orbitals from the donor and acceptor molecular fragments, respectively, and $\hat{P}^{\rm active}$ is the vector space spanned by active orbitals. Eq. 11 defines a constrained optimization problem and one can write the Lagrangian formally as

[TABLE]

where $\lambda$ is the Lagrange multiplier. To solve this problem efficiently, in Ref. 19, we applied the sequential quadratic programming (SQP) method and later introduced a DIIS-SQP algorithm in Ref. 20 to further and dramatically reduce the computational cost to only twice the cost of Restricted Open-shell Hartree Fock (ROHF) calculation.

In the present paper, our goal will now be to extend this theory to study systems with more than two charge centers and active orbitals. In a companion paper, we will further derive the relevant analytic nuclear gradient of adiabatic energies and the derivative coupling.

Before proceeding, we remind the reader that the primary purpose of eDSCn/hDSCn is to study electron/hole transfer in a radical system and we consider only electronic configurations that are mutually single excitations of each other. In other words, we seek the absolutely most efficient and streamlined active space for coupled nuclear-electronic charge transfer simulations, and generating such a platform does necessarily leave a risk that we will ignore important electron-electron correlations. That being said, for large scale simluations, strong approximations are imperative; all of the configurations for eDSCn/hDSCn are defined explicitly in Sec. 3.

3 Extension to $M^{f}$ Fragments with $M^{c}$ Configurations

To probe $M^{c}$ configurations in an odd-number electron system, as above, it is natural define the total energy as the dynamical weighted sum of energies from all configurations, namely,

[TABLE]

where $E_{j}$ is the energy of the $j$th electronic configuration, $\ket{\Psi_{j}}$, and $w_{j}$ is the corresponding weighting factor. Similar to eDSC/hDSC, we will also make a distinction between electron transfer and hole transfer and refer electron transfer as eDSCn and hole transfer as hDSCn. Importantly, the electron configurations considered depend on whether we study the electron transfer or hole transfer. Specifically, for eDSCn, let the system have $2N+1$ electrons, the $j$th electronic configuration is ($j=1\cdots M^{c}$)

[TABLE]

which consists the complete active space for one electron in $M^{c}$ orbitals, often referred to as CAS$(1,M^{c})$ in literature. For hDSCn, let the system have $2N-1$ electrons, the electronic configurations we considered are (again, $j=1\cdots M^{c}$)

[TABLE]

which consists the complete active space for $2M^{c}-1$ electron in $M^{c}$ orbitals, or CAS$(2M^{c}-1,M^{c})$. For both eDSCn and hDSCn, we will refer to the singly occupied orbital in each configuration as the active orbital.

The next step is to develop the appropriate constraints for systems with more than 2 molecular fragments. To suppress local excitations in a system where the charge can travel through $M^{f}$ centers, we find that the best approach is to study $M^{c}=M^{f}$ configurations and apply the constraints that the vector space spanned by all active orbitals projects equally onto each molecular fragments, i.e.,

[TABLE]

where $\hat{P}^{j}$ is the vector space spanned by atomic orbitals from the $j$th molecular fragment and $\hat{P}^{\rm active}$ is the vector space spanned by all active orbitals. Thus, for a system with $M^{f}$ fragments (i.e., charge centers), we have $M^{f}-1$ constraints:

[TABLE]

Finally, the weighting factors $w_{j}$ in Eq. 13 need to be defined in a slightly different way than Eqs. 9-10. We propose the following form (assuming $\ket{\Psi_{1}}$ has the lowest energy)

[TABLE]

where

[TABLE]

Note that $a_{1}=\lim_{\Delta E\rightarrow 0}\frac{1-e^{-\Delta E/T}}{\Delta E/T}=1$. With this definition, the first differential of $E_{\rm tot}$ can be written as

[TABLE]

where

[TABLE]

At this point, one can write a formal Lagrangian for the constrained optimization problem

[TABLE]

4 Solving Constrained Optimization Problem with DIIS-SQP Framework

At the heart of a eDSCn/hDSCn calculation is a constrained optimization problem represented by the Lagrangian in Eq. 24. The framework applied here is an extension to the approach in Ref. 20. For maximum clarity, we will deal with eDSCn in Sec. 4.1; we will review necessary revisions for solving hDSCn problems in Sec. 4.2.

4.1 Case of eDSCn

Let $\bm{P}^{0}$ be the density matrix of orbitals that are always doubly occupied, i.e. the one-electron molecular orbitals $\phi_{1},\phi_{2},\cdots,\phi_{N}$. Let $\bm{P}^{j>0}$ be the density matrix of singly occupied orbitals in configuration $\ket{\Psi_{j}}$, i.e., $\phi_{N+j}$. The energy of the $j$th configuration can be expressed as

[TABLE]

where $\bm{H}^{0}$ is the one-electron Hamiltonian, $\bm{F}^{j\alpha}$, $\bm{F}^{j\beta}$ are spin up and down Fock matrices built from spin densities matrices that $\bm{P}^{j\alpha}\equiv\bm{P}^{0}$, $\bm{P}^{j\beta}\equiv\bm{P}^{0}+\bm{P}^{j}$, respectively. Next, let $\bm{Q}^{j}$ represents the projection operator difference between fragment $j$ and $j+1$, i.e.,

[TABLE]

such that our constraints can be represented most simply as

[TABLE]

where

[TABLE]

One can rewrite the Lagrangian in Eq. 24 using the density matrices $\bm{P}^{j}$ (here we assume an orthonormal basis)

[TABLE]

where $\bm{\epsilon}^{j}$, $\bm{\sigma}^{ij}$, and $\lambda_{j}$ are Lagrangian multipliers. Note that $\bm{\epsilon}^{j}$ and $\bm{\sigma}^{ij}$ are symmetric matrices and with the same dimension as $\bm{P}^{j}$. In Eq. 29, $\lambda_{j}$ enforces the charge transfer constraints we introduced in Eq. 16, $\bm{\epsilon}^{j}$ ensures that the density matrices $\bm{P}^{j}$ are idempotent, and $\bm{\sigma}^{ij}$ ensures that $\bm{P}^{i}$ and $\bm{P}^{j}$ are orthogonal for $i\neq j$ (see Ref. 20). We will now outline how to solve the constrained optimization problem in Eq. 29 in a series of seven steps.

To begin with, one can write

[TABLE]

Note that in the equations above, we define $\bm{\sigma}^{k>j}=\bm{\sigma}^{jk}$ to keep the equations formally symmetric. ($\bm{\sigma}^{kj}$ is defined only for $k<j$ in Eq. 29) 2. 2.

Let us define

[TABLE]

At this point, if we add Eqs. 30 and 31 for all $k$ together, we find:

[TABLE] 3. 3.

Since $\bm{P}^{j}$, $\bm{\epsilon}^{j}$, and $\bm{\sigma}^{kj}$ are all symmetric matrices, the R.H.S. of Eq. 35 is also symmetric, which suggests that the L.H.S. of Eq. 35 must be symmetric at the solution of the constrained optimization problem. In other words, one can use this property to measure the deviation or error from the solution and define a DIIS error vector accordingly:

[TABLE]

Here, $[\bm{A},\bm{B}]\equiv\bm{AB}-\bm{BA}$ is the commutator. 4. 4.

Following Ref. 20, one can establish a DIIS procedure for the optimization of Eq. 29 provided that Eqs. 30-31 can be solved for fixed quantities of $w_{j}^{\prime},\bm{F}^{j\alpha},\bm{F}^{j\beta}$, or equivalently, $\bm{M}^{j}$. In particular, it can be verified that solving Eqs. 30-31 with fixed $\bm{M}^{j}$ is equivalent to minimizing the auxiliary function

[TABLE]

under the same set of constraints in Eq. 29. 5. 5.

To work with the smallest parameter space possible, one can absorb the orthonormality constraints into the parametrization by writing

[TABLE]

where $\bm{C}$ represents molecular orbital coefficients (in an orthonormal basis, $\bm{C}$ is a unitary matrix) and $\bm{N}^{j}$ is an occupation matrix, and parameterize $\bm{C}$ locally by an anti-symmetric matrix $\bm{A}$ (rotation generator), i.e.,

[TABLE]

With this transformation, one can build an auxiliary Lagrangian

[TABLE]

where $\lambda_{k}$ are Lagrangian multipliers. 6. 6.

The constrained minimization problem in Eq. 40 can be solved with SQP, which requires the gradient of the target function and the constraints w.r.t. $\bm{A}$. The energy term gives (note that $\bm{A}$ is anti-symmetric)

[TABLE]

where

[TABLE]

and the constraint term gives

[TABLE]

where $\tilde{\bm{Q}}^{k}=\bm{C}_{0}^{\top}\bm{Q}^{k}\bm{C}_{0}$. Since Eqs 41 and 43 are only valid at $\bm{A}=\bm{0}$, $\bm{C}_{0}$ needs to be updated in each iteration. 7. 7.

For practical optimization problems, it is often desirable to use the diagonal of the Hessian matrix as a preconditioner. One can verify that

[TABLE]

Similar to the algorithm in Ref. 20, let $df_{pq}$ and $dg_{pq}^{k}$ be the gradient of the target function and the $k$th constraint that are required in SQP algorithm, the preconditioner suggests that instead of using Eq. 41 and 43 as the expression for $df_{pq}$ and $dg_{pq}^{k}$, one can instead use

[TABLE]

The analysis above highlights the key conceptual points needed for optimization; the remaining procedure is identical to the algorithm proposed in Ref. 20. For the convenience of the reader, we provide a flowchart in Sec. 4.4.

4.2 Case of hDSCn

The same procedure as above can be applied to hDSCn after redefining some quantities. To be specific, different from Eq. 25, the energy expressions for configurations in the case of hDSCn can be written as

[TABLE]

where $\bm{P}^{0}$ still represents the density matrix from molecular orbitals $\phi_{1},\phi_{2},\cdots,\phi_{N}$. Note that, though, we define $N$ by counting $2N-1$ total electrons in hDSCn (see Eq. 15) rather than $2N+1$ total electrons in eDSCn (see Eq. 14), such that $\bm{P}^{0}$ differs when applying eDSCn or hDSCn even to the same molecule with the same number of electrons. $\bm{P}^{j>0}$ still represents the density matrix from singly occupied orbitals in configuration $\ket{\Psi_{j}}$, but this molecular orbital is now $\phi_{N+1-j}$. $\bm{H}^{0}$ is still the one-electron Hamiltonian, $\bm{F}^{j\alpha}$, $\bm{F}^{j\beta}$ are spin up and down Fock matrices built from the spin densities $\bm{P}^{j\alpha}\equiv\bm{P}^{0}$, $\bm{P}^{j\beta}\equiv\bm{P}^{0}-\bm{P}^{j}$. As a result, Eqs. 32-33 should be replaced by

[TABLE]

With these revised definition of $\bm{M}^{j}$, the rest of the procedure remains the same as eDSCn in Sec. 4.1.

4.3 Avoiding the saddle point

As discussed in Ref. 20, there is one extra symmetry that needs to be addressed to avoid converging to the saddle points during the optimization procedure. Specifically, it can be shown that for eDSCn, the following identity holds:

[TABLE]

where

[TABLE]

Similarly, for hDSCn

[TABLE]

To avoid the saddle point caused by this symmetry, one can flip the sign of half of these matrix elements. Specifically, for eDSCn, for every $0<j<k\leq M^{c}$, we apply

[TABLE]

and for hDSCn, for every $0<k<j\leq M^{c}$, we apply

[TABLE]

This procedure should be done before building $\bm{M}^{j}$ in Eqs. 32-33 and Eqs. 48-49.

4.4 Flowchart

For the convenience of the reader, let us now provide a flowchart of the final procedure. This flowchart is practically identical to the flowchart in Ref. 20

1. Step 1:

Compute $\{\bm{Q}^{j>0}\}$ from Eq. 26, choose an initial guess $\bm{C}_{\rm ini}$ for $\bm{C}$ and initial guesses for $\{\lambda_{j}\}$, and choose a DIIS threshold $T_{D}$. 2. Step 2:

Build density matrices $\{\bm{P}^{j}\}$ from $\bm{C}$ using Eq. 38. 3. Step 3:

Build $\{w_{j}^{\prime}\}$, $\{\bm{F}^{j\alpha}\}$, and $\{\bm{F}^{j\beta}\}$ from $\{\bm{P}^{j}\}$ using Eqs. 22-23 and 25, then build $\{\bm{M}^{j}\}$ using Eqs. 32-33. 4. Step 4:

Compute the initial DIIS error vector $\bm{V}^{0}$ using Eq. 36, i.e.,

[TABLE] 5. Step 5:

If norm $(\bm{V}^{0})<T_{D}$, jump to Step 15. Otherwise, set the initial DIIS variable $\bm{A}^{0}$ to $\bm{0}$. 6. Step 6:

Solve the nSCF problem for the new $\bm{C}$ and $\{\lambda_{j}\}$ using SQP with Eqs. 45-46 given $\{\bm{M}^{j}\}$, and $\{\bm{Q}^{j>0}\}$. The nSCF threshold for the norm of the gradient and the constraint is set to norm $(\bm{V}^{0})/100$. 7. Step 7:

(DIIS iteration starts. Let the iteration counter $n$ starts from 1.) Build density matrices $\{\bm{P}^{j}\}$ from the new $\bm{C}$ using Eq. 38. 8. Step 8:

Build $\{w_{j}^{\prime}\}$, $\{\bm{F}^{j\alpha}\}$, and $\{\bm{F}^{j\beta}\}$ from $\{\bm{P}^{j}\}$ using Eqs. 22-23 and 25, then build $\{\bm{M}^{j}\}$ using Eqs. 32-33. 9. Step 9:

Calculate the DIIS error vector at iteration $n$:

[TABLE] 10. Step 10:

Solve the nSCF problem for the new $\bm{C}$ and $\lambda$ using SQP with Eqs. 45-46 given $\{\bm{M}^{j}\}$, and $\{\bm{Q}^{j>0}\}$. The nSCF threshold for the norm of the gradient and the constraint is set to norm $(\bm{V}^{n})/100$. 11. Step 11:

Calculate the DIIS variable at iteration $n$:

[TABLE] 12. Step 12:

If norm $(\bm{V}^{n})<T_{D}$, exit the iteration and jump to Step 15. Otherwise, use DIIS algorithm to build the interpolated DIIS variable $\bm{A}^{\rm next}$ from DIIS errors $\{\bm{V}^{i=0\cdots n}\}$ and previous DIIS variables $\{\bm{A}^{i=0\cdots n}\}$. 13. Step 13:

Calculate the new $\bm{C}$ as

[TABLE] 14. Step 14:

Return to Step 7. 15. Step 15:

DIIS convergence achieved.

5 A More General Form of Constraints

Before presenting results for the method above, we note that the framework presented above around the constraints in Eqs. 16 and 17, can be applied to a more general set of constraints. In particular, for constraints written in the form of Eq. 27, the full procedure discussed in Sec. 4.1 and 4.2 remains unchanged if one generalizes the kernel of the $j$th constraint to be

[TABLE]

Note that Eq. 59 reduces to the constraints defined in Sec. 4.1 and 4.2 if one fixes:

[TABLE]

6 Numerical results

We have applied the method above to a finite polyacetylene CH2(CH)n(CH2)+ cation chain, designed to mimic the famous SSH model problem. Here, we define each fragment to be a carbon dimer together with the relevant connected hydrogen atoms. A scheme of such a partition for the SSH chain with 14 carbon atoms (which will be referred to as ssh-c14) is shown in Fig. 1, where we have 7 fragments in total. For a SSH chain with $2M$ carbon atoms (ssh-2M), we partition the molecule into $M$ fragments in the similar way, and we include $M$ active orbitals within a hDSCn (hole transfer) approach.

Though perhaps not surprising to those who have worked with polyacetylene chains21, 22, 23, we have found that working with polyacetylene chains (as in Fig. 1) requires a bit of care. In particular, we have found that Eq. 16 is not in fact appropriate for this model: if we impose a constraint dictating that the vector space spanned by all active orbitals projects equally onto each molecular fragments, we find that hDSCn will ultimately converge to unwanted active orbitals due to boundary effect on the finite chain. Nevertheless, a meaningful, well-balanced set of active orbitals can be defined by requiring that vector space spanned by all active orbitals projects equally onto the first fragment and the last fragment, then projects equally onto the second fragment and the second last fragment, and so forth. For example, for ssh-c14, we impose three constraints of the form (using Eq. 59):

[TABLE]

Below, we report results for analogous SSH chains with $8,10,12,...,26,28$ carbon atoms. Fig. 2a-g shows the density plot of the active orbitals.

Next, in order to connect the present calculations to a true SSH model Hamiltonian, we applied Boys diabatization scheme 24 to the active orbitals. Figs. 2h-n portrays the diabatized orbitals. The localized orbitals in Figs. 2h-n represent the carbon-carbon $\pi$-bond, as expected. Finally, the most essential quantities to compute are the diabatic couplings between diabatic orbitals. To that end, our results are presented in Fig. 3. In Fig. 3a, we show the couplings between the nearest and the second nearest diabatic orbitals as a function of the length of the SSH chain; as one would hope, both quantities reach a constant plateau as the length of the SSH chain increases. In Fig. 3b, we further plot the logarithm of the diabatic coupling versus the rank of nearest neighbor for ssh-c28. The reasonably well linear relation shown here suggests the expected exponential decay of the diabatic coupling as a function of the distance. The $\beta$ value is 0.36 Å*-1*.

At the end of the day, it is clear that the present approach has delivered a meaningful means to parameterize an SSH model. Moreover, all of the data here arises from a smooth parameterization; one can easily begin to deform the nuclear geometry and investigate the resulting changes in diabatic couplings. To that end, in a companion paper, we will derive the relevant nuclear gradients and derivative couplings and begin exploring eDSC/hDSC potential energy surfaces.

7 Discussion: A Trick for Convergence

Before concluding this manuscript, a few words are appropriate regarding the flipping sign technique introduced in Sec. 4.3. One cannot emphasize enough that the current algorithm quickly loses its efficiency without such a sign flip and, as such, the physical interpretation of this procedure deserves further elaboration.

From a bird’s eye point of view, the need for sign flip arises from the specific form of the configuration interaction (CI) Hamiltonian within eDSCn/hDSCn. Note that, for a standard CASSCF procedure, one usually optimizes the orbitals first and then diagonalizes the CI Hamiltonian to find the correct CI coefficients. For our purposes, however, we have no such two-step procedure, but rather a one-shot optimization of the orbitals. Let us index the active orbital of the $j$th configuration as $a_{j}$. Then, the off-diagonal of the CI Hamiltonian gives

[TABLE]

At the same time, however, Eqs. 50 and 52 can be written together as

[TABLE]

Thus, the correct (converged) CASSCF solution therefore requires that several matrix elements of the Fock matrices be zero upon convergence.

However, consider now the DIIS error vector in Eq. 36; this vector is essentially the gradient of the dynamically weighted sum of the energies (plus the Lagrange multiplier weighted sum of the constraints) w.r.t. the rotation generator. Upon convergence, the DIIS error vector should satisfy:

[TABLE]

Given Eq. 65, the above equation has two solutions, namely either $w_{j}^{\prime}=w_{k}^{\prime}$ or $\left(\tilde{\bm{F}}^{j\beta}\right)_{a_{j},a_{k}}=0$. Note that the latter solution is the correct CASSCF solution while the former solution is the saddle point that we aim to avoid. By flipping the sign for one of the term in each pair, the saddle point solution is removed and we always reach the correct CASSCF solution.

8 Conclusion

We have successfully extended the eDSC/hDSC methodology to treat charge transfer through multiple molecular fragments, developing a so-called eDSCn/hDSCn for systems with $n$ (i.e. more than two) charge centers. The method employs dynamically-weighted state averaging with constraints in a very general form. The DIIS-SQP algorithm developed previously20 can be extended and remains extremely efficient for this constrained optimization problem. Moreover, the flexible constraint formulation allows adaptation to different molecular topologies and boundary conditions, an application to SSH chains demonstrates the method’s ability to treat systems with strong couplings between charge localized diabatic states. This extension opens new possibilities for studying complex charge transfer processes in extended systems. In particular, by prioritizing computational efficiency, the present approach is suitable for nonadiabatic dynamics simulations of CT across complex extended systems, and in a companion paper , we develop analytic gradients and derivative couplings for dynamics applications.

Looking forward, beyond applications, given the current intense focus on chiral induced spin selectivity25 and spin-dependent electron transfer, future theoretical developments should focus on extensions of eDSC/hDSC to include complex orbitals with spin-orbital couplings (SOC) and a phase space Hamiltonian26, 27, 28; in particular, if one can successfully incorporate SOC and a phase space couplings, one hopes to learn a great deal about angular momentum conservation in the context of electron transfer. Finally, while the present approach is limited to one electron or one hole, one would like to explore closed-shell systems with double excitations. While the computational cost would presumably increase dramatically when considering double excitations, some of our preliminary work has suggested that one can dramatically speed up constrained CASSCF(2,2) calculations by folding CI diagonalization into the DIIS-SQP orbital optimization procedure for one-shot optimization cycle (as above). There are indeed still many open avenues for exploration as we seek the fastest means possible to coupled coupled nuclear-electronic charge transfer dynamics.

9 Acknowledgment