---
title: 'E-SHAKE: Constrained-Dynamics for Diabatic Seam Sampling'
url: https://www.emergentmind.com/topics/e-shake
type: topic
---

# E-SHAKE: Constrained-Dynamics for Diabatic Seam Sampling

E-SHAKE is a constrained-dynamics method for sampling the nuclear-geometry seam on which two diabatic electronic states are degenerate. Introduced in the context of revisiting Marcus theory and the Condon approximation for the Closs triplet energy-transfer systems, it adapts the classical SHAKE/RATTLE paradigm from geometric constraints to an electronic constraint, typically the diabatic energy-gap condition $E_\mathrm{D}(\mathbf R)=E_\mathrm{A}(\mathbf R)$ or $\sigma(\mathbf R)=0$. In this form, E-SHAKE enables ab initio molecular dynamics directly on a $(3N-1)$-dimensional seam, so that the geometry dependence of the diabatic coupling $H_\mathrm{DA}$ can be sampled explicitly rather than inferred from a single crossing-point geometry [2510.11810].

## 1. Definition and theoretical setting

E-SHAKE is described as an **electronic** version of the classical **SHAKE/RATTLE** constrained-dynamics algorithms used in molecular dynamics [2510.11810]. Its defining operation is the imposition of an electronic constraint on nuclear motion, rather than a conventional holonomic constraint such as a bond length or bond angle. The principal constraint used in the original work is the diabatic degeneracy condition
\[
E_\mathrm{D}(\mathbf R)=E_\mathrm{A}(\mathbf R),
\qquad
\sigma(\mathbf R)=E_\mathrm{D}(\mathbf R)-E_\mathrm{A}(\mathbf R)=0.
\]

For a system with $3N$ nuclear degrees of freedom, this condition defines a seam of dimension $(3N-1)$. The method was introduced because prior strategies for exploring such seams usually target minimum-energy crossing points or rely on biasing potentials, metadynamics, or NEB-type methods, whereas E-SHAKE is designed to constrain ab initio molecular dynamics directly onto the seam and thereby sample a broad range of seam geometries [2510.11810].

The immediate motivation is Marcus theory. The high-temperature Marcus rate expression used in the study is
\[
k_{\mathrm D\to \mathrm A} =
\frac{2\pi}{\hbar}
|H_{\mathrm{DA}}|^2
(4\pi \lambda k_\mathrm B T)^{-1/2}
\exp\!\left[
-\frac{(\lambda+\Delta G^0)^2}{4\lambda k_\mathrm B T}
\right].
\]
In the formulation under discussion, one crucial assumption is the **Condon approximation**, namely that $H_\mathrm{DA}$ is constant and independent of nuclear geometry. E-SHAKE was introduced precisely to test that assumption by sampling $H_\mathrm{DA}$ across the full seam rather than at a single representative point [2510.11810].

## 2. Diabatic representation and seam definition

The method is built on a diabatic description of donor and acceptor states. Localized diabatic states are constructed by rotating adiabatic states into a diabatic basis,
\[
\ket{\Xi_i}=\sum_{j=1}^{N_\mathrm{states}} U_{ji}\ket{\Phi_j},
\]
where $U$ is the adiabatic-to-diabatic transformation matrix [2510.11810].

For excitation and energy transfer, the work employs localized diabatization, especially the **BoysOV** scheme. The Boys criterion maximizes the separation of dipole moments,
\[
f_\mathrm{Boys}(U)=
\sum_{i,j=1}^{N_\mathrm{states}}
\left|
\langle \hat{\mu}\rangle_{\Xi_i}
-
\langle \hat{\mu}\rangle_{\Xi_j}
\right|^2,
\]
and the BoysOV variant separates occupied and virtual contributions,
\[
f_\mathrm{BoysOV}(U)=
\sum_{i,j=1}^{N_\mathrm{states}}
\left|
\langle \hat{\mu}^{\mathrm{occ}}\rangle_{\Xi_i}
-
\langle \hat{\mu}^{\mathrm{occ}}\rangle_{\Xi_j}
\right|^2
+
\left|
\langle \hat{\mu}^{\mathrm{virt}}\rangle_{\Xi_i}
-
\langle \hat{\mu}^{\mathrm{virt}}\rangle_{\Xi_j}
\right|^2.
\]

