# Response to “Response to ‘Comment on Theoretical examination of QED
Hamiltonian in relativistic molecular orbital theory”’ [J. Chem. Phys. 160,
187102 (2024)]

Wenjian Liu<EMAIL_ADDRESS>Qingdao Institute for Theoretical and
Computational Sciences, School of Chemistry and Chemical Engineering, Shandong
University, Qingdao, Shandong 266237, P. R. China

The Response of Inoue and coworkers[1] to my Comment[2] on their Paper I[3]
clarifies some points, but also raises some new issues that require further
clarifications.

To begin with, it should be noted that any second quantization must start with
a definition of the fermion field operator $\hat{\phi}$ in terms of a vacuum
and corresponding annihilation and creation operators, such that the resulting
Hamiltonian takes the whatsoever vacuum as its ground state of zero energy.
For instance, if we start with

$\displaystyle\hat{\phi}$ $\displaystyle=a_{p}\psi_{p},\quad
a_{p}|\mathrm{vac}\rangle=0,\quad p\in\mbox{PES, NES},$ (1)

we will have the following Hamiltonian (in the notation of the Comment[2])

$\displaystyle\mathcal{H}$
$\displaystyle=D_{p}^{q}a^{p}_{q}+\frac{1}{2}g_{pr}^{qs}a^{pr}_{qs},\quad
p,q,r,s\in\mbox{PES, NES},$ (2)

which is already normal ordered with respect to $|\mathrm{vac}\rangle$.
However, this Hamiltonian does not distinguish the empty from the filled Dirac
picture. Since the empty Dirac picture must be abandoned (for it implies that
no atom would be stable), we ought to incorporate the field Dirac picture, by
introducing a reference $|0;\tilde{N}\rangle$
($=\Pi_{\tilde{i}}^{\tilde{N}}a^{\tilde{i}}|0;\tilde{0}\rangle$) built up with
zero positive energy states (PES) and $\tilde{N}$ ($\rightarrow\infty$)
negative energy states (NES). The Hamiltonian (2) must then be normal ordered
with respect to $|0;\tilde{N}\rangle$, so as to obtain the physical
Hamiltonian $H_{n}^{\mathrm{QED(MO)}}$[4, 5]

$\displaystyle\mathcal{H}$
$\displaystyle=H_{n}^{\mathrm{QED(MO)}}+C_{n}^{\mathrm{QED(MO)}},$ (3)
$\displaystyle H_{n}^{\mathrm{QED(MO)}}$
$\displaystyle=H^{\mathrm{FS(MO)}}_{n}+Q_{p}^{q}\\{a^{p}_{q}\\}_{n},\quad
p,q,r,s\in\mbox{PES, NES},$ (4) $\displaystyle H^{\mathrm{FS(MO)}}_{n}$
$\displaystyle=D_{p}^{q}\\{a^{p}_{q}\\}_{n}+\frac{1}{2}g_{pr}^{qs}\\{a^{pr}_{qs}\\}_{n},$
(5) $\displaystyle Q_{p}^{q}$
$\displaystyle=\bar{g}^{qs}_{pr}\acontraction[0.5ex]{}{a^{r}}{}{a_{s}}a^{r}a_{s},\quad\bar{g}_{pr}^{qs}=g^{qs}_{pr}-g^{sq}_{pr},$
(6) $\displaystyle C_{n}^{\mathrm{QED(MO)}}$ $\displaystyle=\langle
0;\tilde{N}|H|0;\tilde{N}\rangle$
$\displaystyle=D_{p}^{q}\acontraction[0.5ex]{}{a^{p}}{}{a_{q}}a^{p}a_{q}+\frac{1}{2}\bar{g}_{pr}^{qs}\acontraction[0.5ex]{}{a^{p}}{}{a_{q}}a^{p}a_{q}\acontraction[0.5ex]{}{a^{r}}{}{a_{s}}a^{r}a_{s},$
(7)

where the contraction can be one of the three cases

$\displaystyle\acontraction[0.5ex]{}{a^{p}}{}{a_{q}}a^{p}a_{q}$
$\displaystyle=\begin{cases}0\quad\mbox{ (constantly null contraction
(CNC))};\\\ \langle 0,\tilde{N}|a^{p}_{q}|0,\tilde{N}\rangle\quad\mbox{
(conventional contraction (CC))};\\\ \frac{1}{2}\langle
0,\tilde{N}|[a^{p},a_{q}]|0,\tilde{N}\rangle\quad\mbox{ (charge-conjugated
contraction (CCC))}.\end{cases}$ (8)

Note that the CC is unphysical[6] and is included here only for later use. By
making use[1] of the elementary anticommutation relation (ACR)

$a^{p}a^{q}+a^{q}a^{p}=0,$ (9a) $a_{p}a_{q}+a_{q}a_{p}=0,$ (9b)
$a^{p}a_{q}+a_{q}a^{p}=\delta_{pq},$ (9c)

the (generic) Hamiltonian[4, 5] $H_{n}^{\mathrm{QED(MO)}}$ (4) was rewritten
in Paper I[3] in terms of unordered operators

$O_{p}^{q}\\{a^{p}_{q}\\}_{n}=O_{p}^{q}a^{p}_{q}-O_{\tilde{i}}^{\tilde{i}},\quad
O=D,Q,$ (10a)
$\frac{1}{2}g_{pr}^{qs}\\{a^{pr}_{qs}\\}_{n}=\frac{1}{2}g_{pr}^{qs}a^{pr}_{qs}-\bar{g}_{p\tilde{j}}^{q\tilde{j}}a^{p}_{q}+\frac{1}{2}\bar{g}_{\tilde{i}\tilde{j}}^{\tilde{i}\tilde{j}},$
(10b)

so as to obtain

$\displaystyle H^{\mathrm{QED(MO)}}_{n}$
$\displaystyle=\mathcal{H}+(Q_{p}^{q}-\bar{g}_{p\tilde{j}}^{q\tilde{j}})a^{p}_{q}+\bar{C}_{n}^{\mathrm{QED(MO)}},$
(11) $\displaystyle\bar{C}_{n}^{\mathrm{QED(MO)}}$
$\displaystyle=-(Q+D)_{\tilde{i}}^{\tilde{i}}+\frac{1}{2}\bar{g}_{\tilde{i}\tilde{j}}^{\tilde{i}\tilde{j}},$
(12)

where

$\displaystyle Q_{p}^{q}$ $\displaystyle=\begin{cases}0\quad\mbox{ (CNC)},\\\
\bar{g}_{p\tilde{j}}^{q\tilde{j}}\quad\mbox{ (CC)},\\\
-\frac{1}{2}\bar{g}_{ps}^{qs}{\rm{sgn}}({\epsilon_{s}})\quad\mbox{
(CCC)}.\end{cases}$ (13)

A close inspection reveals that the Hamiltonian $\mathcal{H}$ in Eq. (11) must
adopt the CC [which is in line with Eq. (9)] in order to be consistent.
Therefore, Eq. (11) had better be rewritten as

$\displaystyle H^{\mathrm{QED(MO)}}_{n(u)}$
$\displaystyle=\mathcal{H}_{\mathrm{CC}}+(Q_{p}^{q}-\bar{g}_{p\tilde{j}}^{q\tilde{j}})a^{p}_{q}+\bar{C}_{n}^{\mathrm{QED(MO)}}$
(14)

in order to be precise. As such, Eq. (11)/(14) is indeed a hybrid[2] of “two
pictures”, CNC/CCC vs. CC. Note in passing that Eq. (11) can also be
understood as a reexpression of $H^{\mathrm{QED(MO-X)}}_{n}$ (X = CNC, CC,
CCC) in terms of $H^{\mathrm{QED(MO-CC)}}_{n}$, i.e.,

$\displaystyle H^{\mathrm{QED(MO-X)}}_{n}=H^{\mathrm{QED(MO-
CC)}}_{n}+(Q-\bar{g}_{p\tilde{j}}^{q\tilde{j}})\\{a^{p}_{q}\\}_{n},$ (15)

which is then converted to the formally unordered form (14) (indicated by the
subscript $n(u)$).

Since quantities like $D_{\tilde{i}}^{\tilde{i}}$,
$\bar{g}_{\tilde{i}\tilde{j}}^{\tilde{i}\tilde{j}}$, and especially
$\bar{g}_{p\tilde{j}}^{q\tilde{j}}a^{p}_{q}$ are all divergent (due to the
existence of infinitely many NES $\\{\tilde{i},\tilde{j}\\}$), they cannot be
evaluated numerically. Rather, they must be canceled out analytically when
calculating the physical energy of an $N$-electron system according to

$\displaystyle E$
$\displaystyle=\langle\Psi(N;\tilde{N})|H^{\mathrm{QED(MO)}}_{n(u)}|\Psi(N;\tilde{N})\rangle-\langle\Psi(0;\tilde{N})|H^{\mathrm{QED(MO)}}_{n(u)}|\Psi(0;\tilde{N})\rangle.$
(16)

Here, the (normalized) wave function $|\Psi(0;\tilde{N})\rangle$ represents a
polarizable vacuum, with $|0;\tilde{N}\rangle$ as its zeroth order[5]. As
shown by Inoue and coworkers[1], the divergent terms appearing in Eqs. (13),
(16), (33), (81), and (99) in Paper I[3] do disappear after working out the
matrix elements explicitly. Therefore, the criticisms on these equations made
in the Comment[2] should be eliminated. Still, however, this twisted form (14)
of the Hamiltonians is unpleasant, for it looks like that the divergent terms,
which are not present in the elegant expressions of $H^{\mathrm{QED(MO-
CNC)}}_{n}$/$H^{\mathrm{FS(MO)}}_{n}$ (5) and $H^{\mathrm{QED(MO-CCC)}}_{n}$
(4), are first extracted out deliberately and then canceled out carefully. In
our opinion, it can be viewed at most as an intermediate quantity (instead of
a proper form of the Hamiltonians) for the evaluation of the energy (16).
Actually, it is not needed even for the calculation of the energy, for which
can actually be calculated as[7, 5]

$\displaystyle E$
$\displaystyle=\langle\Psi(N;\tilde{N})|H^{\mathrm{QED(MO)}}_{n}|\Psi(N;\tilde{N})\rangle-\langle\Psi(0;\tilde{N})|H^{\mathrm{QED(MO)}}_{n}|\Psi(0;\tilde{N})\rangle$
(17)
$\displaystyle=\langle\Psi(N;\tilde{N})|\mathcal{H}-C^{\mathrm{QED(MO)}}_{n}|\Psi(N;\tilde{N})\rangle-\langle\Psi(0;\tilde{N})|\mathcal{H}-C^{\mathrm{QED(MO)}}_{n}|\Psi(0;\tilde{N})\rangle$
(18)
$\displaystyle=\langle\Psi(N;\tilde{N})|\mathcal{H}|\Psi(N;\tilde{N})\rangle-\langle\Psi(0;\tilde{N})|\mathcal{H}|\Psi(0;\tilde{N})\rangle.$
(19)

It is just that the same CNC/CCC contraction should be followed when going
from Eq. (17) to (19). To the best of our knowledge, this ansatz was first
proposed in Ref. 7 and presented more clearly in Ref. 5. It can readily
checked that the CC even does not yield the correct first order energy
$E^{(1)}_{\mathrm{CC}}$, viz.,

$\displaystyle E^{(1)}_{\mathrm{CC}}$ $\displaystyle=\langle
N;\tilde{N}|H^{\mathrm{QED(MO-CC)}}_{n}|N;\tilde{N}\rangle$
$\displaystyle=D_{i}^{i}+\frac{1}{2}\bar{g}_{ij}^{ij}+\bar{g}_{i\tilde{j}}^{i\tilde{j}}\rightarrow+\infty,\quad
i,j\in\mbox{ PES},\tilde{j}\in\mbox{ NES},$ (20)

and should hence be rejected from the outset.

At variance with Eq. (17)/(19), the energy $E$ of an $N$-electron system can
also be calculated either as[4]

