Organizations: Department of Chemistry, The University of Hong Kong, Pok Fu Lam, Hong Kong SAR, China · Hong Kong Quantum AI Lab, Pak Shek Kok, Hong Kong SAR, China · MattVerse Limited, Pak Shek Kok, Hong Kong SAR, China · Shenzhen Institute for Advanced Study, University of Electronic Science and Technology of China, Shenzhen, China
Nonadiabatic molecular dynamics provides mechanistic insight into light-driven processes and informs the design of molecules and materials for solar energy conversion, photocatalysis and photo switching. Accurately describing these processes requires a representation that respects electronic symmetry and consistently relates energies to interstate couplings. Here we introduce PEACE, which combines a parity-equivariant latent Hamiltonian with a learned electronic connection. Controlled ablations reveal the complementary roles of symmetry-allowed state mixing and electronic-frame variation in reproducing crossing structures and relaxation dynamics. PEACE closely reproduces excited-state population dynamics from first-principles simulations, while its extension to spin-orbit coupling enables simulations of intersystem crossing. These results demonstrate that a more complete incorporation of the underlying physics into learned electronic representations leads to more accurate predictions of nonadiabatic dynamics.
Figures & tables
Fig. 1: The overall architecture of PEACE framework. a , The atomic species Z and Cartesian coordinates R define a molecular graph G=(V,E) . Nodes V are represented by one-hot encodings of atomic species, while edges E are described by radial basis functions of interatomic distances and spherical harmonics of interatomic directions. b, A shared O (3)-equivariant tensor-product encoder combines message aggregation and equivariant updates. c , Two independently parameterized equivariant self-attention branches feed the Hamiltonian and connection heads. In the selected electronic basis, pi∈{+1,−1} are latent-basis parity labels, not labels fixed to energy-sorted adiabatic states. Same-parity diabatic Hamiltonian ( H ) entries use even scalar 0 e features, and opposite-parity entries use odd pseudoscalar 0 o features. Pseudoscalars construct Hamiltonian entries and their coordinate derivatives enter the physical readout. Connection same or opposite-parity pairs use polar 1 o or axial 1 e channels, respectively. Bμ is antisymmetric in electronic indices and combines directly learned rigid and internal motion components using the implemented regularized projections. d , The adiabatic energies E are eigenvalues of H , while the forces F and smoothed non-adiabatic couplings (NACs) b are calculated through covariant derivatives.
Fig. 2: Couplings, crossing topology and dynamics of CH 2 NH 2 + . Curves compare PEACE with models lacking the connection, the opposite-parity Hamiltonian block or both. a , b , smoothed nonadiabatic coupling (NAC) and raw NAC norms, and phase-insensitive absolute collinearity of the predicted and reference coupling vectors, along an S 0 /S 1 torsional scan and an S 1 /S 2 C–N stretching scan. Undefined directions for vanishing predicted couplings are omitted. c , S 1 /S 2 energy-gap maps in common even and odd rigid-motion-free displacement directions; the reference and each model are centred on their own mirror-plane crossing. The MR-CISD map uses a 15 × 15 grid, while model maps use 41 × 41 grids. White contours mark gaps of 0.02, 0.06 and 0.12 eV. d , Mean electronic populations from trajectories initialized in S 2 . Shading shows pointwise 95% intervals from trajectory-bootstrap resamples. The reference and four model ensembles contain 145, 919, 926, 901 and 960 independently retained trajectories in legend order with maximum energy drift less than 0.25 eV.
Fig. 3: Ensemble electronic populations in four molecules . Solid curves show PEACE and dotted curves show the reference. The maximum energy drift of all trajectories is less than 0.25 eV. a - c C 2 H 4 , C 3 H 6 , and C 4 H 8 over 100 fs. The reference ensemble contains 200 trajectories for each system; PEACE ensembles contain 919, 906 and 920, respectively. Shading is the archived pointwise bootstrap interval. d , CH 2 S over 3 ps, with 954 PEACE trajectories and a pointwise interval standard errors. The reference curves were extracted from the reference 20 .
Supplementary Fig. 1: Learning curves for PEACE and the ablation models. a , Training loss. b , Validation loss used for checkpoint selection. c – e , Validation root-mean-square errors (RMSEs) for energies, Cartesian force components and smoothed nonadiabatic coupling (SNAC) components, respectively. PEACE is compared with models excluding the connection contribution (w/o connection), opposite-parity Hamiltonian blocks (w/o H-odd), or both (w/o both). Each model was trained for 1,000 epochs using 2,500 training and 250 validation geometries. Stars mark the validation-selected checkpoints.
Supplementary Fig. 2: Static test predictions for the CH 2 NH 2 + ablations. Rows a – d correspond to PEACE, w/o connection, w/o H-odd and w/o both, respectively. Columns show predicted versus MR-CISD/aug-cc-pVDZ reference energies, Cartesian force components and SNAC components 4 , 5 , 6 . All four models use validation-selected checkpoints and the same 1,250 held-out test geometries. Insets report mean absolute errors (MAE) and root mean squared errors (RMSE).
Supplementary Fig. 3: Static test predictions of SPaiNN models for CH 2 NH 2 + . SPaiNN 1 models use a , SchNet 12 and b , PaiNN 13 architectures. Columns show predicted versus MR-CISD/aug-cc-pVDZ reference energies, Cartesian force components and SNAC components. Both models use validation-selected checkpoints and the same 1,250 held-out test geometries as Supplementary Fig. 2 . Insets report MAE and RMSE.
Supplementary Fig. 4: CH 2 NH 2 + population dynamics. Reference S 0 , S 1 and S 2 populations are compared with PEACE, SPaiNN–SchNet, SPaiNN–PaiNN 1 and Exciting DeePMD Dyad 14 over 100 fs. The four machine-learning methods use the same 1,000 initial geometries and velocities, with trajectories initiated in S 1 . Independent selection for complete trajectories with a maximum absolute total-energy drift below 0.25 eV retains 947, 418, 758 and 453 model trajectories, respectively. Shading indicates pointwise 95% confidence intervals from trajectory-level bootstrap resampling 16 .
Supplementary Fig. 5: Unseen 770-geometry CH 2 NH 2 + benchmark. Energy, force and SNAC MAE are evaluated on the same 770 geometries from the SchNarc benchmark 7 . Comparisons include PEACE, Exciting DeePMD 14 , SchNarc 7 and SchNet- and PaiNN-based SPaiNN models 1 . Force and SNAC errors are averaged over Cartesian components. Each predicted SNAC vector was sign-aligned independently for each geometry and state pair by choosing the sign that minimized its error relative to the reference vector. The public SchNarc 7 checkpoint has no NAC prediction head, so no SNAC error is reported.
Supplementary Fig. 6: Static predictions for C 2 H 4 . Predicted versus SA(3)-CASSCF(2,2)/cc-pVDZ reference energies, Cartesian force components and SNAC components for a , PEACE; b , SPaiNN–SchNet; and c , SPaiNN–PaiNN. 1 Dashed lines indicate equality. Insets report MAE and RMSE for the test set.
Supplementary Fig. 7: Electronic-state populations and S 1 relaxation dynamics of C 2 H 4 . a , b , SPaiNN populations using SchNet ( a ) and PaiNN ( b ) 1 . Solid lines show model predictions and dotted lines the SA(3)-CASSCF(2,2)/cc-pVDZ reference. c , S 1 populations from the reference (black), PEACE (blue), SPaiNN–SchNet (red) and SPaiNN–PaiNN (yellow). Solid lines show ensemble means; dashed lines show fits to the effective one-step S1→S0 model, PS1(t)=exp(−t/τ) , over 0–100 fs. Fitted decay times are 68.8, 68.5, 59.3 and 69.6 fs, respectively. Shading indicates pointwise 95% confidence intervals from trajectory-level bootstrap resampling 16 . Trajectories start in S 1 ; only complete trajectories with a maximum absolute total-energy drift below 0.25 eV are retained.
Supplementary Fig. 8: Static predictions for C 3 H 6 . Predicted versus SA(3)-CASSCF(2,2)/cc-pVDZ reference energies, Cartesian force components and SNAC components for a , PEACE; b , SPaiNN–SchNet; and c , SPaiNN–PaiNN. 1 Dashed lines indicate equality. Insets report MAE and RMSE for the test set.
Supplementary Fig. 9: Electronic-state populations and S 1 relaxation dynamics of C 3 H 6 . a , b , SPaiNN populations using SchNet ( a ) and PaiNN ( b ) 1 . Solid lines show model predictions and dotted lines the SA(3)-CASSCF(2,2)/cc-pVDZ reference. c , S 1 populations from the reference (black), PEACE (blue), SPaiNN–SchNet (red) and SPaiNN–PaiNN (yellow). Solid lines show ensemble means; dashed lines show fits to the effective one-step S1→S0 model, PS1(t)=exp(−t/τ) , over 0–100 fs. Fitted decay times are 53.3, 51.4, 41.2 and 107.5 fs, respectively. Shading indicates pointwise 95% confidence intervals from trajectory-level bootstrap resampling 16 . Trajectories start in S 1 ; only complete trajectories with a maximum absolute total-energy drift below 0.25 eV are retained.
Supplementary Fig. 10: Static predictions for C 4 H 8 . Predicted versus SA(3)-CASSCF(2,2)/cc-pVDZ reference energies, Cartesian force components and SNAC components for a , PEACE; b , SPaiNN–SchNet; and c , SPaiNN–PaiNN 1 . Dashed lines indicate equality. Insets report MAE and RMSE for the test set.
Supplementary Fig. 11: Electronic-state populations and S 1 relaxation dynamics of C 4 H 8 . a , b , SPaiNN populations using SchNet ( a ) and PaiNN ( b ) 1 . Solid lines show model predictions and dotted lines the SA(3)-CASSCF(2,2)/cc-pVDZ reference. c , S 1 populations from the reference (black), PEACE (blue), SPaiNN–SchNet (red) and SPaiNN–PaiNN (yellow). Solid lines show ensemble means; dashed lines show fits to the effective one-step S1→S0 model, PS1(t)=exp(−t/τ) , over 0–100 fs. Fitted decay times are 66.9, 67.8, 99.6 and 129.6 fs, respectively. Shading indicates pointwise 95% confidence intervals from trajectory-level bootstrap resampling 16 . Trajectories start in S 1 ; only complete trajectories with a maximum absolute total-energy drift below 0.25 eV are retained.
Supplementary Fig. 12: Static predictions for CH 2 S including spin–orbit coupling. Predicted versus independently recomputed OpenMolcas reference values 8 for a , energies; b , Cartesian force components; c , complex spin–orbit couplings (SOCs); and d , SNAC components on the test set. In c , real and imaginary parts are plotted separately, whereas the inset errors use the modulus of the complex difference.
Supplementary Fig. 13: CH 2 S properties along a C–S bond-length scan. PEACE predictions and independent CASSCF(6,5)/def2-SVP OpenMolcas reference calculations 8 are compared for a , electronic energies; b , selected signed SOC components; and c , SNAC norms. Panel b shows Re(SOC) + Im(SOC) for the M s = +1 triplet component. Panel c shows Euclidean norms over all 12 Cartesian components. The 34 plotted geometries span C–S distances of 1.52–2.18 Å.
Supplementary Fig. 14: Per-geometry electronic-evaluation times for PEACE and CASSCF. FP64 PEACE inference at batch size one on a single NVIDIA H100 GPU is compared with SA(3)-CASSCF(2,2)/cc-pVDZ calculations using six CPU workers on an Intel Xeon Platinum 8480C system. Each alkene uses 4,000 common configurations from 20 complete trajectories, with three PEACE repeats after warm-up. Error bars and shading indicate 95% confidence intervals from paired trajectory-block bootstrap resampling 16 . Ratios are computed from the mean wall times and are specific to these hardware and workflow settings.
Supplementary Fig. 15: Numerical tests of symmetry consistency. Relative errors in energies, forces, and smoothed nonadiabatic couplings (NACs) were evaluated in FP64 across 32 molecular geometries subjected to rotations, reflections, and translations. Boxes span the interquartile range, with horizontal lines indicating medians. Whiskers extend to the most extreme observations within 1.5 times the interquartile range; points denote outliers.
Dataset
Quantity
SPaiNN (SchNet)
SPaiNN (PaiNN)
PEACE
C 2 H 4
Energies (eV)
0.026
0.026
0.008
Forces (eV/Å)
0.146
0.113
0.036
NACs (Bohr -1 )
0.096
0.099
0.028
C 3 H 6
Energies (eV)
0.031
0.016
0.013
Forces (eV/Å)
0.149
0.091
0.042
NACs (Bohr -1 )
0.081
0.044
0.016
Supplementary Table 1: Static mean absolute errors across four molecules. Errors compare PEACE with SchNet- and PaiNN-based SPaiNN models 1 . NAC rows report raw derivative-coupling errors in atomic units.
Dataset
Training
Validation
Test
C 2 H 4
6000
500
1000
C 3 H 6
6000
500
1000
C 4 H 8
12000
2000
1000
CH 2 NH 2 +
2500
250
1250
CH 2 S
4000
200
503
Supplementary Table 2: Dataset partitions for training and evaluation.
We describe a neural network architecture and training procedure designed to model electronic ground and excited states of arbitrary molecular systems. By indirectly learning a latent, implicit basis representation of the electronic-state Hamiltonian, the model offers a unified treatment of multiple electronic states, conical intersections, and non-adiabatic couplings. The formalism can be further extended to learn consistent latent representations of additional operators such as transition dipole moments, for example. To demonstrate the general capabilities of our architecture, we train and evaluate networks on two realistic photochemical systems, thymine and azobenzene. The resulting models accurately reproduce energies and oscillator strengths for the ground- and low-lying excited states relevant to the photochemistry of these systems. We highlight the performance of the trained networks by studying critical molecular geometries, including conical intersections and excited state minima. By construction, the proposed framework also recovers the emergence of Berry phase accumulation around conical intersections. By pairing key mathematical structure from quantum chemistry with the representation learning power of transformers, the presented architecture offers a qualitatively new path toward fast and accurate ground- and excited-state simulations.
David Juergens, Martin Stöhr, Andreas E. Hillers-Bendtsen +2
Department of Chemistry and The PULSE Institute, Stanford University, Stanford 4305, California, USA · SLAC National Accelerator Laboratory, Menlo Park, 94025, California, USA
Accurate ab initio molecular dynamics (AIMD) simulations of complex, fluxional chemical systems are severely limited by the high computational scaling of correlated electronic structure methods. To overcome this bottleneck, we present a robust, graph-theoretic molecular fragmentation framework integrated with machine learning to directly model post-Hartree-Fock nuclear forces at coupled cluster accuracy. Bypassing the limitations of automatic differentiation on learned energy surfaces that may struggle with link-atom Jacobians, our approach directly predicts nuclear force vectors. By projecting these vectors onto fragment-fixed principal axes of inertia, we establish co-variant descriptors that naturally preserve rotational, translational, and permutational invariance. The methodology achieves exceptional high parameter efficiency through a vector-valued training protocol that reduces trainable parameters by over an order of magnitude, while an unsupervised mini-batch k-means space tessellation algorithm constructs highly representative training databases using only 10% to 20% of reference configurations. We rigorously validated this framework on the highly fluxional solvated Zundel cation H_{13}O_6^+ ). Our fully machine-learning-predicted AIMD trajectories successfully reproduced complex dynamical signatures and key structural characteristics, including radial distribution functions and the velocity autocorrelation power spectrum. Ultimately, this scalable, systematically improvable framework bridges the gap between high-level correlated wavefunction theories and long-timescale reactive sampling, laying the foundation for advanced, LLM-inspired transfer learning in modern chemical dynamics simulations.
Xiao Zhu, Srinivasan S. Iyengar
Department of Chemistry, Department of Physics, and the Indiana University Quantum Science and Engineering Center (IU-QSEC), Indiana University, 800 E. Kirkwood Ave, Bloomington, IN-47405
Machine-learning interatomic potentials (MLIPs) have enabled molecular dynamics at near ab initio accuracy, yet remain limited to energies and forces by construction, leaving electronic observables such as dipole moments and polarizabilities inaccessible. We introduce DenSNet, a density-first approach to machine-learned electronic structure that learns the Hohenberg--Kohn map from nuclear configurations to the ground-state electron density. Our approach employs an SE(3)-equivariant neural network to predict density coefficients of a flexible atom-centered Gaussian basis, combined with a Δ-learning strategy that uses superposed atomic densities as a prior to accelerate training. A second equivariant network then maps the predicted density to the total energy, providing a unified framework for molecular dynamics and electronic structure. We validate DenSNet on ethanol, ethanethiol, and resorcinol, where infrared spectra from machine-learned trajectories show excellent agreement with experimental gas-phase measurements. To test scalability, we train on polythiophene oligomers with 1--6 monomers and extrapolate to chains of up to 12 monomers, generating stable long-time trajectories whose infrared spectra agree with reference density functional theory calculations. Here, we show that reinstating the electron density as the central learned quantity opens a practical route to transferable prediction of spectroscopic and electronic observables in large-scale molecular simulations.
Mihail Bogojeski, Muhammad R. Hasyim, Leslie Vogt-Maranto +3
BIFOLD and Machine Learning Group, Technische Universität Berlin, Franklinstr. 28/29, 10587 Berlin, Germany · Department of Chemistry, New York University, New York, NY 10003, USA · Department of Artificial Intelligence, Korea University, Seoul 02841, Korea +5