Once the diabatic basis is obtained, the two-state Hamiltonian is written as
\[
U^\ddagger
\begin{bmatrix}
E_1 & 0\\
0 & E_2
\end{bmatrix}
U
=
\begin{bmatrix}
E_\mathrm D & H_\mathrm{DA}\\
H_\mathrm{DA} & E_\mathrm A
\end{bmatrix}.
\]
This representation isolates the donor and acceptor diabatic energies and the diabatic coupling $H_\mathrm{DA}$ used in Marcus theory [2510.11810].

Within this framework, the seam sampled by E-SHAKE is not defined by adiabatic degeneracy alone. It is defined by **diabatic** degeneracy. That distinction is central, because the method is intended to monitor how the coupling varies across the seam. A plausible implication is that E-SHAKE is best understood not merely as a crossing-point locator, but as a seam-resolved probe of the validity of constant-coupling rate models.

## 3. Constrained-dynamics formulation

The unconstrained nuclear equations of motion are
\[
\ddot{\mathbf R} = -\mathbf M^{-1}\nabla V(\mathbf R),
\]
with $\mathbf M$ the mass matrix. The underlying integration scheme is velocity Verlet,
\[
\mathbf R(t+h)=\mathbf R(t)+h\dot{\mathbf R}(t)+\frac{h^2}{2}\mathbf M^{-1}\mathbf F(t),
\]
\[
\dot{\mathbf R}(t+h)=\dot{\mathbf R}(t)+\frac{h}{2}\mathbf M^{-1}\big[\mathbf F(t)+\mathbf F(t+h)\big],
\]
where $\mathbf F=-\nabla V$ [2510.11810].

To impose the seam condition $\sigma(\mathbf R)=0$, the method introduces Lagrange multipliers $\lambda_{RR}$ and $\lambda_{RV}$ and rewrites the updates as
\[
\mathbf R(t+h)=\mathbf R(t)+h\mathbf Q,
\]
\[
\dot{\mathbf R}(t+h)=\mathbf Q+\frac{1}{2}\mathbf M^{-1}\!\left[h\mathbf F(t+h)+k\mathbf G(t+h)\right],
\]
with
\[
g=h\lambda_{RR},
\qquad
k=h\lambda_{RV},
\]
\[
\mathbf Q=\dot{\mathbf R}(t)+\frac12 \mathbf M^{-1}[h\mathbf F(t)+g\mathbf G(t)],
\]
\[
\mathbf G(t)=-\nabla \sigma(\mathbf R(t)).
\]

The constraint is enforced in two stages. First, a **position correction** solves for $g$ from
\[
\sigma\!\left(\tilde{\mathbf R}(t+h)+g\frac{h}{2}\mathbf M^{-1}\mathbf G(t)\right)=0.
\]
The original implementation solves this numerically by bisection with a Newton-Raphson bracketing step. Second, a **velocity correction** determines $k$ from
\[
k =
2\,
\frac{
\tilde{\dot{\mathbf R}}(t+h)^{\mathrm T}\mathbf G(t+h)
}{
\mathbf G(t+h)^{\mathrm T}\mathbf M^{-1}\mathbf G(t+h)
}.
\]
Here $\tilde{\dot{\mathbf R}}(t+h)$ is the partially unconstrained velocity after the position correction [2510.11810].

The authors also note that the added cost is essentially **two extra diabatic gradient evaluations per time step**. They further state that the diabatic gradients are approximate, using the **strictly diabatic approximation**
\[
\langle \Xi_i|\nabla_R|\Xi_j\rangle = 0,
\]
but regard this as sufficient because the constraint keeps the trajectory on the seam [2510.11810].

## 4. Seam-sampling protocol

In the reported application, the dynamics are initialized near a minimum-energy crossing between the first two triplet surfaces and propagated on the $T_1$ surface while enforcing
\[
|E_\mathrm D(\mathbf R)-E_\mathrm A(\mathbf R)| < 4\times 10^{-5}\ \text{Hartree} \approx 1\ \text{meV}.
\]
This operational condition defines the seam-sampling window used in practice [2510.11810].