$\displaystyle E$ $\displaystyle=\langle
N;\tilde{N}|H^{\mathrm{QED(MO)}}_{n}|\bar{\Psi}(N;\tilde{N})\rangle,\quad|\bar{\Psi}(N;\tilde{N})\rangle=\Omega_{n}|N;\tilde{N}\rangle,$
(21)

or as[5]

$\displaystyle E$
$\displaystyle=\langle\bar{\Psi}(N;\tilde{N})|H^{\mathrm{QED(MO)}}_{n}|\bar{\Psi}(N;\tilde{N})\rangle$
(22)
$\displaystyle\equiv\langle\bar{\Psi}(N;\tilde{N})|\mathcal{H}|\bar{\Psi}(N;\tilde{N})\rangle-C^{\mathrm{QED(MO)}}_{n}.$
(23)

The two differ only in normalization (intermediate vs. unitary normalization).
Again, the same CNC/CCC contraction should be followed when going from Eq.
(22) to (23). A subtle but important point lies in that the wave function
$|\bar{\Psi}(N;\tilde{N})\rangle$ is here generated by the action of a normal-
ordered wave operator $\Omega_{n}$ on the non-interacting reference
$|N;\tilde{N}\rangle$. As such, a polarizable vacuum is not needed here (for
detailed derivations, see Section IIB.2 and Appendix A in Ref. 4). The so-
calculated energy is identical to that obtained by Eq. (19)/(16) and is in
termwise agreement with that by the S-matrix formulation of QED[4]. With this
more transparent presentation (as compared to Eq. (27) in Ref. 5 or Eqs. (48)
and (49) in the Comment[2]), it should be clear that the remark (on Eq. (23))
that ”the terms subtracted from the referenced Hamiltonians are single Slater
determinants and cannot remove total energy divergence caused by the
generalized electron correlation”[3] does not hold true.

As far as the nonrelativistic limit of $H^{\mathrm{QED(MO-CCC)}}_{n}$ (4) is
concerned, Inoue and coworkers obviously did not recognize the physical
meaning of the $Q$-potential (6)/(13) arising from the CCC (8). As emphasized
from its advent[4] (see also Ref. 6 for a thorough discussion), the direct and
exchange parts of $Q$ (13) represent the vacuum polarization (VP) and electron
self-energy (ESE), respectively (NB: the frequency-dependent transverse-photon
interaction has to be further added to the exchange part to make the ESE
complete[5]). It is then textbook knowledge that both the VP and ESE have to
be regularized/renormalized before they can actually be used (see Ref. 8 for
numerical demonstrations). After regularization/renomalization, the
$Q$-potential is at most of $\mathcal{O}(Z^{3}\alpha^{3})$ and therefore
vanishes in the nonrelativistic limit. That is, the two terms in Eq. (10) in
the Response[1] are divergent individually but tend to cancel each other, and
can hence be renomalized away. This is an essential feature of QED. Speaking
of “QED Hamiltonian” in the absence of the VP-ESE (_genuine QED effect_[4]) is
hardly meaningful. What were reported in Paper I[3] are actually relativistic
configuration interaction and many-body perturbation theories of composite
systems of $N$ electrons and $\tilde{M}$ (real) positrons governed by the
instantaneous Coulomb interaction, instead of the initially anticipated[2]
contribution of the NES (virtual positrons) to the correlation of $N$
electrons, a kind of _derived QED effect_[4]. For instance, Eqs. (B1), (B4),
and B(6) in the Response [1] all go back to the corresponding no-pair
expressions for $N$ electrons by setting the number of positrons to zero. It
should be restated that Eq. (95) in Paper I[3] misses by construction the one-
body term $E_{1}^{(2)}$ in Eq. (52) in the Comment[2]. Likewise, their “QED-
based” Dirac-Hartree-Fock (DHF) theory (see Section IIC.1 in Paper I) is also
merely a relativistic but non-QED mean-field theory[9] of $N$ electrons and
$\tilde{M}$ positrons. The genuine QED mean-field theories for electrons only
and for both electrons and positrons were presented in Ref. 5 (see Eq.
(54)/(81) or (85) therein) and Ref. [6, 10], respectively. In such genuine QED
mean-field theories, the PES and NES, whether occupied or not, are all
coupled, thereby fundamentally different from the relativistic but non-QED
counterparts.

Inoue and coworkers also made an interesting observation[1]: in contrast to
the CC, both the CNC and CCC do not satisfy the basic ACR of fermion
operators(9). At first glance, this might lead to that the physically
incorrect CC is good, but the physically correct CCC/CNC is bad! To resolve
this apparent self-contradiction, they separated the CNC/CCC from the ACR (or
equivalently the CC), and took the former only as a formal-theoretical process
for constructing the Hamiltonians. This route is certainly not incorrect.
However, it is not a proper interpretation. In sharp contrast, the
introduction of the CCC[4] is firmly founded, not merely a formal
(mathematical) step. As scrutinized in depth in Ref. 6, the introduction of a
filled Dirac sea of electrons (which is not part of the Dirac equation itself)
is a must but not yet complete. By charge conjugation symmetry (which is
indeed a property of the Dirac equation), there exists also a filled Dirac sea
of positrons (which is again not part of the Dirac equation itself). The two
seas are coexistent and equivalent, and should hence be averaged with an equal
weight of one half, thereby leading naturally to the CCC[4]. Just like that
the ACR is dictated by fermi statistics of nonrelativistic or relativistic
fermions, the CCC is dictated by charge conjugation symmetry of relativistic
fermions. Both are fundamental laws of physics that have to be imposed from
the outside of quantum mechanical equations, so as to render the latter
physically correct (recalling that the Schrödinger equation can also describe
bosons, provided that boson statistics is mposed). Note in passing that the
CCC is just the time-independent analog of the symmetric-in-time, equal-time
contraction (ETC) of time-dependent fermion operators[11], which is embodied
automatically in the Feynman fermion propagator[12], again not merely a formal
ingredient. The correct four-current for electrons (and positrons), and hence
the VP-ESE represented by the $Q$-potential (6)/(13), can only be obtained by
the CCC/ETC (see Eqs. (36)–(58) in Ref. 6). Yet, it should be kept in mind
that the CCC is to be applied only when normal ordering with respect to the
NES (virtual positrons). In contrast, the CC (or equivalently the ACR) should
still be applied when further normal ordering with respect to the occupied
PES[5], where charge conjugation is irrelevant. It is in this sense that the
$Q$-potential arising from the CCC (zero from the CNC) should be interpreted
as an integral part of the Hamiltonian, as implied in Eqs. (11), (14), and
(15). In other words, being imposed from the outset, the CCC
$\acontraction[0.5ex]{}{a^{p}}{}{a_{q}}a^{p}a_{q}$ over the vacuum
$|0;\tilde{N}\rangle$ should be viewed as an integral part of the normal
ordering process. It can also be shown that the CCC arises naturally from the
averaging of the Hamiltonians for electrons and positrons moving in the same
external field (which is a must, as shown in Appendix A).

Since the Hamiltonian $H^{\mathrm{QED(MO-CNC)}}$ misses by construction the
_genuine QED effect_[4], it had better go back to its original name, “no-
photon (non-QED) Fock space Hamiltonian”, as advocated by Kutzelingg more than
a decade ago[13].

In summary, the major conclusion in the Comment[2] remains unaltered: it is
the CCC-based Hamiltonian $H^{\mathrm{QED(MO-CCC)}}$, instead of the CNC-based
Hamiltonian $H^{\mathrm{QED(MO-CNC)}}$/$H^{\mathrm{FS(MO)}}$, that is a
genuine QED Hamiltonian, and should hence be recommended as the basis of the
emerging field of “molecular QED”[5].

## Acknowledgement

This work was supported by the National Natural Science Foundation of China
(Grant No. 22373057) and Mount Tai Scholar Climbing Project of Shandong
Province.

## Data Availability Statement

The data that supports the findings of this study is available within the
article.

##  Appendix A Hamiltonian quantized with charge conjugated fields

The charge conjugation is defined as

$\hat{C}=\mathbf{C}_{0}\hat{K}_{0},\quad\hat{C}^{\dagger}=\hat{C}^{-1}=\hat{C},\quad\mathbf{C}_{0}=-\mathbbm{i}\alpha_{y}\beta=\begin{pmatrix}\mathbf{0}_{2}&i\sigma_{y}\\\
-i\sigma_{y}&\mathbf{0}_{2}\end{pmatrix}=\mathbf{C}_{0}^{\dagger}=\mathbf{C}_{0}^{-1},$
(24)

which transforms the one-electron Dirac equation

$\displaystyle D\psi_{p}$ $\displaystyle=\epsilon_{p}\psi_{p},$ (25)
$\displaystyle D$ $\displaystyle=D_{0}+q\phi_{\mathrm{ext}}(\bm{r}),\quad
q=-1,$ (26) $\displaystyle D_{0}$ $\displaystyle=c\bm{\alpha}\cdot\bm{p}+\beta
mc^{2}$ (27)

to that of a positron

$\displaystyle D^{C}\psi^{C}_{p}$
$\displaystyle=\epsilon_{p}^{C}\psi^{C}_{p},$ (28) $\displaystyle D^{C}$
$\displaystyle=-\hat{C}D\hat{C}^{-1}=-\mathbf{C}_{0}^{\dagger}D^{*}\mathbf{C}_{0}=D_{0}-q\phi_{\mathrm{ext}}(\bm{r}),$
(29) $\displaystyle\psi_{p}^{C}$
$\displaystyle=\hat{C}\psi_{p}=\mathbf{C}_{0}\psi^{*}_{p},\quad\epsilon_{p}^{C}=-\epsilon_{p}.$
(30)

Likewise, the four-component Dirac field operator $\hat{\phi}(\bm{r})$ (in the
particle-hole picture)

$\displaystyle\hat{\phi}_{\sigma}(\bm{r})$
$\displaystyle=b_{p}\psi_{p\sigma}(\bm{r})+b^{\tilde{p}}\psi_{\tilde{p}\sigma}(\bm{r}),\quad
p\in\mbox{ PES},\quad\tilde{p}\in\mbox{ NES},\quad\sigma\in[1,4]$ (31)

will be transformed by charge conjugation to

$\displaystyle\hat{\phi}^{C}(\bm{r})$
$\displaystyle=\mathbf{C}_{0}\hat{\phi}^{{\dagger}T}(\bm{r})$ (32)
$\displaystyle=b_{\tilde{p}}\psi^{C}_{\tilde{p}}(\bm{r})+b^{p}\psi_{p}^{C}(\bm{r}),$
(33)

where the first and second terms annihilate a positive-energy positron (NB:
$\epsilon_{\tilde{p}}^{C}>0$) and create a positive-energy electron,
respectively, in accordance with the fact that the first and second terms of
$\hat{\phi}(\bm{r})$ (31) annihilate a positive-energy electron and create a
positive-energy positron, respectively. That is, $\hat{\phi}^{C}(\bm{r})$ (33)
would be the starting quantized Dirac field, had we lived in the world of
antiparticles.

Both $\hat{\phi}(\bm{r})$ (31) and $\hat{\phi}^{C}(\bm{r})$ (33) are
associated with the vacuum $|\mathrm{vac}\rangle=|0;\tilde{0}\rangle$. To
expedite algebraic manipulations, we can rewrite them as[6]

$\displaystyle\hat{\phi}(\bm{r})$
$\displaystyle=a_{p}\psi_{p}(\bm{r})+a_{\tilde{p}}\psi_{\tilde{p}}(\bm{r}),\quad\epsilon_{p}>0,\quad\epsilon_{\tilde{p}}<0,$
(34) $\displaystyle\hat{\phi}^{C}(\bm{r})$
$\displaystyle=a_{\tilde{p}}\psi^{C}_{\tilde{p}}(\bm{r})+a_{p}\psi_{p}^{C}(\bm{r}),\quad\epsilon_{\tilde{p}}^{C}>0,\quad\epsilon_{p}^{C}<0,$
(35)

by taking $|0;\tilde{N}\rangle$ (Dirac sea of electrons) and
$|0_{e^{+}};\tilde{N}_{e^{+}}\rangle$ (Dirac sea of positrons) as the vacua,
respectively.

