E-SHAKE: Constrained-Dynamics for Diabatic Seam Sampling

• E-SHAKE is a constrained-dynamics method that samples diabatic seams in nuclear geometries, enabling explicit mapping of electronic coupling variations.
• It adapts the classical SHAKE/RATTLE algorithm to enforce the diabatic energy-gap condition via Lagrange multipliers in ab initio molecular dynamics simulations.
• By resolving seam fluctuations and coupling distributions, E-SHAKE challenges the Condon approximation and refines Marcus theory estimates for energy-transfer systems.

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 (Cofer-Shabica et al., 13 Oct 2025).

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 (Cofer-Shabica et al., 13 Oct 2025). 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 (Cofer-Shabica et al., 13 Oct 2025).

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 (Cofer-Shabica et al., 13 Oct 2025).

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,

$\sigma(\mathbf R)=0$0

where $\sigma(\mathbf R)=0$1 is the adiabatic-to-diabatic transformation matrix (Cofer-Shabica et al., 13 Oct 2025).

For excitation and energy transfer, the work employs localized diabatization, especially the BoysOV scheme. The Boys criterion maximizes the separation of dipole moments,

$\sigma(\mathbf R)=0$2

and the BoysOV variant separates occupied and virtual contributions,

$\sigma(\mathbf R)=0$3

Once the diabatic basis is obtained, the two-state Hamiltonian is written as

$\sigma(\mathbf R)=0$4

This representation isolates the donor and acceptor diabatic energies and the diabatic coupling $\sigma(\mathbf R)=0$5 used in Marcus theory (Cofer-Shabica et al., 13 Oct 2025).

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

$\sigma(\mathbf R)=0$6

with $\sigma(\mathbf R)=0$7 the mass matrix. The underlying integration scheme is velocity Verlet,

$\sigma(\mathbf R)=0$8

$\sigma(\mathbf R)=0$9

where $(3N-1)$0 (Cofer-Shabica et al., 13 Oct 2025).

To impose the seam condition $(3N-1)$1, the method introduces Lagrange multipliers $(3N-1)$2 and $(3N-1)$3 and rewrites the updates as

$(3N-1)$4

$(3N-1)$5

with

$(3N-1)$6

$(3N-1)$7

$(3N-1)$8

The constraint is enforced in two stages. First, a position correction solves for $(3N-1)$9 from

$H_\mathrm{DA}$0

The original implementation solves this numerically by bisection with a Newton-Raphson bracketing step. Second, a velocity correction determines $H_\mathrm{DA}$1 from

$H_\mathrm{DA}$2

Here $H_\mathrm{DA}$3 is the partially unconstrained velocity after the position correction (Cofer-Shabica et al., 13 Oct 2025).

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

$H_\mathrm{DA}$4

but regard this as sufficient because the constraint keeps the trajectory on the seam (Cofer-Shabica et al., 13 Oct 2025).

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 $H_\mathrm{DA}$5 surface while enforcing

$H_\mathrm{DA}$6

This operational condition defines the seam-sampling window used in practice (Cofer-Shabica et al., 13 Oct 2025).

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}$7 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 (Cofer-Shabica et al., 13 Oct 2025).

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 (Cofer-Shabica et al., 13 Oct 2025).

When the couplings were recomputed near the $H_\mathrm{DA}$8 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 (Cofer-Shabica et al., 13 Oct 2025). 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}$9 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 (Cofer-Shabica et al., 13 Oct 2025).

The same seam-sampling data also provide evidence for a conical intersection. When the two conditions

$$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.$$0

are simultaneously met, the diabatic Hamiltonian becomes

$$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.$$1

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 (Cofer-Shabica et al., 13 Oct 2025).

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 (Cofer-Shabica et al., 13 Oct 2025).

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

$$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.$$2

Within that decomposition, the dominant fluctuation was traced to a single hydrogen atom in the equatorial site below the donor, labeled $$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.$$3 (Cofer-Shabica et al., 13 Oct 2025).

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 (Cofer-Shabica et al., 13 Oct 2025).

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 $$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.$$4 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.