The purpose of the resulting constrained trajectory is not simply to remain at one crossing geometry. Rather, the dynamics explore geometries along the diabatic intersection seam, which allows monitoring of how $H_\mathrm{DA}$ varies with energy above the minimum seam point. This is the key methodological distinction between E-SHAKE and approaches that only locate a minimum-energy crossing point [2510.11810].

Within the logic of the study, seam sampling serves as a direct test of whether the Condon approximation is tenable. If the coupling remains tightly distributed around a finite mean value across the seam, the constant-coupling picture remains plausible. If the coupling changes strongly over the seam, then the Condon approximation fails, and the simplest Marcus picture becomes questionable. The method therefore turns seam sampling into a diagnostic for the rate-theory assumptions themselves.

## 5. Application to the Closs triplet energy-transfer systems

The method was introduced to revisit one disagreement between theory and experiment in the naphthalene-bridge-biphenyl and naphthalene-bridge-benzophenone systems studied by Piotrowiak, Miller, and Closs, with particular focus on the molecule **C-13-ae** [2510.11810].

When the couplings were recomputed near the $T_1/T_2$ crossing, most systems behaved similarly, but **C-13-ae** was identified as a strong outlier: its computed coupling was **about two orders of magnitude smaller than C-13-ea** [2510.11810]. For **C-13-ea** and **C-13-ee**, the coupling is described as fairly tightly distributed around a finite mean value across the seam, which is consistent with the Condon approximation. For **C-13-ae**, by contrast, $H_\mathrm{DA}$ is **centered near zero** and varies by roughly **two orders of magnitude** over the seam. This was taken as direct evidence that the coupling depends strongly on nuclear geometry and that the Condon approximation fails for that system [2510.11810].

The same seam-sampling data also provide evidence for a conical intersection. When the two conditions
\[
E_\mathrm D = E_\mathrm A
\qquad \text{and} \qquad
H_\mathrm{DA}\to 0
\]
are simultaneously met, the diabatic Hamiltonian becomes
\[
\begin{bmatrix}
E & 0\\
0 & E
\end{bmatrix},
\]
which the work identifies as the hallmark of a **conical intersection** rather than a generic seam. The distribution of couplings around zero for C-13-ae was interpreted accordingly as evidence for an accessible conical intersection in the seam [2510.11810].

On this basis, the authors predict that **C-13-ae should exhibit a much slower triplet-triplet energy transfer rate than C-13-ea**, and that a slow rate scale likely exists that was not resolved experimentally. The analysis also argues that, although a conical intersection can accelerate downhill relaxation in photoexcited systems, in the thermally activated regime relevant here it can act like a barrier or funnel that requires trajectories to reach a special geometry where coupling becomes nonzero, thereby slowing transfer [2510.11810].

## 6. Structural interpretation and broader significance

The study goes beyond aggregate seam statistics and analyzes the structural origin of the coupling fluctuations through the per-atom decomposition
\[
\frac{dV}{dt}^{(i)} =
\sum_{\mu=x,y,z}
\frac{dV}{dR_{i\mu}}
\frac{dR_{i\mu}}{dt}.
\]
Within that decomposition, the dominant fluctuation was traced to a **single hydrogen atom** in the equatorial site below the donor, labeled $\mathrm{H_3^e}$ [2510.11810].

Because that atom dominates the seam fluctuations, the authors propose that replacing H by D at that site should produce a large isotopic effect for C-13-ae. Their stated expectation is that deuteration would alter seam dynamics and should make the C-13-ae transfer rate faster relative to H, providing an experimental signature consistent with conical-intersection-mediated transfer [2510.11810].

In methodological terms, E-SHAKE establishes a way to treat the seam itself as a sampled object rather than as a single optimized geometry. This shifts the analysis of donor-acceptor transfer away from point estimates of $H_\mathrm{DA}$ and toward seam-resolved distributions. The reported results suggest that the method is most consequential when rate predictions hinge on the assumption that the coupling is effectively constant. In that setting, E-SHAKE functions as a direct test of whether the seam supports a stable finite coupling or instead contains regions where the coupling approaches zero and the Marcus/Condon picture breaks down.

Source: https://www.emergentmind.com/topics/e-shake