With the above background, we can calculate

$\displaystyle h^{C}$
$\displaystyle=\int\hat{\phi}^{C{\dagger}}D^{C}\hat{\phi}^{C}d\tau=\int(\mathbf{C}_{0}\hat{\phi}^{{\dagger}T})^{\dagger}D^{C}(\mathbf{C}_{0}\hat{\phi}^{{\dagger}T})d\tau=\int\hat{\phi}^{T}(\mathbf{C}_{0}^{\dagger}D^{C}\mathbf{C}_{0})\hat{\phi}^{{\dagger}T}d\tau$
$\displaystyle=-\int\hat{\phi}^{T}D^{*}\hat{\phi}^{{\dagger}T}d\tau=-\int\hat{\phi}^{T}D^{T}\hat{\phi}^{{\dagger}T}d\tau$
$\displaystyle=-\int\hat{\phi}_{\rho}D_{\sigma\rho}\hat{\phi}^{{\dagger}}_{\sigma}=-\int
D\hat{\phi}\hat{\phi}^{{\dagger}}d\tau$
$\displaystyle=\int\\{\hat{\phi}^{{\dagger}}D\hat{\phi}\\}_{n}d\tau-\int\langle\mathrm{vac}|D\hat{\phi}\hat{\phi}^{\dagger}|\mathrm{vac}\rangle
d\tau$ (36) $\displaystyle=h_{n}-D_{p}^{q}\langle
0;\tilde{N}|a_{q}a^{p}|0;\tilde{N}\rangle,\quad p,q\in\mbox{ PES, NES},$ (37)
$\displaystyle h_{n}$ $\displaystyle=D_{p}^{q}\\{a^{p}_{q}\\}_{n},$ (38)

which is to be compared to the usual second-quantized Dirac operator

$\displaystyle\tilde{h}$
$\displaystyle=\int\hat{\phi}^{{\dagger}}D\hat{\phi}d\tau=\int\\{\hat{\phi}^{{\dagger}}D\hat{\phi}\\}_{n}d\tau+\int\langle\mathrm{vac}|\hat{\phi}^{\dagger}D\hat{\phi}|\mathrm{vac}\rangle
d\tau$ (39) $\displaystyle=h_{n}+D_{p}^{q}\langle
0;\tilde{N}|a^{p}a_{q}|0;\tilde{N}\rangle,\quad p,q\in\mbox{ PES, NES}.$ (40)

Since charge conjugation is an inherent symmetry, $\tilde{h}$ (39)/(40) and
$h^{C}$ (36)/(37) should be averaged with an equal weight, leading to

$\displaystyle h$
$\displaystyle=\frac{1}{2}(\tilde{h}+h^{C})=\int\\{\hat{\phi}^{{\dagger}}D\hat{\phi}\\}_{n}d\tau+\frac{1}{2}\int\langle\mathrm{vac}|[\hat{\phi}^{\dagger},D\hat{\phi}]|\mathrm{vac}\rangle
d\tau$ (41) $\displaystyle=h_{n}+C_{n1},$ (42) $\displaystyle C_{n1}$
$\displaystyle=D_{p}^{q}\acontraction[0.5ex]{}{a^{p}}{}{a_{q}}a^{p}a_{q},\quad\acontraction[0.5ex]{}{a^{p}}{}{a_{q}}a^{p}a_{q}=\frac{1}{2}\langle
0;\tilde{N}|a^{p},a_{q}|0;\tilde{N}\rangle=-\frac{1}{2}\delta_{pq}{\rm{sgn}}(\epsilon_{p}),$
(43)

where $\acontraction[0.5ex]{}{a^{p}}{}{a_{q}}a^{p}a_{q}$ is just the CCC (8)
and arises here naturally from the averaging process.

The two-body operator can be calculated in the same way, viz.,

$\displaystyle G^{C}$
$\displaystyle=\frac{1}{2}\int\int\hat{\phi}^{C{\dagger}}(1)\hat{\phi}^{C{\dagger}}(2)V(1,2)\hat{\phi}^{C}(2)\hat{\phi}^{C}(1)d\tau_{1}d\tau_{2}$
$\displaystyle=\frac{1}{2}\int\int\hat{\phi}^{T}(1)\hat{\phi}^{T}(2)[\mathbf{C}_{0}^{\dagger}(1)\mathbf{C}_{0}^{\dagger}(2)V(1,2)\mathbf{C}_{0}(2)\mathbf{C}_{0}(1)]\hat{\phi}^{{\dagger}T}(2)\hat{\phi}^{{\dagger}T}(1)$
$\displaystyle=\frac{1}{2}\int\int\hat{\phi}^{T}(1)\hat{\phi}^{T}(2)V^{*}(1,2)\hat{\phi}^{{\dagger}T}(2)\hat{\phi}^{{\dagger}T}(1)$
$\displaystyle=\frac{1}{2}\int\int\hat{\phi}_{\rho}(1)\hat{\phi}_{\sigma}(2)[V(1,2)]_{\delta\rho,\gamma\sigma}\hat{\phi}^{\dagger}_{\gamma}(2)\hat{\phi}^{\dagger}_{\delta}(1)$
$\displaystyle=\frac{1}{2}\int\int
V(1,2)\hat{\phi}(1)\hat{\phi}(2)\hat{\phi}^{\dagger}(2)\hat{\phi}^{\dagger}(1)$
(44) $\displaystyle=\frac{1}{2}g_{pr}^{qs}a_{q}a_{s}a^{r}a^{p}$ (45)
$\displaystyle=\frac{1}{2}g_{pr}^{qs}\\{a^{pr}_{qs}\\}_{n}-\bar{g}_{pr}^{qs}\langle
0;\tilde{N}|a_{s}a^{r}|0;\tilde{N}\rangle\\{a^{p}_{q}\\}_{n}+\frac{1}{2}\bar{g}_{pr}^{qs}\langle
0;\tilde{N}|a_{q}a^{p}|0;\tilde{N}\rangle\langle
0;\tilde{N}|a_{s}a^{r}|0;\tilde{N}\rangle,$ (46)

which should be compared to the usual two-body operator

$\displaystyle\tilde{G}$
$\displaystyle=\frac{1}{2}\int\int\hat{\phi}^{{\dagger}}(1)\hat{\phi}^{{\dagger}}(2)V(1,2)\hat{\phi}(2)\hat{\phi}(1)d\tau_{1}d\tau_{2}$
$\displaystyle=\frac{1}{2}g_{pr}^{qs}a^{pr}_{qs}$ (47)
$\displaystyle=\frac{1}{2}g_{pr}^{qs}\\{a^{pr}_{qs}\\}_{n}+\bar{g}_{pr}^{qs}\langle
0;\tilde{N}|a^{r}a_{s}|0;\tilde{N}\rangle\\{a^{p}_{q}\\}_{n}+\frac{1}{2}\bar{g}_{pr}^{qs}\langle
0;\tilde{N}|a^{p}a_{q}|0;\tilde{N}\rangle\langle
0;\tilde{N}|a^{r}a_{s}|0;\tilde{N}\rangle.$ (48)

The average of $\tilde{G}$ and $G^{C}$ leads to

$\displaystyle G$ $\displaystyle=\frac{1}{2}(\tilde{G}+G^{C})=G_{n}+C_{n2},$
(49) $\displaystyle G_{n}$
$\displaystyle=\frac{1}{2}g_{pr}^{qs}\\{a^{pr}_{qs}\\}_{n}+Q_{p}^{q}\\{a^{p}_{q}\\}_{n},\quad
Q_{p}^{q}=\bar{g}_{pr}^{qs}\acontraction[0.5ex]{}{a^{r}}{}{a_{s}}a^{r}a_{s},\quad\quad\acontraction[0.5ex]{}{a^{r}}{}{a_{s}}a^{r}a_{s}=\frac{1}{2}\langle
0;\tilde{N}|a^{r},a_{s}|0;\tilde{N}\rangle,$ (50) $\displaystyle C_{n2}$
$\displaystyle=\frac{1}{2}\bar{g}_{pr}^{qs}[\langle
0;\tilde{N}|a^{p}a_{q}|0;\tilde{N}\rangle\langle
0;\tilde{N}|a^{r}a_{s}|0;\tilde{N}\rangle+\langle
0;\tilde{N}|a_{p}a^{q}|0;\tilde{N}\rangle\langle
0;\tilde{N}|a_{s}a^{r}|0;\tilde{N}\rangle].$ (51)

It can be seen that the sum of $h_{n}$ (38) and $G_{n}$ (50) is just
$H_{n}^{\mathrm{QED(MO)}}$ (4) along with the CCC. In particular, in sharp
contrast to the $Q$-potential in Eq. (50), the second term of $G^{C}$ (46) or
that of $\tilde{G}$ (48) is divergent and cannot be rgularized/renormalized.
As such, the averaging of the Hamiltonian $\tilde{h}+\tilde{G}$ for electrons
and the Hamiltonian $h^{C}+G^{C}$ for positrons (moving in the same external
field $\phi_{\mathrm{ext}}(\bm{r})$ as electrons) is a must rather than merely
a formal step. Note in passing that the constant $C_{n2}$ (51) is different
from the second term of Eq. (7) (which is more symmetric due to the direct use
of the CCC). However, the difference does doe not matter at all, for such
constants will be renormalized away.

The above manipulation reveals that the genuine QED Hamiltonian
$H_{n}^{\mathrm{QED(MO-CCC)}}$ (4), especially the Q-potential in Eq. (50) [or
Eq. (6)/(13)], does stem from the symmetric treatment of electrons and
positrons, as dictated by charge conjugation symmetry of relativistic
fermions. In particular, Eqs. (37), (40), (46), and (48) all follow strictly
the ACR (9).

## References

  * Inoue  _et al._ [2024] N. Inoue, Y. Watanabe, and H. Nakano, J. Chem. Phys. 160, 187102 (2024).
  * Liu [2024] W. Liu, J. Chem. Phys. 160, 187101 (2024).
  * Inoue  _et al._ [2023] N. Inoue, Y. Watanabe, and H. Nakano, J. Chem. Phys. 159, 054105 (2023).
  * Liu and Lindgren [2013] W. Liu and I. Lindgren, J. Chem. Phys. 139, 014108 (2013), (E)144, 049901 (2016).
  * Liu [2014] W. Liu, Phys. Rep. 537, 59 (2014).
  * Liu [2020] W. Liu, J. Chem. Phys. 152, 180901 (2020).
  * Liu [2012] W. Liu, Phys. Chem. Chem. Phys. 14, 35 (2012).
  * Sapirstein and Cheng [2003] J. Sapirstein and K. Cheng, Phys. Rev. A 68, 042111 (2003).
  * Dyall [2012] K. G. Dyall, Chemi. Phys. 395, 35 (2012).
  * Liu [2017] W. Liu, With-pair relativistic Hamiltonians, in  _Handbook of Relativistic Quantum Chemistry_ , edited by W. Liu (Springer-Verlag, Berlin, 2017) pp. 345–373.
  * Schwinger [1951] J. Schwinger, Phys. Rev. 82, 664 (1951).
  * Liu [2015] W. Liu, Int. J. Quantum Chem. 115, 631 (2015), (E)116, 971 (2016).
  * Kutzelnigg [2012] W. Kutzelnigg, Chem. Phys. 395, 16 (2012).


#  Ab-initio Symmetric Quasi-Classical Approach to Investigate Molecular Tully
Models

Braden M. Weight  Department of Physics and Astronomy, University of
Rochester, Rochester, NY 14627, U.S.A.  Arkajit Mandal  Department of
Chemistry, University of Rochester, Rochester, NY 14627, U.S.A.  Pengfei Huo
<EMAIL_ADDRESS>Department of Chemistry, University of Rochester,
Rochester, NY 14627, U.S.A.  The Institute of Optics, Hajim School of
Engineering, University of Rochester, Rochester, New York, 14627

###### Abstract

