Nonadiabatic molecular dynamics is essential for simulating photoinduced processes, but its application to large molecular systems is often limited by the cost of repeatedly evaluating excited-state energies, forces, and nonadiabatic couplings. Here, we present a GPU-accelerated implementation of trajectory surface hopping within PYSEQM (https://github.com/lanl/PYSEQM/), a PyTorch-based semiempirical quantum chemistry framework. The module combines semiempirical model Hamiltonians with configuration-interaction singles excited states and introduces the nonadiabatic coupling vectors and time-derivative couplings required for surface-hopping dynamics. The implementation further includes treatments of trivial crossings and electronic decoherence, providing a complete framework for propagating coupled electronic and nuclear dynamics. To improve GPU utilization, independent trajectories are propagated concurrently in batches. Benchmark calculations show that this strategy substantially increases aggregate simulation throughput as the batch size is increased for small molecules. We demonstrate the applicability of the method to large molecular systems by simulating the excited-state population dynamics of methyl-and tert-butyl substituted molecular nanostar dendrimers, which contain dense manifolds of electronically excited states. These calculations show that semiempirical surface-hopping dynamics can be extended to molecular systems that remain challenging for conventional ab-initio approaches. In addition, the differentiable PyTorch formulation provides a data-science route for optimizing semiempirical Hamiltonian parameters against reference excited-state properties and, ultimately, dynamical observables. PYSEQM therefore provides a platform for high-throughput nonadiabatic molecular dynamics and for the development of data-driven semiempirical models for photophysics and photochemistry.

We present a multitask machine-learning strategy for improving the prediction of molecular dipole moments by simultaneously training on quantum dipole magnitudes and inexpensive Mulliken atomic charges. With dipole magnitudes as the primary target and assuming only scalar dipole values are available without vector components, we examine whether incorporating lower-quality labels that do not quantitatively reproduce the target property can still enhance model accuracy. Mulliken charges were chosen intentionally as an auxiliary task since they lack quantitative accuracy yet encode qualitative physical information about charge distribution. Our results show that including Mulliken charges with a small weight in the loss function yields up to a 30% improvement in the dipole prediction accuracy. This multitask approach enables the model to learn a more physically grounded representation of charge distributions, thereby improving both the accuracy and the consistency of dipole magnitude predictions. These findings highlight that even auxiliary data of limited quantitative reliability can provide valuable qualitative physical insights, ultimately strengthening the predictive power of machine-learning models for molecular properties.

We present several finite-temperature recursive Fermi-operator expansion schemes based on the second-order spectral projection (SP2) method. Our approach builds on a previous observation that the electronic structure problem, as formulated through a recursive SP2 expansion, can be mapped onto the architecture of a deep neural network. Using this perspective, we generalize SP2 to finite electronic temperatures and construct machine learning models to determine optimized expansion coefficients. These coefficients are trained for a specified chemical potential and electronic temperature and are not available in closed analytical form. However, by employing an appropriate affine rescaling strategy to the Hamiltonian matrix, we eliminate the need to retrain the model during a simulation if the temperature and chemical potential change. Our approach avoids explicit diagonalization and relies solely on highly optimized matrix-matrix multiplication kernels. Compared to state-of-the-art diagonalization, we achieve an order-of-magnitude speedup in the single-particle finite-temperature density matrix calculation for small and moderately sized matrices on modern GPUs and dense matrix multiply units.

This Letter introduces excited-state molecular dynamics in PYSEQM, a GPU-accelerated semiempirical quantum chemistry engine implemented in PyTorch. The new module enables Born-Oppenheimer molecular dynamics (BOMD) using configuration-interaction singles and random phase approximation for excited states, allowing long trajectories and large statistical ensembles to be simulated efficiently on a single GPU. We also implement an extended Lagrangian excited-state BOMD (XL-ESMD) scheme that propagates auxiliary electronic variables, enabling relaxed ground and excited-state convergence thresholds without compromising energy conservation. The excited-state BOMD implementation scales smoothly from small chromophores to a nearly 900-atom dendrimer (taking 6.5 s per MD step). PYSEQM also supports batched execution, allowing many geometries or trajectories to be evaluated in a single GPU launch, substantially increasing throughput and making ensemble-based protocols routine. As a demonstration, we compute absorption, emission, and infrared spectra from trajectories propagated on the ground and first excited states. The XL-ESMD scheme yields identical spectra at significantly lower computational cost, establishing the role of extended Lagrangian based dynamics for efficient excited-state BOMD simulations. Beyond raw performance, PYSEQM's PyTorch foundation provides automatic differentiation for forces, efficient GPU batching, and seamless interfacing with machine learning models. These capabilities position PYSEQM as a practical platform for machine learning-augmented excited-state dynamics and lay the foundation for future data-driven nonadiabatic excited-state dynamics modeling of ultrafast spectroscopic probes.

We present an extension of the shadow extended Lagrangian Born-Oppenheimer molecular dynamics method to excited state molecular dynamics (ESMD) in the context of ΔSCF Kohn-Sham density functional theory, with demonstrations performed using self-consistent charge density functional tight binding (SCC-DFTB) theory. In this shadow ESMD approach, the approximate iterative solution to the exact potential in conventional ESMD is replaced by an exact single-step solution to an approximate shadow excited-state potential. The energy functional that defines this shadow excited-state potential as a stationary (non-aufbau) solution is obtained from a linearization about an approximate excited state density, which would become a progressively worse approximation as the dynamics ensue if it were static. To avoid this, we propagate the approximate excited-state (charge) density as an additional dynamical variable in an extended Lagrangian approach. We show that, in addition to offering significant improvement in computational cost relative to direct ESMD, our shadow ESMD method provides enhanced stability and robustness relative to its "exact" counterpart. Our implementation is carried out in the context of SCC-DFTB theory but should be broadly generalizable, both to ab initio electronic structure methods and to other semi-empirical quantum chemistry approaches.

Shadow molecular dynamics provide an efficient and stable atomistic simulation framework for flexible charge models with long-range electrostatic interactions. Shadow molecular dynamics simulations are driven by approximate "shadow" Born-Oppenheimer potentials for which the exact charges and forces are directly accessible without relying on costly (and approximate) iterative solvers. While previous implementations have been limited to atomic monopole charge distributions, we extend this approach to flexible multipole models. We derive detailed expressions for the shadow energy functions, potentials, and force terms, explicitly incorporating monopole-monopole, dipole-monopole, and dipole-dipole interactions. In our formulation, both atomic monopoles and atomic dipoles are treated as extended dynamical variables alongside the propagation of the nuclear degrees of freedom. We demonstrate that introducing the additional dipole degrees of freedom preserves the stability and accuracy previously seen in monopole-only shadow molecular dynamics simulations. In addition, we present a shadow molecular dynamics scheme where the monopole charges are held fixed while the dipoles remain flexible. Our extended shadow dynamics provide a framework for stable, computationally efficient, and versatile molecular dynamics simulations involving long-range interactions between flexible multipoles. This is of particular current interest in combination with machine-learned interatomic potentials, including long-range electrostatic interactions.

Accurate and efficient gradients of molecular energy with respect to nuclear degrees of freedom are essential for geometry optimization and molecular dynamics, including simulations that go beyond the Born-Oppenheimer regime. A common approach involves deriving analytical formulas for new electronic structure methods, which is often conceptually difficult and requires tedious coding. Here, we implement analytical, semi-numerical, and automatic differentiation (AD)-based gradient pathways for semiempirical Hamiltonian models in the PYSEQM software package, leveraging both graphics processing unit (GPU) and central processing unit (CPU) architectures. We further extend these capabilities to excited states calculated using the configuration interaction singles and time-dependent Hartree-Fock ansätze. We benchmark wall time, peak memory usage, and accuracy across three molecular families of varying chemical complexity, including systems of up to a thousand atoms. For ground-state simulations, analytical and AD gradients achieve near-identical GPU runtimes, while semi-numerical gradients are slower on GPU but remain competitive on CPU. For excited states, both analytical and custom AD approaches using implicit differentiation show similar performance and low memory requirements, whereas gradients with full AD are memory-limited. AD gradients match analytical ones in accuracy across all tested systems, aided by a quaternion-based diatomic frame rotation for two-center quantities that ensures smooth energy surfaces. Overall, automatic differentiation emerges as a practical alternative to analytical gradients in semiempirical quantum chemistry, offering high accuracy while allowing seamless integration in AI-driven workflows and popular packages, such as PyTorch and JAX. Our results provide actionable guidance for selecting optimal gradient strategies in large-scale ground- and excited-state molecular dynamics simulations.

Graph-based linear-scaling electronic-structure theory provides a scalable framework for parallel quantum-mechanical molecular dynamics (QMD) simulations by exploiting the nearsightedness of the non-local electronic connectivity in non-metallic systems. When combined with recent shadow molecular dynamics in an extended-Lagrangian formulation, it enables stable long-time simulations of large, chemically active systems. This article introduces the Scalable Ecosystem, Driver, and Analyzer for Complex Chemistry Simulations (SEDACS), which integrates all these advances within a modular, Python-based software package for large-scale QMD simulations driven by external electronic-structure codes. SEDACS provides a tunable, adaptive graph construction in which edges encode the non-local electronic overlap between atoms. This graph is then decomposed into a set of smaller, overlapping subgraphs, where the electronic structure of each of these subgraphs is solved for independently and in parallel using an external electronic-structure code. SEDACS can be coupled to a variety of external electronic-structure solvers with minimal modifications to their software, enabling rapid adoption of the graph-based QMD approach. In this way, SEDACS can greatly extend the capability of existing electronic-structure packages by enabling stable QMD simulations of systems that were previously computationally inaccessible. We demonstrate highly efficient and stable QMD simulations for chemically active systems with tens of thousands of atoms by interfacing SEDACS with an external Fortran-based electronic-structure code based on self-consistent-charge density functional tight-binding theory.

A complete understanding of enzyme mechanisms requires atomistic details of chemical reactions. Quantum-based molecular dynamics simulations (QMD) are a potential source of this information, but tradeoffs between accuracy and computational cost have limited their use. Here, we develop a reactive QMD approach to investigate mechanisms of isocyanide hydratase (ICH) catalysis. In QMD simulations, molecular analogs of ICH active site residues reacted with para-nitrophenol isocyanide, forming a thioimidate intermediate. Analysis of simulated atomic configurational and charge dynamics revealed a pathway where protonation of the isocyanide carbon occurs prior to thioimidate formation. X-ray crystallography and functional assays of ICH mutants suggest this order of events might occur during enzyme catalysis. Mobile protons play essential roles in many enzymes, yet they are difficult to observe experimentally, making the ordering of proton-dependent steps ambiguous in many enzyme mechanisms. The ability to directly simulate reactions relevant to enzyme catalysis involving mobile protons demonstrates the significance of our reactive QMD approach and motivates further biological applications.

A complete understanding of enzyme mechanisms requires atomistic details of chemical reactions. Quantum-based molecular dynamics simulations (QMD) are a potential source of this information, but trade-offs between accuracy and computational cost have limited their use. We previously developed extended Lagrangian Born-Oppenheimer molecular dynamics (XL-BOMD) methods that leverage a negligible compromise in accuracy to substantially decrease the cost of QMD simulations. Here, we develop a reactive QMD approach using the latest XL-BOMD formulation, which enables efficient simulations of highly reactive systems, and use it to investigate mechanisms of intermediate formation in isocyanide hydratase (ICH) catalysis. In QMD simulations, molecular analogs of ICH active site residues reacted with para-nitrophenyl isocyanide, forming a thioimidate. Analysis of simulated atomic configurational and charge dynamics revealed a pathway where protonation of the isocyanide carbon occurs prior to thioimidate formation and suggested a possible role of Asp17 as a proton donor in the early phase of ICH catalysis. To test whether the pathway seen using the reactive QMD approach might be relevant to ICH catalysis, we performed X-ray crystallography and pre-steady-state enzyme kinetics studies of wild-type and D17N mutant ICH. Both the structure and kinetics are sensitive to the D17N mutation in a manner that is consistent with the order of the reaction steps seen in the simulations. Mobile protons play essential roles in many enzymes, yet they are difficult to observe experimentally, making the ordering of proton-dependent steps ambiguous in many enzyme mechanisms. The ability to directly simulate model reactions for the design of experiments that provide information about enzyme mechanisms involving mobile protons demonstrates the significance of our reactive QMD approach and motivates further biological applications.

Graph-based electronic structure theory offers a scalable approach to study large, complex atomistic systems using distributed and hybrid computational platforms. We demonstrate the coupling of graph-based linear scaling electronic structure theory, as implemented in the Scalable Ecosystem, Driver, and Analyzer for Complex Chemistry Simulations (SEDACS), with semiempirical quantum chemistry methods as implemented in the PySEQM code, with Graphics Processing Unit (GPU) acceleration. This powerful combination enables efficient, scalable electronic structure calculations over many nodes, significantly reducing computational cost while naturally harnessing parallelism. Detailed analyses of parallelization efficiency, computational accuracy, and communication overheads are provided, highlighting an order-of-magnitude speedup for systems of up to 10,000 atoms.

We present a framework for atomistic simulations of surface catalysis under electrochemical bias. The framework makes use of extended Lagrangian Born-Oppenheimer quantum-based molecular dynamics (XL-BOMD) simulations, which provide the speed and accuracy required for explicit atomistic treatment of both electrode and electrolyte. Simulations of solvated O_2 near nitrogen-doped graphene (NG) were performed to gain insight into the oxygen reduction reaction (ORR). Different mechanisms were observed, depending on the applied bias. Under high bias ORR occurred by an outer sphere mechanism, without adsorption of O_2 to NG. In this mechanism, electron transfer between the catalyst and the O_2 was mediated by the solvent. Under low bias ORR occurred by an inner sphere mechanism involving adsorption of O_2 to NG, leading to direct electron transfer. Combining quantum accuracy with explicit solvation and bias, XL-BOMD opens a route to predictive, atomistic insight into electrocatalytic processes beyond the reach of traditional methods.

We report the development and implementation of electronic excited-state capabilities for semiempirical quantum chemical methods at both the Configuration Interaction Singles and Time-Dependent Hartree-Fock levels of theory, integrated within the PYSEQM 2.0 software package (https://github.com/lanl/PYSEQM). PYSEQM is a Python-based package designed for efficient and scalable quantum chemical simulations. Leveraging the PyTorch framework enables PYSEQM to benefit from automatic differentiation and GPU acceleration, leading to substantial performance gains in molecular property evaluations. In particular, our implementation enables efficient calculation of excited-state properties for large molecular systems. Benchmarking on systems with up to a thousand atoms demonstrates that excited-state computations can be completed in under a minute on modern GPUs, making this approach particularly suitable for high-throughput screening, real-time feedback in interactive simulations, and large-scale dynamical studies. Additionally, PYSEQM includes a machine learning interface that supports Hamiltonian parameter reoptimization and neural network training. These capabilities open new avenues for data-driven excited-state dynamics simulations, offering a path toward combining quantum chemical rigor with machine learning efficiency. Overall, this work facilitates access to excited-state quantum chemistry for large systems, while laying the foundation for future hybrid quantum-machine-learning approaches in photochemistry, photophysics, and materials discovery.

We present an extended Lagrangian shadow molecular dynamics scheme with an interatomic Born-Oppenheimer potential determined by the relaxed atomic charges of a second-order charge equilibration model. To parametrize the charge equilibration model, we use machine learning with neural networks to determine the environment-dependent electronegativities and chemical hardness parameters for each atom, in addition to the charge-independent energy and force terms. The approximate shadow molecular dynamics potential in combination with the extended Lagrangian formulation improves the numerical stability and reduces the number of Coulomb potential calculations required to evaluate accurate conservative forces. We demonstrate efficient and accurate simulations with excellent long-term stability of the molecular dynamics trajectories. The significance of choosing fixed or environment-dependent electronegativities and chemical hardness parameters is evaluated. Finally, we compute the infrared spectrum of molecules via the dipole autocorrelation function and compare to experiments to highlight the accuracy of the shadow molecular dynamics scheme with a machine learned flexible charge potential.

We introduce a shadow molecular dynamics (MD) approach based on the Atom-Condensed Kohn-Sham second-order (ACKS2) charge-potential equilibration model. In contrast to regular flexible charge models, the ACKS2 model includes both flexible atomic charges and potential fluctuation parameters that allow for physically correct charge fragmentation and improved scaling of the polarizability. Our shadow MD scheme is based on an approximation of the ACKS2's flexible charge-potential energy function, in combination with extended Lagrangian Born-Oppenheimer MD. Utilizing this shadow charge-potential equilibration approach mitigates the costly overhead and stability problems associated with finding well-converged iterative solutions to the charges and potential fluctuations of the ACKS2 model in an MD simulation. Our work provides a robust and versatile framework for efficient, high-fidelity MD simulations of diverse physical phenomena and applications.

With recent advancements in machine learning for interatomic potentials, Python has become the go-to programming language for exploring new ideas. While machine-learning potentials are often developed in Python-based frameworks, existing molecular dynamics software is predominantly written in lower-level languages. This disparity complicates the integration of machine learning potentials into these molecular dynamics libraries. Additionally, machine learning potentials typically focus on local features, often neglecting long-range electrostatics due to computational complexities. This is a key limitation as applications can require long-range electrostatics and even flexible charges to achieve the desired accuracy. Recent charge equilibration models can address these issues, but they require iterative solvers to assign relaxed flexible charges to the atoms. Conventional implementations also demand very tight convergence to achieve long-term stability, further increasing computational cost. In this work, we present a scalable Python implementation of a recently proposed shadow molecular dynamics scheme based on a charge equilibration model, which avoids the convergence problem while maintaining long-term energy stability and accuracy of observable properties. To deliver a functional and user-friendly Python-based library, we implemented an efficient neighbor list algorithm, Particle Mesh Ewald, and traditional Ewald summation techniques, leveraging the GPU-accelerated power of Triton and PyTorch. We integrated these approaches with the Python-based shadow molecular dynamics scheme, enabling fast charge equilibration for scalable machine learning potentials involving systems with hundreds of thousands of atoms.

Co-design across the Exascale Computing Project has been critical for both enabling science applications and bringing disparate communities together. Developing and porting applications to the various high-performance computing architectures on pre-exascale and exascale computers has been quite challenging due to the diversity of hardware features and software stacks. The Co-design Center for Particle Applications (CoPA) has developed and enhanced the Cabana and Parallel, Rapid O(N), and Graph-Based Recursive Electronic Structure Solver (PROGRESS)/Basic Matrix Library (BML) libraries to facilitate the creation of new particle applications, make existing particle applications exascale capable, and allow teams to explore new capabilities. Particle methods from the atomistic, mesoscale, and continuum through cosmological scales have been built with Cabana, along with new possibilities for application coupling. Similarly, the PROGRESS/BML library has enabled quantum particle applications with linear algebra solvers to use advanced hardware. Across these CoPA-developed libraries, the co-design abstraction layer combines performance portability with math library support to facilitate the separation of concerns and directly support science runs.

This paper introduces a novel framework for matrix diagonalization, recasting it as a sequential decision-making problem and applying the power of Decision Transformers (DTs). Our approach determines optimal pivot selection during diagonalization with the Jacobi algorithm, leading to significant speedups compared to the traditional max-element Jacobi method. To bolster robustness, we integrate an epsilon-greedy strategy, enabling success in scenarios where deterministic approaches fail. This work demonstrates the effectiveness of DTs in complex computational tasks and highlights the potential of reimagining mathematical operations through a machine learning lens. Furthermore, we establish the generalizability of our method by using transfer learning to diagonalize matrices of smaller sizes than those trained.

In this paper, we introduce innovative approaches for accelerating the Jacobi method for matrix diagonalization, specifically through the formulation of large matrix diagonalization as a Semi-Markov Decision Process and small matrix diagonalization as a Markov Decision Process. Furthermore, we examine the potential of utilizing scalable architecture between different-sized matrices. During a short training period, our method discovered a significant reduction in the number of steps required for diagonalization and exhibited efficient inference capabilities. Importantly, this approach demonstrated possible scalability to large-sized matrices, indicating its potential for wide-ranging applicability. Upon training completion, we obtain action-state probabilities and transition graphs, which depict transitions between different states. These outputs not only provide insights into the diagonalization process but also pave the way for cost savings pertinent to large-scale matrices. The advancements made in this research enhance the efficacy and scalability of matrix diagonalization, pushing for new possibilities for deployment in practical applications in scientific and engineering domains.

In recent years, a new kind of accelerated hardware has gained popularity in the artificial intelligence (AI) community which enables extremely high-performance tensor contractions in reduced precision for deep neural network calculations. In this article, we exploit Nvidia Tensor cores, a prototypical example of such AI-hardware, to develop a mixed precision approach for computing a dense matrix factorization of the inverse overlap matrix in electronic structure theory, S-1. This factorization of S-1, written as ZZT = S-1, is used to transform the general matrix eigenvalue problem into a standard matrix eigenvalue problem. Here we present a mixed precision iterative refinement algorithm where Z is given recursively using matrix-matrix multiplications and can be computed with high performance on Tensor cores. To understand the performance and accuracy of Tensor cores, comparisons are made to GPU-only implementations in single and double precision. Additionally, we propose a nonparametric stopping criteria which is robust in the face of lower precision floating point operations. The algorithm is particularly useful when we have a good initial guess to Z, for example, from previous time steps in quantum-mechanical molecular dynamics simulations or from a previous iteration in a geometry optimization.