We perform on-the-fly non-adiabatic molecular dynamics simulations using the
symmetrical quasi-classical (SQC) approach with the recently suggested
molecular Tully models: ethylene and fulvene. We attempt to provide benchmarks
of the SQC methods using both the square and the triangle windowing schemes as
well as the recently proposed electronic zero-point-energy correction scheme
(so-called the $\gamma$ correction). We use the quasi-diabatic propagation
scheme to directly interface the diabatic SQC methods with adiabatic
electronic structure calculations. Our results showcase the drastic
improvement of the accuracy by using the trajectory-adjusted
$\gamma$-corrections, which outperform the widely used trajectory surface
hopping method with decoherence corrections. These calculations provide useful
and non-trivial tests to systematically investigate the numerical performance
of various diabatic quantum dynamics approaches, going beyond simple diabatic
model systems that have been used as the major workhorse in the quantum
dynamics field. At the same time, these available benchmark studies will also
likely foster the development of new quantum dynamics approaches based on
these techniques.

##  I Introduction

Simulating on-the-fly non-adiabatic quantum dynamics in molecular systems
remains a central challenge in modern theoretical chemistry, despite the
impressive progress made in the past several decades[1, 2, 3, 4, 5, 6, 7, 8,
9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 22, 24, 25, 26, 27,
28]. The two main components for performing an on-the-fly quantum dynamics
simulation are (i) obtaining accurate electronic structure information and
(ii) using it to propagate the coupled motion of nuclear and electronic
degrees of freedom (DOF) in an efficient manner[29]. Mixed quantum classical
(MQC) approaches such as the fewest-switches surface hopping[1] (FSSH) and the
mean-field Ehrenfest[30] (MFE) approach, which uses the outputs of electronic
structure methods to evolve the electronic subsystem quantum mechanically and
nuclear DOFs classically, has remained popular for simulating on-the-fly
quantum dynamics. However, the inherent mixed-quantum classical approximation
in these approaches can lead to the break-down of detail balance[31], the
artificial creation of electronic coherence[18] or incorrect chemical
kinetics[18].

In response to these deficiencies, a wide range of non-adiabatic dynamics
approaches have been developed in the diabatic representation, some of which
include partial linearized density matrix[7, 32] (PLDM), symmetrical quasi-
classical[11, 12] (SQC), state-dependent ring polymer molecular dynamics[15,
13, 23], quantum-classical path integral (QCPI) approach[33, 34, 35, 36], and
the quantum classical Liouville equation (QCLE) dynamics.[37, 14] In
particular, the recently developed $\gamma$-SQC has been shown[27] to provide
impressively accurate non-adiabatic photo-dissociation quantum dynamics with
coupled Morse potentials through the adjusted zero point energy parameter of
the mapping variables, thus appearing to be a promising method to simulate on-
the-fly quantum dynamics of complex molecular systems. Testing these
approaches with simple model systems becomes the major workhorse in the
quantum dynamics field. What has been largely missing, on the other hand, are
the calculations that go beyond simple diabatic models. However, reformulating
these approaches from diabatic to adiabatic representation often requires non-
trivial theoretical tasks and introduces new numerical complications.

In our recent works, we have developed the quasi-diabatic (QD) propagation
scheme[38, 39, 40, 41, 42] as a general strategy to seamlessly combine a
diabatic quantum dynamics approach with the adiabatic outputs of an electronic
structure method. The QD propagation scheme uses the adiabatic states with a
reference nuclear geometry (the so-called “crude adiabatic” states) as the
local diabatic states during a short-time propagation and continuously updates
the QD basis at each consecutive nuclear propagation step. In this propagation
scheme, one does not construct a global diabatic representation but uses a
sequence of local diabatic representations for each short-time segment to
propagate quantum dynamics. Note that the quasi-diabatic propagation scheme
[38, 39, 40, 41, 42] should not be confused with the approximate diabatic
representation which are also referred often as the QD representation in the
literature.[43, 44, 45]

In this work we use the QD propagation scheme to combine the $\gamma$-SQC
approach with adiabatic outputs of the complete active space self consistent
field (CAS-SCF) approach to perform on-the-fly non-adiabatic quantum dynamics
in molecular systems. We directly simulate photo-excited non-adiabatic
dynamics in two molecular systems, ethylene and fulvene. These molecular
systems have been recently proposed[46] to be analogues of the Tully curve-
crossing models,[1] I and III, respectively. The original Tully models[1]
explore nuances in excited state dynamics that are ubiquitous in ‘real’
molecular systems to varying degrees of complexity, while only involving one
degree of freedom, and have been extensively used to benchmark quantum
dynamics approaches.[6, 47, 7, 11, 48, 49] The molecular analogues of the
Tully models, on the other hand, capture the basic physics of the original
Tully models while concurrently showcasing complex dynamical features due to
the coupled motion between the electronic and multiple nuclear DOF. These
molecular systems serve as robust benchmarks offering complex non-adiabatic
dynamics beyond one-dimensional, overly simplified model systems and are
representative of typical molecular systems.

Our numerical results demonstrate that the zero-point energy (ZPE) corrected
SQC ($\gamma$-SQC ) improves the population dynamics compared to the original
SQC approach, benchmarked against to the accurate but expensive ab-initio
multiple spawning (AIMS) approach. Our numerical results also show that
$\gamma$-SQC can outperform the state-of-the-art decoherence-corrected surface
hopping (dTSH) approach. Overall, our results demonstrate the accuracy and
applicability of the $\gamma$-SQC approach for ab-initio on-the-fly simulation
enabled by the QD propagation scheme, opening up future opportunities for
simulating on-the-fly quantum dynamics of complex molecular systems.

##  II Theory

The SQC approach[11, 12] uses symmetrical window functions to sample
electronic DOF at initial time and provides an estimate of the reduced density
matrix at later times. It relies on the Meyer-Miller-Stock-Thoss (MMST)
mapping Hamiltonian approach, which transforms the electronic degrees of
freedom onto an effective set of singly-excited and fictitious classical
harmonic oscillators. It has been shown to provide accurate non-adiabatic
dynamics in a wide range of model systems[11, 50, 51, 52]. Recently, the
original SQC method was also combined with the QD propagation scheme[41] for
direct quantum dynamics simulation seamlessly using adiabatic electronic
structure outputs.[42] Here, we briefly discuss the essential idea of the
mapping Hamiltonian, the $\gamma$-SQC approach and the QD propagation scheme.

###  II.1 Mapping Hamiltonian Formalism

The Meyer-Miller-Stock-Thoss (MMST) mapping representation [53, 54, 55]
transforms the discrete electronic DOFs onto an effective set of fictitious,
singly-excited classical harmonic oscillators, thus mapping the electronic
non-adiabatic dynamics onto these oscillators’ phase space motion.

The total molecular Hamiltonian in the diabatic representation is expressed as
follows

$\hat{H}=\hat{T}+\sum_{ij}V_{ij}(\hat{\bf R})|i\rangle\langle j|,$ (1)

where $V_{ij}(\hat{\bf R})=\langle i|\hat{V}(\hat{\bf r},\hat{\bf
R})|j\rangle$ are the matrix elements of the electronic Hamiltonian in the
diabatic basis $\\{|i\rangle\\}$. Using the Meyer-Miller-Stock-Thoss[56, 54,
55] mapping representation, the discrete electronic states are transformed
into continuous phase-space variables

$|i\rangle\langle j|\rightarrow{\hat{a}}_{i}^{\dagger}{\hat{a}}_{j},$ (2)

where ${\hat{a}}^{\dagger}_{i}=({\hat{q}}_{i}-i{\hat{p}}_{i})/\sqrt{2}$ and
${\hat{a}}_{j}=({\hat{q}}_{j}+i{\hat{p}}_{j})/\sqrt{2}$. With this
transformation, the molecular Hamiltonian in Eq. 1 is transformed into the
following MMST mapping Hamiltonian

$\hat{H}_{\mathrm{m}}=\hat{T}+{\frac{1}{2}}\sum_{ij}V_{ij}(\hat{R})\left(\hat{p}_{i}\hat{p}_{j}+\hat{q}_{i}\hat{q}_{j}-2\gamma\delta_{ij}\right),$
(3)

where $\gamma=0.5$ is the ZPE for the mapping harmonic oscillators.
Historically, it is recognized as the Langer correction by Meyer and
Miller[56] for the quasi-classical description. Note, until Eq. 3, no
approximations have been made.

In the SQC approach, instead of evolving all DOF quantum mechanically, the
coupled electronic-nuclear dynamics are propagated using the following
classical Hamiltonian[52]

${H}_{\mathrm{m}}={{{\bf
P}^{2}}\over{2M}}+{\frac{1}{2}}\sum_{ij}V_{ij}(R)\left(p_{i}p_{j}+q_{i}q_{j}-2\gamma\delta_{ij}\right),$
(4)

where $\gamma$ is viewed as a parameter[52] which specifies the ZPE of the
mapping oscillators.[57]

Classical trajectories are generated based on Hamilton’s equations of motion

$\displaystyle\dot{q}_{j}$ $\displaystyle=$ $\displaystyle\partial
H_{\mathrm{m}}/\partial p_{j};~{}~{}\dot{p}_{i}=-\partial
H_{\mathrm{m}}/\partial q_{i}$ (5) $\displaystyle\dot{\bf R}$ $\displaystyle=$
$\displaystyle\partial H_{\mathrm{m}}/\partial{\bf P};~{}~{}\dot{\bf
P}=-\partial H_{\mathrm{m}}/\partial{\bf R}={\bf F},$ (6)

with the nuclear force expressed as

${\bf F}=-{1\over 2}\sum_{ij}\nabla V_{ij}({\bf
R})\big{(}p_{i}p_{j}+q_{i}q_{j}-2\gamma\delta_{ij}\big{)}.$ (7)

Overall, the MMST mapping Hamiltonian provides a consistent classical footing
for both electronic and nuclear DOFs, and the non-adiabatic transitions
between electronic states are captured through the classical motion of the
fictitious harmonic oscillators.

###  II.2 Symmetric Quasi-Classical (SQC) Approach

To sample the electronic initial condition and estimate the population, the
SQC approach uses the action-angle variables, $\\{e_{j},\theta_{j}\\}$, which
are related to the canonical mapping variables $\\{p_{j},q_{j}\\}$ through

$e_{j}=\frac{1}{2}\left(p_{j}^{2}+q_{j}^{2}\right);~{}~{}~{}\theta_{j}=-\tan^{-1}\left(\frac{p_{i}}{q_{i}}\right),$
(8)

and the inverse relations are

$q_{j}=\sqrt{2e_{j}}\cos(\theta_{j});~{}~{}p_{j}=-\sqrt{2e_{j}}\sin(\theta_{j}),$
(9)

where $e_{j}$ is a positive-definite action variable introduced by Cotton and
Miller that is directly proportional to the mapping variables’ radius in
action-space,[27] which allows for conceptual simplification (compared to the
$n_{j}=e_{j}-\gamma$ action variable used in previous works[41, 11, 58]) as it
is independent of the ZPE parameter $\gamma$, which will be allowed to be
state-dependent in subsequent sections of this work.

The SQC approach allows for the population of electronic state $|j\rangle$ to
be evaluated as[52]

$\displaystyle\rho_{jj}(t)$ $\displaystyle=\mathrm{Tr}_{\bf
R}\left[\hat{\rho}_{R}|i\rangle\langle i|e^{i\hat{H}t/\hbar}|j\rangle\langle
j|e^{-i\hat{H}t/\hbar}\right]$ (10)
$\displaystyle\approx\frac{1}{\left(2\pi\hbar\right)^{N+M}}\int
d\bm{\tau}\rho_{\mathrm{W}}({\bf P},{\bf
R})W_{i}(\mathbf{e}(0))W_{j}(\mathbf{e}(t)),$

where $\hat{\rho}(0)=|i\rangle\langle i|\otimes\hat{\rho}_{R}$ is the initial
density operator, $\rho_{\mathrm{W}}({\bf P},{\bf R})$ is the Wigner density
of $\hat{\rho}_{R}$ operator that contains $\mathcal{N}$ nuclear DOFs,
$\mathbf{e}=\\{e_{1},e_{2},...,e_{\mathcal{F}}\\}$ is the positive-definite
action variable vector for $\mathcal{F}$ electronic states,
$W_{i}(\mathbf{e})=\delta(e_{i}-(1+\gamma))\prod_{i\neq
j}\delta(e_{j}-\gamma)$ are the Wigner transformed action variables,[59] and
$d\bm{\tau}\equiv d{\bf P}\cdot d{\bf R}\cdot d\mathbf{e}\cdot d\bm{\theta}$.

For practical reasons, the above delta functions are artificially broadened
using two well-explored distribution functions (i.e., square and triangle)
that can be used to bin the resulting electronic action variables in action-
space, depicted for any two-state projection in Fig. 1a,b.[52] The square
distribution for an $\mathcal{F}$-state system is defined as[27]

$W_{j}({\bf e})=w_{1}(e_{j})\prod_{j^{\prime}\neq
j}^{\mathcal{F}}w_{0}(e_{j^{\prime}}),$ (11)

where the function $w_{N}(e)$ is expressed as

$w_{N}(e)=\begin{cases}1,&0<e-N<2\ \times\ 0.366\\\
0,&\mathrm{else}\end{cases}$ (12)

and $\gamma=\sqrt{3}/2-1\approx 0.366$ is the optimal width parameter of the
square window.[11, 28]

Figure 1: (a) Symmetric quasi-classical (a) square $\square$ and (b) triangle
$\Delta$ window distributions depicted for a two-state projection involving
state 1 (blue) and 2 (red). The action-space window distributions are depicted
using the positive-definite action variables $\\{e_{k}\\}$, which are shifted
quantities according to the corresponding zero-point energy parameter
$\gamma$.

The triangle window[58, 27] is expressed as

$W_{j}({\bf e})=w_{1}(e_{j})\prod_{j^{\prime}\neq
j}^{\mathcal{F}}w_{0}(e_{j},e_{j^{\prime}}),$ (13)

where

$w_{1}(e)=\begin{cases}(2-e)^{2-\mathcal{F}},&1<e<2\\\
0,&\mathrm{else}\end{cases}$ (14)

and

$w_{0}(e,e^{\prime})=\begin{cases}1,&e^{\prime}<2-e\\\
0,&\mathrm{else},\end{cases}$ (15)

and trajectories are assigned to state $j$ at time $t$ if $e_{j}\geq 1$ and
$e_{j^{\prime}}<1$ for all $j^{\prime}\neq j$. The ZPE parameter in the
triangle scheme is $\gamma=1/3$. The triangle window scheme for a 2-state
system is depicted in Fig. 1b.

The time-dependent population at time $t$ is then calculated by applying the
window function estimator to action variables $\\{e_{j}(t)\\}$ for an ensemble
of trajectories. Starting from the initial diabatic state $|i\rangle$, the
time-dependent population of the states $|j\rangle$ is computed with Eqn. 10.
However, by using the window function estimator, the total population is no
longer properly normalized due to the fraction of trajectories that are
outside of any window region at any given time.[11] Thus, the total population
must be normalized[11] with the following procedure:

${{\rho}_{jj}(t)}/{\sum_{i=1}^{N}{\rho}_{ii}(t)}\rightarrow{\rho}_{jj}(t).$
(16)

This SCQ approach provides a dramatic improvement to Ehrenfest dynamics, even
though they utilize the same equation of motion for the coupled electronic-
nuclear DOFs.[52, 60] The SQC method allows for the elimination of known
issues present in Ehrenfest dynamics, including preserving detailed
balance,[31, 61, 60], and this method has shown to be quite accurate in an
assortment of model systems [50, 51, 58, 62] while only needing a few thousand
trajectories for convergence.[52, 11, 51] As such, the SQC method is well-
suited for use in on-the-fly non-adiabatic simulations of real molecular
systems.

###  II.3 The $\gamma$-Correction Approach

It was recently proposed that the mapping zero-point energy should be chosen
in such a way as to constrain the initial force to be composed purely from the
initially occupied state,[27] which was not previously enforced using a fixed
$\gamma$ in Eq. 7. This new scheme has shown to provide a significant
improvement for photo-dissociation problems with coupled Morse potentials[27]
and has been combined with the kinematic momentum approach[63] to carry out
on-the-fly simulations of the methaniminium cation.[64]

The basic logic of this scheme is to choose an appropriate $\gamma_{j}$ for
each state $|j\rangle$ in a given individual trajectory, such that the initial
population is forced to respect the initial electronic excitation focused onto
a single excited state. If the initial electronic state is $|i\rangle$, then

$\gamma_{j}=e_{j}-\delta_{ji},$ (17)

or equivalently,

$\delta_{ji}=e_{j}-\gamma_{j},$ (18)

where the $\\{e_{j}\\}$ are uniformly sampled, and then the $\gamma_{j}$ are
chosen to satisfy Eq. 18.

These $\gamma_{j}$ will be explicitly used in the EOMs in Eqs. 5-7, and in
particular, the nuclear forces are now

${\bf F}=-{1\over 2}\sum_{ij}\nabla V_{ij}({\bf
R})\big{(}p_{i}p_{j}+q_{i}q_{j}-2\gamma_{j}\delta_{ji}\big{)},$ (19)

ensuring the initial forces (at $t=0$) are simply ${\bf F}=-\nabla V_{ii}(R)$.
Previously, without any adjustments to $\gamma_{k}$, the chosen values for
$\gamma_{k}$ were only dependent on the windowing function itself, i.e.,
$\gamma_{k}=0.366$ for the square Windows and $\gamma_{k}=1/3$ for the
triangle windows. With the above $\gamma$-correction method,[27] each
individual trajectory will has its own state-specific $\gamma_{j}$ for state
$|j\rangle$ that is completely independent of the choice of window function.
The choices of the $\gamma$ parameter for different SQC approaches are
summarized in Table 1.

Table 1: Choice of $\gamma$ for state $|j\rangle$ with initial excitation on state $|i\rangle$ of each SQC method employed in this work. Method |  $\gamma$-correction | $\gamma_{j}$
---|---|---
$\square$-SQC | False | 0.366
$\Delta$-SQC | False | 1/3
$\square$-$\gamma$-SQC | True | $e_{j}-\delta_{ji}$
$\Delta$-$\gamma$-SQC | True | $e_{j}-\delta_{ji}$

Note that reformulating $\gamma$-SQC in the adiabatic representation (such as
the kinematically transformed SQC[63]) has been done recently to perform on-
the-fly simulations.[64] However, formulating a quantum dynamics approach in
the adiabatic representation introduces additional numerical issues as the
molecular Hamiltonian in the adiabatic representation involves first and
second derivative coupling elements that are typically sharply peaked around
avoided crossings and become singular at a conical intersection (CI). In our
previous work,[41] we showed that the kinematically transformed SQC (KT-SQC),
formulated in the adiabatic representation, may require substantially small
time step to converge dynamics due to the presence of sharply peaked first
derivative coupling elements, even though the second derivative coupling
elements are removed through mathematical transformation in this
formalism.[63] To this end, we use the Quasi-Diabatic Propagation Scheme to
directly interface the diabatic $\gamma$-SQC approach with adiabatic
electronic structure information.

###  II.4 Quasi-Diabatic Propagation Scheme

In this work, we combine the $\gamma-$SQC approach formulated in the diabatic
representation with adiabatic outputs of an electronic structure approach
using the QD propagation scheme.[38, 41, 42]

Despite recent theoretical progress,[65, 66, 67, 68, 69] strict diabatic
states $\\{|i\rangle,|j\rangle\\}$ are neither uniquely defined nor routinely
available for ‘real’ molecular systems. In contrast, it is convenient to
obtain adiabatic states by solving the following eigenequation

$\hat{V}(\hat{\bf r},{\bf R})|\Phi_{\alpha}({\bf
R})\rangle=E_{\alpha}(\mathbf{R})|\Phi_{\alpha}({\bf R})\rangle,$ (20)

where $\hat{V}(\hat{\bf r};{\bf R})$ is the electronic part of the molecular
Hamiltonian at a nuclear configuration ${\bf R}$, and $|\Phi_{\alpha}({\bf
R})\rangle$ is the adiabatic state which is the eigenstate of
$\hat{V}(\hat{\bf r};{\bf R})$, with the corresponding eigenvalue
$E_{\alpha}(\mathbf{R})$ referred to as the adiabatic potential energy.

Consider a short-time propagation of the nuclear DOFs during
$t\in[t_{0},t_{1}]$, where the nuclear positions evolve from ${\bf R}(t_{0})$
to ${\bf R}(t_{1})$, and the corresponding adiabatic states are
$\\{|\Phi_{\alpha}({\bf R}(t_{0}))\rangle\\}$ and $\\{|\Phi_{\lambda}({\bf
R}(t_{1}))\rangle\\}$. The QD scheme uses the nuclear geometry at time $t_{0}$
as a reference geometry, ${{\bf R}_{0}}\equiv{\bf R}(t_{0})$, and uses the
adiabatic basis $\\{|\Phi_{\alpha}({\bf R}(t_{0}))\rangle\\}$ as the quasi-
diabatic basis during this short-time propagation, such that

$|\Phi_{\alpha}({{\bf R}_{0}})\rangle\equiv|\Phi_{\alpha}({\bf
R}(t_{0}))\rangle,~{}~{}\mathrm{for}~{}t\in[t_{0},t_{1}].$ (21)

With the above QD basis defined independently of $R(t)$ within each
propagation segment (or nuclear time step), $\hat{V}(\hat{\bf r},{\bf R})$ in
the QD basis becomes off-diagonal, while the derivative couplings vanish.
Within this basis, all of the necessary diabatic quantities can be evaluated
and used to propagate quantum dynamics during $t\in[t_{0},t_{1}]$.

The electronic Hamiltonian operator $\hat{V}(\hat{\bf r},{\bf R}(t))$ in the
QD basis is evaluated as

$V_{\alpha\beta}({\bf R}(t))=\langle\Phi_{\alpha}({{\bf R}_{0}})|\hat{V}({\bf
R}(t))|\Phi_{\beta}({{\bf R}_{0}})\rangle.$ (22)

For on-the-fly simulations, this quantity is obtained from a linear
interpolation[70] between $V_{\alpha\beta}({\bf R}(t_{0}))$ and
$V_{\alpha\beta}({\bf R}(t_{1}))$ as follows

$\displaystyle V_{\alpha\beta}({\bf R}(t))$
$\displaystyle=V_{\alpha\beta}({\bf
R}_{0})+\frac{(t-t_{0})}{(t_{1}-t_{0})}\big{[}V_{\alpha\beta}({\bf
R}(t_{1}))-V_{\alpha\beta}({\bf R}_{0})\big{]},$ (23)

where $V_{\alpha\beta}({\bf R}_{0})=\langle\Phi_{\alpha}({{\bf
R}_{0}})|\hat{V}({\bf R}(t_{0}))|\Phi_{\beta}({\bf
R_{0}})\rangle=E_{\alpha}({\bf R}(t_{0}))\delta_{\alpha\beta}$. The matrix
elements $V_{\alpha\beta}({\bf R}(t_{1}))$ are computed as follows

$V_{\alpha\beta}({\bf
R}(t_{1}))=\sum_{\lambda\nu}S_{\alpha\lambda}{V}_{\lambda\nu}({\bf
R}(t_{1}))S^{\dagger}_{\beta\nu},$ (24)

where ${V}_{\lambda\nu}({\bf R}(t_{1}))=\langle\Phi_{\lambda}({\bf
R}(t_{1}))|\hat{V}({\bf R}(t_{1}))|\Phi_{\nu}({\bf
R}(t_{1}))\rangle=E_{\lambda}({\bf R}(t_{1}))\delta_{\lambda\nu}$, and the
overlap matrix between two adiabatic electronic states (at two different
nuclear geometries) are $S_{\alpha\lambda}=\langle\Phi_{\alpha}({{\bf
R}_{0}})|\Phi_{\lambda}({\bf R}(t_{1}))\rangle$ and
$S^{\dagger}_{\beta\nu}=\langle\Phi_{\nu}({{\bf R}}(t_{1}))|\Phi_{\beta}({\bf
R}_{0})\rangle$. These overlap matrices are computed based on the approach
outlined in Ref. 71.

The nuclear gradients $\nabla V_{\alpha\beta}({\bf R}(t_{1}))\equiv\partial
V_{\alpha\beta}({\bf R}(t_{1}))/\partial{\bf R}$ are evaluated as

$\displaystyle\nabla V_{\alpha\beta}({\bf
R}(t_{1}))=\nabla\langle\Phi_{\alpha}({\bf R}_{0})|\hat{V}({\bf
R}(t_{1}))|\Phi_{\beta}({\bf R}_{0})\rangle$
$\displaystyle=\langle\Phi_{\alpha}({\bf R}_{0})|\nabla\hat{V}({\bf
R}(t_{1}))|\Phi_{\beta}({\bf R}_{0})\rangle$ (25)
$\displaystyle=\sum_{\lambda\nu}S_{\alpha\lambda}\langle\Phi_{\lambda}({\bf
R}(t_{1}))|\nabla\hat{V}({\bf R}(t_{1}))|\Phi_{\nu}({\bf R}(t_{1}))\rangle
S^{\dagger}_{\beta\nu}.$

To derive the above expression, we have used the fact that
$\\{|\Phi_{\alpha}({\bf R}_{0})\rangle\\}$ is a diabatic basis during the
$[t_{0},t_{1}]$ propagation, allowing us to move the gradient operator to
bypass $\langle\Phi_{\alpha}({\bf R}_{0})|$. We have also inserted the
resolution of identity $\sum_{\lambda}|\Phi_{\lambda}({\bf
R}(t_{1}))\rangle\langle\Phi_{\lambda}({\bf R}(t_{1}))|=1$ in the second line
of the above equation and assume that the QD basis at nuclear position ${\bf
R}(t_{1})$ is complete. It should be noted that Eq. II.4 includes derivatives
with respect to all possible sources of the nuclear dependence which include
those from the adiabatic potentials and the adiabatic states[39, 40].

During the next short-time propagation segment $t\in[t_{1},t_{2}]$, the QD
scheme adopts a new reference geometry ${{\bf R}^{\prime}_{0}}\equiv{\bf
R}(t_{1})$ and new diabatic basis $|\Phi_{\mu}({{\bf
R}^{\prime}_{0}})\rangle\equiv|\Phi_{\mu}({\bf R}(t_{1}))\rangle$. Between
$[t_{0},t_{1}]$ propagation and $[t_{1},t_{2}]$ propagation segments, all of
these quantities will be transformed from $\\{|\Phi_{\alpha}({\bf
R}_{0})\rangle\\}$ to $\\{|\Phi_{\mu}({\bf R}^{\prime}_{0})\rangle\\}$ basis,
using the relation

$|\Phi_{\mu}({\bf R}(t_{1}))\rangle=\sum_{\alpha}\langle\Phi_{\alpha}({\bf
R}(t_{0}))|\Phi_{\mu}({\bf R}(t_{1}))\rangle|\Phi_{\alpha}({\bf
R}(t_{0}))\rangle.$ (26)

Since the mapping relation between the physical state and the singly excited
oscillator state is $|\Phi_{\mu}({\bf
R}(t_{1}))\rangle=a_{\mu}^{\dagger}|0\rangle=\frac{1}{\sqrt{2}}(\hat{q}_{\mu}+i\hat{p}_{\mu})|0\rangle$
and $|\Phi_{\alpha}({\bf
R}(t_{0}))\rangle=a_{\alpha}^{\dagger}|0\rangle=\frac{1}{\sqrt{2}}(\hat{q}_{\alpha}+i\hat{p}_{\alpha})|0\rangle$,
the relations for the mapping variables associated with two bases are

$\displaystyle\frac{1}{\sqrt{2}}(\hat{q}_{\mu}+i\hat{p}_{\mu})|0\rangle=\sum_{\alpha}$
$\displaystyle\langle\Phi_{\alpha}({\bf R}(t_{0}))|\Phi_{\mu}({\bf
R}(t_{1}))\rangle$ (27)
$\displaystyle\times\frac{1}{\sqrt{2}}(\hat{q}_{\alpha}+i\hat{p}_{\alpha})|0\rangle.$

For molecular systems, one can always find a suitable choice for the basis set
in order to make $\langle\Phi_{\alpha}({\bf R}(t_{0}))|\Phi_{\mu}({\bf
R}(t_{1}))\rangle$ real, which guarantees that the mapping variables are
transformed with the same relations as the bases. Based on this relation in
Eq. 27, we transform the time-dependent mapping variables between the two
consecutive QD bases as follows

$\displaystyle\sum_{\alpha}q_{\alpha}\langle\Phi_{\alpha}({\bf
R}(t_{0}))|\Phi_{\mu}({\bf R}(t_{1}))\rangle\rightarrow q_{\mu}$ (28a)
$\displaystyle\sum_{\alpha}p_{\alpha}\langle\Phi_{\alpha}({\bf
R}(t_{0}))|\Phi_{\mu}({\bf R}(t_{1}))\rangle\rightarrow p_{\mu}$ (28b)

When performing the transformation in Eq. 26 and Eq. 28, the eigenvectors
maintain their mutual orthogonality subject to a very small error when they
are expressed in terms of the previous basis due to the incompleteness of the
basis.[72, 8] Nevertheless, the orthogonality remains to be well satisfied
among $\\{|\Phi_{\alpha}({\bf R}(t_{0}))\rangle\\}$ or
$\\{|\Phi_{\lambda}({\bf R}(t_{1}))\rangle\\}$. This small numerical error
generated from each step can, however, accumulate over many steps and cause a
significant error at longer times, leading to non-unitary dynamics.[72, 8]
This potential issue can be easily resolved by using orthonormalization
procedure among the vectors of the overlap matrix $\langle\Phi_{\alpha}({\bf
R}(t_{0}))|\Phi_{\mu}({\bf R}(t_{1}))\rangle$, as been done in our previous
work[39] for simulating photoinduced charge transfer dynamics. Here, we
perform the Löwdin orthogonalization procedure[73] as commonly used in the the
local diabatization approach[72] to ensure this.

As the nuclear geometry closely follows the reference geometry throughout the
propagation, the QD basis forms a convenient and compact basis. Note that, in
principle, one needs infinite crude adiabatic states $\\{|\Phi_{\alpha}({\bf
R}_{0})\rangle\\}$ to represent the time-dependent electronic wavefunction,
because the electronic wavefunction could change rapidly with the motion of
the nuclei, and the crude adiabatic basis is only convenient when the
reference geometry ${\bf R}_{0}$ is close to the nuclear geometry ${\bf R}$.
By dynamically updating the basis in the QD scheme, the time-dependent
electronic wavefunction is expanded with the “moving crude adiabatic
basis”[74] that explores the most relevant and important parts of the Hilbert
space, thus requiring few states for quantum dynamics propagation.

Thus, the QD representation provides several unique advantages over the strict
diabatic or adiabatic representation for quantum dynamics propagation. On one
hand, the QD basis is constructed from the crude adiabatic basis, which can be
easily obtained from any commonly used electronic structure calculation. On
the other hand, the diabatic nature of the QD basis makes derivative couplings
explicitly vanish and allows using any diabatic dynamics approaches to perform
on-the-fly propagation. Further, the QD scheme ensures a stable propagation of
the quantum dynamics compared to directly solving it in the adiabatic
representation. This is due to the fact that directly solving electronic
dynamics in the adiabatic state requires the non-adiabatic coupling
$\langle\Phi_{\beta}({\bf R}(t))|{\partial\over{\partial t}}\Phi_{\alpha}({\bf
R}(t))\rangle={\bf d}_{\beta\alpha}({\bf R})\cdot\dot{\bf R}$, which might
exhibit highly peaked values and cause large numerical errors[75, 20] when
using a linear interpolation scheme.[76] The QD scheme explicitly alleviates
this difficulty by using the well behaved transformation matrix elements
$\langle\Phi_{\beta}({\bf R}(t_{1}))|\Phi_{\alpha}({\bf R}(t_{2}))\rangle$
instead of $\langle\Phi_{\beta}({\bf R}(t))|{\partial\over{\partial
t}}\Phi_{\alpha}({\bf R}(t))\rangle$.

Note that the SQC approach has been derived in the adiabatic
representation,[63] which contains the derivative of the derivative coupling
$\nabla{\bf d}_{\alpha\beta}({\bf R})$ in the equation of motion. With
kinematic momentum transform, it explicitly eliminates the presence of
$\nabla{\bf d}_{\alpha\beta}({\bf R})$ term in the nuclear force (instead of
ignoring it), and the nuclear EOM of the KM-SQC approach is equivalent to the
nuclear forces in QD-SQC.[41] On the other hand, the KM-SQC EOMs still
explicitly contain $\mathbf{d}_{\beta\alpha}(\mathbf{R})$ (through the
presence of $\langle\Phi_{\beta}({\bf R}(t))|{\partial\over{\partial
t}}\Phi_{\alpha}({\bf R}(t))\rangle$, which could lead to numerical
instabilities when these derivative couplings are highly peaked. This has been
extensively investigated in Ref. 41 using Tully’s avoid crossing model, which
has a very narrow derivative coupling, such that it can drastically change on
a time-scale that is shorter than the nuclear time step $dt$. When using a
large $dt$ in the KM-SQC approach, the nuclear position can encounter
drastically different values of the derivative coupling from one step to
another that allow for a discontinuous spike at a CI or even completely step
over it,[17] resulting in different long-time populations and oscillatory
behavior of errors. Thus, the approaches that explicitly require derivative
couplings (and those which use a simple linear interpolation scheme for
obtaining them, as we have previously implemented for the KM scheme) either
encounter numerical challenges or start to accumulate numerical errors.[17]
The QD scheme, on the other hand, provides more accurate results even when
using a relatively larger time step $dt$, simply because the QD schemes only
requires the well-behaved transformation matrix elements $\langle\Phi_{1}({\bf
R}(t_{0}))|\Phi_{2}({\bf R}(t_{1}))\rangle$ instead of the highly peaked
derivative coupling ${\bf d}_{12}({\bf R})$. That being said, there may be
good alternative approaches to achieve the same attractive features for
dynamics propagation, such as those recently developed norm-preserving
interpolation schemes.[17, 20] The QD scheme is perhaps, still the most
straightforward one that allows robust dynamical propagation and enables a
seamless interface between the diabatic quantum dynamics approach (such as
SQC) and adiabatic electronic structure calculations.

###  II.5 Computational Details

Non-adiabatic molecular dynamics simulations based on the QD-$\gamma$-SQC
approach are performed using an in-house modified version[42] of the SHARC
non-adiabatic molecular dynamics code, interfaced to the MOLPRO electronic
structure package.[77, 78] On-the-fly electronic structure calculations are
performed at the level of complete active space self-consistent field (CASSCF)
approach, with 3SA-CASSCF(2,2)/6-31G* and 2SA-CASSCF(6,6)/6-31G* level of
theory for ethylene and fulvene, respectively.[46] The CAS self-consistent
calculation is performed over three lowest adiabatic states for ethylene and
over two lowest adiabatic states for fulvene, whereas the quantum dynamics for
both molecules are only confined in the $\\{\mathrm{S}_{0},\mathrm{S}_{1}\\}$
subspace. All of the energies and gradients are computed at this level of
electronic structure theory. The default accuracy for both the nuclear
gradients and non-adiabatic vectors is $10^{-7}$ a.u.; when this criterion is
not satisfied, a maximum of 900 additional wavefunction optimization
iterations are used to make sure the convergence of $10^{-4}$ a.u. is reached.
All of the electronic structure calculations performed during our quantum
dynamics simulations converge successfully under the above criteria.

Figure 2: Adiabatic potentials for (a) Tully’s model I (a single avoided
crossing) and (b) Tully’s model III (an extended region of coupling with a
reflective barrier). The molecular structures of the ab-initio Tully models
are depicted in the insets. Along a single QD-$\gamma$-SQC trajectory, the
population of the $\mathrm{S}_{1}$ (blue) state and the energies of the
$\mathrm{S}_{0}$ (black) and $\mathrm{S}_{1}$ (red) states as functions of
time for (c) ethylene and (d) fulvene are presented.

The initial Wigner distribution is sampled from the ground vibrational state
$\nu=0$ on the ground electronic state $\mathrm{S}_{0}$, where the normal mode
frequencies (in the harmonic approximation) are calculated based on the
approach outlined in Ref. 79, 80, as implemented in the SHARC package. The
normal mode frequencies are computed at the level of MP2/6-31++G** with the
MOLPRO package, with the optimized structure obtained at the same level of
electronic structure theory for the ground state $|S_{0}({\bf R})\rangle$. In
particular, the nuclear density $\rho_{W}(\tilde{\bf{R}},\tilde{\bf{P}})$ in
terms of the molecular normal-mode frequencies $\\{\tilde{\omega}_{k}\\}$ and
phase space variables $\\{\tilde{\bf{R}},\tilde{\bf{P}}\\}$ is given as[81]

$\rho_{\mathrm{W}}(\tilde{\bf{R}},\tilde{\bf{P}})\propto\prod_{k=1}^{\mathcal{N}}\exp[{-\tanh(\frac{\beta\hbar\tilde{\omega}_{k}}{2})(\frac{m\tilde{\omega}_{k}}{\hbar}\tilde{{R}}_{k}^{2}+\frac{1}{m_{k}\tilde{\omega}_{k}\hbar}\tilde{P}_{k}^{2})}].$
(29)

The initial distribution $\\{{\bf{R}},{\bf{P}}\\}$ is then obtained by
transforming $\\{\tilde{\bf{R}},\tilde{\bf{P}}\\}$ from the normal mode
representation to the primitive coordinates using the unitary transformation
that diagonalizes the Hessian matrix.

A total number of 500 trajectories are used in the QD-$\gamma$-SQC simulations
to achieve converged population dynamics. A rough convergence of the
population dynamics can be already achieved within 50-100 trajectories for
both molecules. The nuclear time step used in the QD-$\gamma$-SQC is $dt=0.1$
fs, with 200 electronic time steps for the mapping variables’ integration
during each nuclear time step. The electronic structure calculations are
performed only at the nuclear time step.

The overlap matrix of CAS wavefunctions between two successive nuclear time
steps is calculated by using the approach outlined in Ref. 71, as implemented
in the program wfoverlap. The random phases generated from electronic
structure calculations for the eigenfunctions are carefully calculated and
accumulated. To ensure the orthonormalization among the vectors of the overlap
matrix $\langle\Phi_{\alpha}({\bf R}(t_{0}))|\Phi_{\mu}({\bf
R}(t_{1}))\rangle$, we perform the Löwdin orthogonalization. [72] All of the
above routines were used as implemented in the SHARC program[82, 83].

In this work, we chose two molecular systems (i) ethylene and (ii) fulvene,
that were previously investigated as the ab-initio analogies to Tully’s curve
crossing models,[46] as illustrated in Fig. 2a,b. In panel (c), the time-
dependent adiabatic potential energies of the ethylene molecule along a single
nuclear trajectory for $\mathrm{S}_{0}$ and $\mathrm{S}_{1}$ states are
presented. In the top sub-panel, the excited state population (blue) computed
with the QD-$\gamma$-SQC approach is plotted, which shows an oscillation of
the population in the the avoided crossing region that eventually relaxes down
to the ground state. In panel (d), the time-dependent adiabatic potential
energies of the fulvene molecule are shown and exhibit many instances of
potential energy barriers within regions of near degeneracy between states,
leading to reflection given a sufficiently low momentum, similar to the
Tully’s model III (which is described canonically by a single spatial
coordinate).[1] The previously encountered avoided crossing is then visited
again where the population transfers back into the excited state. This occurs
twice (within the allotted time-frame) in the single trajectory presented in
panel (d), which gives rise to the important features found in the population
dynamics discussed in more detail later. These two molecules provide rich
physics of two vastly different levels of dynamical complexity. It is
important to note that the features which the Tully models present are with
respect to a single nuclear coordinate, but in the ab initio molecular models
the analogy is represented through the time-dependent adiabtic potential
energies, demonstrated in Fig. 2c,d.

##  III Results and Discussion

Figure 3: Population dynamics of the $\mathrm{S}_{1}$ state in ethylene, using
the square (a,b) and triangle (c,d) windowing schemes. Panels (a,c) utilize
the original SQC method with fixed zero-point energy (ZPE) parameter $\gamma$,
while (b,c) were computed using the trajectory-adjusted ZPE parameter
$\gamma$.

Fig. 3 presents the adiabatic population dynamics of ethylene upon photo-
excitation to the $S_{1}$ state, obtained from various SQC approaches (black
solid lines), compared to the ab-initio multiple spawning (AIMS) (thick green
lines), an approximate Gaussian wavepacket-based non-adiabatic method, which
is used as the benchmark result. The decoherence-corrected trajectory surface
hopping (dTSH) approach (blue lines) is also presented for comparison. Both
AIMS and dTSH results are directly adapted from Ref. 46

Fig. 3a presents the results obtained from the original[12, 11] SQC approach
using the square window scheme, where the mapping ZPE $\gamma=0.366$ is kept
as a constant for all states and trajectories. Fig. 3b presents the results
obtained from the SQC approach using the triangle window scheme[58] with a
fixed $\gamma=1/3$. Both approaches provide reasonably accurate non-adiabatic
dynamics compared to the AIMS results, as well as to the dTSH simulations. In
particular, the $\square$-SQC method seems to show an increased relaxation
time compared to AIMS, whereas the triangle window scheme presented in panel
(b) is more accurate for simulating the Ethylene non-adiabatic dynamics,
compared to the square window scheme presented in panel (a). This trend is in
an agreement with the empirical results of the recent numerical tests of both
window schemes with a wide range of diabatic models,[58, 84] especially for
models with weak non-adiabatic coupling.[58]

Fig. 3c,d presents the ZPE-corrected QD-$\gamma$-SQC dynamics, obtained with
the square window (panel c) and the triangle window (panel d). For the square
window scheme, we find that the $\gamma$-SQC approach (panel c) provides much
better agreement with the AIMS benchmark compared to the original SQC method
(panel a). For the triangle window scheme, we find very similar short-time
relaxation curves for the $\Delta$-SQC (panel b) method compared to
$\Delta$-$\gamma$-SQC (panel d). For the particular case of the ethylene
photo-dissociation dynamics, the ZPE correction does not further improve the
results when the triangle window is used, similar to the recent work that
utilized the kinematic momentum formulation of SQC.[64] With both $\Delta$-SQC
(panel b) and $\Delta$-$\gamma$-SQC (panel d), we see a near quantitative
agreement with AIMS, even slightly out-performing the commonly used dTSH
(blue).

Fig. 4 presents the non-adiabatic photo-relaxation dynamics of fulvene, which
has recently been proposed as a molecular example of Tully’s model III.46 In
fulvene, there exists a so-called slanted CI whereby the wavepacket becomes
reflected and re-interacts with this CI at nearly periodic times later ($\sim
20$-fs intervals), leading to the break down of the mixed quantum-classical
methodologies due to wavepacket bifurcation as well as the added effects of
encircling the CI. These non-adiabtic methods assume a single Gaussian
wavepacket basis for describing the nuclear motion – higher-order modes or
multiple Gaussian wavepackets are needed to fundamentally capture these
effects. For the single SQC trajectory presented in Fig. 2d, one can see that
the population (blue solid line) resides in the ground state while traversing
the CI but jumps back to the excited state after the interaction. Eventually,
the CI will lead to the permanent relaxation to the ground state.

Fig. 4a presents the results obtained from $\square$-SQC, and Fig. 4b presents
the results from $\Delta$-SQC. Both are providing accurate dynamics at the
short time compared to AIMS, including the small shoulder in $S_{1}$
population at $t\sim 10$ fs as well as the subsequent plateaus in
${\mathrm{S}_{1}}$ population. dTSH results (blue) also provides a reasonable
description for the short-time dynamics. However, at a longer time for $t>10$
fs, both SQC and dTSH deviate from the AIMS results, although the SQC approach
(both the square and the triangle windows) outperform dTSH. Even using
different decoherence schemes in dTSH seems unable to provide further
improvement, as shown in Ref. 85.

Fig. 4c,d presents the SQC dynamics with the $\gamma$ correction. Contrary to
the ethylene case, the SQC dynamics is significantly improved upon the
$\gamma$-correction for both window schemes. In particular, the
$\square$-$\gamma$-SQC method (panel c) provides the most quantitative
accuracy among all of the SQC-related methods, capturing both the the initial
inversion of population as well as the population plateau around 20 fs and
finally the tail plateau around 40 fs, slightly outperforming the
$\Delta$-$\gamma$-SQC method (panel d). Both $\gamma$-SQC methods greatly
outperform the dTSH simulation, as there is a severe underestimation of the
population transfer in the dTSH method.[46]

Figure 4: $\mathrm{S}_{1}$ Population dynamics for fulvene using the square
(a,b) and triangle (c,d) windowing schemes. Panels (a,c) utilize the original
SQC method with fixed zero-point energy (ZPE) parameter $\gamma$, while (b,c)
were computed using the trajectory-adjusted ZPE parameter $\gamma$.

##  IV Conclusions

In this work, we use the quasi-diabatic propagation scheme to directly
interface the diabatic symmetric quasi-classical (SQC) approach with the
electronic zero-point energy correction (the $\gamma$ correction)[27] and the
CASSCF on-the-fly electronic structure calculations to propagate ab-initio
non-adiabatic dynamics. We have performed simulations for two recently
suggested molecular models, ethylene and fulvene, that are closely related to
the well-known Tully’s simple curve crossing models. We have shown that the
$\gamma$-SQC method based on the trajectory-adjusted electronic zero-point
energy in classical Meyer-Miller vibronic dynamics provide very accurate non-
adiabatic population dynamics when comparing to ab initio multiple spawning
(AIMS), and even outperforms the widely-used adiabatic decoherence-adjusted
surface hopping (dTSH) method. Specifically, for the fulvene molecule (which
is an molecular analog of Tully’s model III, we found that the
$\gamma$-correction significantly improved the accuracy of the original SQC
approach, for both the square and triangle window schemes. These calculations
provide useful and non-trivial tests to systematically investigate the
numerical performance of various diabatic quantum dynamics approaches, going
beyond the simple diabatic model systems that have historically been the major
workhorse in the quantum dynamics field. At the same time, these available
benchmark studies will also likely foster the development of new quantum
dynamics approaches.

##  V Acknowledgments

This work was supported by the National Science Foundation CAREER Award under
the Grant No. CHE-1845747, as well as by a Cottrell Scholar award (a program
by Research Corporation for Science Advancement). B. M. W. appreciates the
support from the graduate fellowship of the Department of Physics and
Astronomy at the University of Rochester. A.M. appreciates the support from
the ACS Chemical Computing Group Excellence Award for Graduate Students. P.H.
appreciates the support from the ACS COMP OpenEye Outstanding Junior Faculty
Award. Computing resources were provided by the Center for Integrated Research
Computing (CIRC) at the University of Rochester.

## Availability of Data

The data that support the findings of this study are available from the
corresponding author upon a reasonable request.

## References

  * Tully [1990] J. C. Tully, J. Chem. Phys. 93, 1061 (1990).
  * M. Ben-Nun and Martinez [2000] J. Q. M. Ben-Nun and T. J. Martinez, J. Phys. Chem. A 104, 5161 (2000).
  * Micha [1999] D. A. Micha, J. Phys. Chem. A 103, 7562 (1999).
  * Bonella and Coker [2001] S. Bonella and D. F. Coker, J. Chem. Phys. 114, 7778 (2001).
  * Worth  _et al._ [2004] G. A. Worth, M. A. Robba, and I. Burghardt, Faraday Discuss. 127, 307 (2004).
  * Ananth  _et al._ [2007] N. Ananth, C. Venkataraman, and W. H. Miller, J. Chem. Phys. 127, 084114 (2007).
  * Huo and Coker [2011] P. Huo and D. F. Coker, J. Chem. Phys. 135, 201101 (2011).
  * Plasser  _et al._ [2012] F. Plasser, G. Granucci, J. Pittner, M. Barbatti, M. Persico, and H. Lischka, J. Chem. Phys. 137, 22A514 (2012).
  * Saita and Shalashilin [2012] K. Saita and D. V. Shalashilin, J. Chem. Phys. 137, 22A506 (2012).
  * Makhov  _et al._ [2014] D. V. Makhov, W. J. Glover, T. J. Martinez, and D. V. Shalashilin, J. Chem. Phys. 141, 054110 (2014).
  * Cotton and Miller [2013a] S. J. Cotton and W. H. Miller, J. Chem. Phys. 139, 234112 (2013a).
  * Cotton and Miller [2013b] S. J. Cotton and W. H. Miller, J. Phys. Chem. A 117, 7190 (2013b).
  * Ananth [2013] N. Ananth, J. Chem. Phys. 139, 124102 (2013).
  * Hsieh and Kapral [2013] C.-Y. Hsieh and R. Kapral, J. Chem. Phys. 138, 134110 (2013).
  * Richardson and Thoss [2013] J. O. Richardson and M. Thoss, J. Chem. Phys. 139, 031102 (2013).
  * Menzeleev  _et al._ [2014] A. R. Menzeleev, F. Bell, and T. F. Miller, J. Chem. Phys. 140, 064103 (2014).
  * Meek and Levine [2014a] G. A. Meek and B. G. Levine, J. Phys. Chem. Lett. 5, 2351 (2014a).
  * Subotnik  _et al._ [2016] J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, and N. Bellonzi, Annu. Rev. Phys. Chem. 67, 387 (2016).
  * Nelson  _et al._ [2016] T. Nelson, A. Naumov, S. Fernandez-Alberti, and S. Tretiak, Chemical Physics 481, 84 (2016).
  * Jain  _et al._ [2016] A. Jain, E. Alguire, and J. E. Subotnik, J. Chem. Theory Comput. 12, 5256 (2016).
  * Walters and Makri [2016] P. L. Walters and N. Makri, J. Chem. Phys. 144, 044108 (2016).
  * Pal  _et al._ [2016] S. Pal, D. J. Trivedi, A. V. Akimov, B. Aradi, T. Frauenheim, and O. V. Prezhdo, J. Chem. Theory Comput. 12, 1436 (2016).
  * Chowdhury and Huo [2017] S. N. Chowdhury and P. Huo, J. Chem. Phys. 147, 214109 (2017).
  * Curchod and Martinez [2018] B. F. E. Curchod and T. J. Martinez, Chem. Rev. 118, 3305 (2018).
  * Crespo-Otero and Barbatti [2018] R. Crespo-Otero and M. Barbatti, Chem. Rev. 118, 7026 (2018).
  * Mulvihill  _et al._ [2019] E. Mulvihill, A. Schubert, X. Sun, B. D. Dunietz, and E. Geva, J. Chem. Phys. 150, 034101 (2019).
  * Cotton and Miller [2019a] S. J. Cotton and W. H. Miller, J. Chem. Phys. 150, 194110 (2019a).
  * Runeson and Richardson [2020] J. E. Runeson and J. O. Richardson, J. Chem. Phys. 152, 084110 (2020).
  * Tully [2012] J. C. Tully, J. Chem. Phys. 137, 22A301 (2012).
  * Käb [2002] G. Käb, Phys. Rev. E 66, 046117 (2002).
  * Parandekar and Tully [2006] P. V. Parandekar and J. C. Tully, J. Chem. Theory Comput. 2, 229 (2006).
  * Huo and Coker [2012] P. Huo and D. F. Coker, Mol. Phys. 110, 1035 (2012).
  * Lambert and Makri [2012a] R. Lambert and N. Makri, J. Chem. Phys. 137, 22A552 (2012a).
  * Lambert and Makri [2012b] R. Lambert and N. Makri, J. Chem. Phys. 137, 22A553 (2012b).
  * Banerjee and Makri [2013] T. Banerjee and N. Makri, J. Phys. Chem. B 117, 13357 (2013).
  * Makri [2015] N. Makri, Int. J. Quantum Chem. 115, 1209 (2015).
  * Hsieh and Kapral [2012] C.-Y. Hsieh and R. Kapral, J. Chem. Phys. 137, 22A507 (2012).
  * Mandal  _et al._ [2018a] A. Mandal, S. S. Yamijala, and P. Huo, J. Chem. Theory Comput. 14, 1828 (2018a).
  * Mandal  _et al._ [2018b] A. Mandal, F. A. Shakib, and P. Huo, J. Chem. Phys. 148, 244102 (2018b).
  * Mandal  _et al._ [2019] A. Mandal, J. S. Sandoval C., F. A. Shakib, and P. Huo, J. Phys. Chem. A 123, 2470 (2019).
  * Sandoval  _et al._ [2018] J. S. Sandoval, A. Mandal, and P. Huo, J. Chem. Phys. 149, 044115 (2018).
  * Zhou  _et al._ [2019] W. Zhou, A. Mandal, and P. Huo, J. Phys. Chem. Lett. 10, 7062 (2019).
  * Guan  _et al._ [2019] Y. Guan, H. Guo, and D. R. Yarkony, J. Chem. Phys. 150, 214101 (2019).
  * Zhu and Yarkony [2012] X. Zhu and D. R. Yarkony, J. Chem. Phys. 137, 22A511 (2012).
  * Wang  _et al._ [2019] Y. Wang, C. Xie, H. Guo, and D. R. Yarkony, J. Phys. Chem. A 123, 5231 (2019).
  * Ibele and Curchod [2020] L. M. Ibele and B. F. E. Curchod, Phys. Chem. Chem. Phys. 22, 15183 (2020).
  * Duke and Ananth [2015] J. R. Duke and N. Ananth, J. Phys. Chem. Lett. 6, 4219 (2015).
  * Shakib and Huo [2017] F. A. Shakib and P. Huo, J. Phys. Chem. Lett. 8, 3073 (2017).
  * Gao  _et al._ [2020] X. Gao, M. A. C. Saller, Y. Liu, A. Kelly, J. O. Richardson, and E. Geva, J. Chem. Theory Comput. 16, 2883 (2020).
  * Cotton  _et al._ [2014] S. J. Cotton, K. Igumenshchev, and W. H. Miller, J. Chem. Phys. 141, 084104 (2014).
  * Cotton and Miller [2016a] S. J. Cotton and W. H. Miller, J. Chem. Theory Comput. 12, 983 (2016a).
  * Miller and Cotton [2016a] W. H. Miller and S. J. Cotton, Faraday Discuss. 195, 9 (2016a).
  * Meyera and Miller [1979] H. Meyera and W. H. Miller, J. Chem. Phys. 70, 3214 (1979).
  * Stock and Thoss [1997] G. Stock and M. Thoss, Phys. Rev. Lett. 78, 578 (1997).
  * Thoss and Stock [1999] M. Thoss and G. Stock, Phys. Rev. A 59, 64 (1999).
  * Meyer and Miller [1979] H. Meyer and W. H. Miller, J. Chem. Phys. 70, 3214 (1979).
  * Müller and Stock [1999] U. Müller and G. Stock, J. Chem. Phys. 111, 77 (1999).
  * Cotton and Miller [2016b] S. J. Cotton and W. H. Miller, J. Chem. Phys. 145, 144108 (2016b).
  * Miller and Cotton [2016b] W. H. Miller and S. J. Cotton, J. Chem. Phys. 145, 081101 (2016b).
  * Bellonzi  _et al._ [2016] N. Bellonzi, A. Jain, and J. E. Subotnik, J. Chem. Phys. 144, 154110 (2016).
  * Miller and Cotton [2015] W. H. Miller and S. J. Cotton, J. Chem. Phys. 142, 131103 (2015).
  * Tao [2014] G. Tao, J. Phys. Chem. C 118, 17299 (2014).
  * Cotton  _et al._ [2017] S. J. Cotton, R. Liang, and W. H. Miller, J. Chem. Phys. 147, 064112 (2017).
  * Hu  _et al._ [2021] D. Hu, Y. Xie, J. Peng, and Z. Lan, J. Chem. Theory Comput. 0, 10.1021/acs.jctc.0c01249 (2021).
  * Voorhis  _et al._ [2010] T. V. Voorhis, T. Kowalczyk, B. Kaduk, L.-P. Wang, C.-L. Cheng, and Q. Wu, Annu. Rev. Phys. Chem. 61, 149 (2010).
  * Subotnik  _et al._ [2015] J. E. Subotnik, E. C. Alguire, Q. Ou, B. R. Landry, and S. Fatehi, Acc. Chem. Res. 48, 1340 (2015).
  * Kubas  _et al._ [2014] A. Kubas, F. Hoffmann, A. Heck, H. Oberhofer, M. Elstner, and J. Blumberger, J. Chem. Phys. 140, 104105 (2014).
  * Zeng  _et al._ [2012] X. Zeng, X. Hu, and W. Yang, J. Chem. Theory Comput. 8, 4960 (2012).
  * Sirjoosingh and Hammes-Schiffer [2011] A. Sirjoosingh and S. Hammes-Schiffer, J. Chem. Theory Comput. 7, 2831 (2011).
  * Webster  _et al._ [1991] F. Webster, P. J. Rossky, and R. A. Friesner, Comput. Phys. Commun. 63, 494 (1991).
  * Plasser  _et al._ [2016] F. Plasser, M. Ruckenbauer, S. Mai, M. Oppel, P. Marquetand, and L. Gonzàlez, J. Chem. Theory Comput. 12, 1207 (2016).
  * Granucci  _et al._ [2001] G. Granucci, M. Persico, and A. Toniolo, J. Chem. Phys. 114, 10608 (2001).
  * Löwdin [1950] P. O. Löwdin, J. Chem. Phys. 18, 365 (1950).
  * Joubert-Doriol and Izmaylov [2018] L. Joubert-Doriol and A. F. Izmaylov, J. Chem. Phys. 148, 114102 (2018).
  * Meek and Levine [2014b] G. A. Meek and B. G. Levine, J. Phys. Chem. Lett. 5, 2351 (2014b).
  * Hammes-Schiffer and Tully [1994] S. Hammes-Schiffer and J. C. Tully, J. Chem. Phys. 101, 4657 (1994).
  * Mai  _et al._ [2019] S. Mai, M. Richter, M. Heindl, M. F. S. J. Menger, A. Atkins, M. Ruckenbauer, F. Plasser, L. M. Ibele, S. Kropf, M. Oppel, P. Marquetand, and L. Gonz’alez, “SHARC2.1: Surface Hopping Including Arbitrary Couplings — Program Package for Non-Adiabatic Dynamics,” (2019), published: sharc-md.org.
  * Werner  _et al._ [2012] H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, and M. Schütz, WIREs Comput Mol Sci 2, 242 (2012).
  * Dahl and Springborg [1988] J. P. Dahl and M. Springborg, J. Chem. Phys. 88, 4535 (1988).
  * Schinke [1995] R. Schinke,  _Photodissociation Dynamics: Spectroscopy and Fragmentation of Small Polyatomic Molecules_ (Cambridge University Press, 1995).
  * Tannor [2007] D. Tannor, University Science books: Mill Valley, U.S.A (2007).
  * Mai  _et al._ [2018] S. Mai, P. Marquetand, and L. Gonzalez, WIREs Comput Mol Sci. , e1370 (2018).
  * Mai  _et al._ [2015] S. Mai, P. Marquetand, and L. Gonzalez, Int. J. Quantum Chem. 115, 1215 (2015).
  * Cotton and Miller [2019b] S. J. Cotton and W. H. Miller, J. Chem. Phys. 150, 104101 (2019b).
  * Vindel-Zandbergen  _et al._ [2021] P. Vindel-Zandbergen, L. M. Ibele, J.-K. Ha, S. K. Min, B. F. E. Curchod, and N. T. Maitra, J.Chem. Theory Comput. (2021).