Lithium diffusion in silicon battery anodes is governed by thermally activated jumps between (meta)stable sites separated by significant energy barriers, making such events rare on ab initio molecular dynamics (AIMD) time scales. To overcome this limitation, we establish a multiscale workflow that links AIMD, machine-learned force fields (MLFFs), and Markov state models (MSMs) to bridge atomistic mechanisms to mesoscale diffusion. Focusing on crystalline Li-Si phases, our MLFFs trained on AIMD data, achieve near-DFT accuracy while enabling large-scale molecular dynamics simulations extending to tens of nanoseconds. From these trajectories, we extract converged lithium-jump statistics to construct MSMs that quantitatively reproduce diffusivities with uncertainties an order of magnitude smaller than those obtained from 100 ps AIMD simulations. Demonstrated here for crystalline LixSiy phases, the AIMD → MLFF → MSM workflow provides a transferable route for quantitative transport modeling in amorphous structures, defect-mediated diffusion, and alternative solid-state anodes.

In this work, a new strategy to truncate high-order terms in the many-body expansion (MBE) is proposed. This new approach, which we call a hierarchical many-body expansion (HMBE), is based on a hierarchical partition of the system into multitier fragments and can in principle be applied to any large molecular system. Numerical tests on a series of (H2O)64 structures are presented, demonstrating satisfactory relative energies between the clusters and binding energies of individual clusters compared with full-cluster calculations, with significantly fewer high-order terms computed than conventional MBE. The hierarchical truncation can be augmented by certain many-body terms for fragments at the interface between the partitions (called "Schengen terms") to further improve accuracy. This work establishes the HMBE scheme as a promising framework to model very large systems (e.g., proteins), which are naturally built on a hierarchical structure.

We report a k-point extension of the second-order co-iterative augmented Hessian (CIAH) algorithm, termed k-CIAH, for Pipek-Mezey (PM) localization of Wannier functions (WFs). By exploiting an efficient evaluation of the Hessian-vector product, k-CIAH achieves O(Nk2n3) scaling in both CPU time and memory, matching that of previously reported first-order k-space approaches while improving upon the O(Nk3n3) scaling of Γ-point CIAH, where Nk denotes the number of k-points sampling the first Brillouin zone and n characterizes the unit-cell size. Benchmark calculations on a diverse set of solids─including insulators, semiconductors, metals, and surfaces─demonstrate the fast and robust convergence of k-CIAH-based PMWF optimization, which yields an overall computational efficiency approximately 2-3-fold higher than first-order k-space methods and orders of magnitude higher than Γ-point CIAH for localizing 1000-5000 orbitals. The quality of the resulting PMWFs is further validated by accurate electronic band structures obtained via PMWF-based Wannier interpolation.

Proton-coupled electron transfers (PCETs) are elementary steps in electrocatalysis. However, accurate calculations of PCET rates remain challenging, especially considering nuclear quantum effects (NQEs) under a constant potential condition. Statistical sampling of reaction paths is an ideal approach for rate calculations; however, it is always limited by the rare-event issue. Here, we develop an electrochemistry-driven quantum dynamics approach enabling realistic enhanced paths sampling under constant potentials without a priori defined reaction coordinates. We apply the method in modeling the Volmer step of the hydrogen evolution reaction and demonstrate that the NQEs exhibit more than 1 order of magnitude impact on the computed rate constant, indicating an essential role of NQEs in electrochemistry.

Precise modeling of the energetic landscape is a prerequisite for predicting the charge transport properties of organic light-emitting diodes (OLEDs). However, a significant gap remains between highly accurate but computationally prohibitive self-consistent field (SCF) calculations and efficient but often oversimplified models. In this work, we propose an accurate and effective electrostatic framework with high computational efficiency that encompasses these complex polarization effects through an anisotropically screened dielectric function augmented by a position-dependent background potential. Optimized for the archetypal host material 4,4'-Bis(N-carbazolyl)-1,1'-biphenyl (CBP), our model accurately reproduces the microscopic details, including the polarization-induced stabilization and the surface-reduced energetic disorder, while maintaining high transferability across independent morphological realizations and film thicknesses down to D ≈ 4 nm. Kinetic Monte Carlo (KMC) simulations further confirm that the model faithfully replicates the reference mean squared displacement (MSD) and current-voltage (JV) characteristics, whereas simple image charge models significantly underestimate the current density by failing to describe the downhill gradient at the interface. This framework offers a practical pathway for generating realistic energy distributions for large-scale device simulations, effectively bridging the tradeoff between physical accuracy and computational efficiency.

The valence intermediate effective Hamiltonian (VIEH) approach is a powerful tool for decomposing an overall magnetic coupling strength into ferromagnetic and antiferromagnetic contributions. Herein, we suggest an adaptation of the VIEH framework for complexes with multiple magnetic coupling pathways. This approach is showcased with a homologous series of [M2(μ-O)(NH3)n]2+ complexes, including Cu(II), Ni(II), and Fe(II). Typically, antiferromagnetic coupling progressively weakens from copper to nickel and iron. We demonstrate that this weakening arises primarily from an increase in the ferromagnetic coupling contribution relative to the antiferromagnetic one.

Nonadiabatic effects play a key role in the photophysics and photochemistry of molecular systems, yet their efficient inclusion in quantum molecular dynamics simulations remains challenging, due to the need to construct accurate representations of the molecular Hamiltonian within the manifold of relevant electronic states. Here, we present the PyVCHAM library, which builds on the established multimode vibronic-coupling framework, enhanced by modern machine learning techniques for efficient parameter optimization. The code interfaces with electronic structure packages to generate potential energy surfaces, enabling the parametrization of diabatic Hamiltonians for quantum dynamics calculations. Leveraging the optimization of specialized loss functions, the use of automatic differentiation to compute their analytical gradient in parameter space, and the availability of a wide range of optimization algorithms, the core engine substantially improves in terms of accuracy, flexibility, and efficiency compared to existing implementations. The PyVCHAM library introduces a standardized format for storing vibronic-coupling Hamiltonians based on the JSON data format. It also introduces the ability to combine an arbitrary number of existing vibronic Hamiltonians into interacting supersystems or aggregates, where the constituents couple through dipole-dipole interactions. Several illustrative examples highlight the program's ability to treat complex, high-dimensional molecular systems that were previously difficult to access.

We generalize the Aufbau-suppressed coupled cluster formalism into the realm of doubly excited states by deriving, implementing, and testing a wave function initialization strategy that allows the zeroth-order wave function to match the largest configurations of a doubly excited reference wave function while maintaining the method's overall asymptotic cost parity with ground-state singles and doubles theory. Starting from state-averaged complete active space self-consistent field references, this approach produces highly accurate excitation energies for states dominated by a single doubly excited determinant, as well as states in glyoxal and similar molecules where two different doubly excited determinants have large weights. Typical excitation energy errors in both types of states are on the order of 0.15 eV, with the largest observed error being 0.3 eV. These errors stand in stark contrast to equation of motion methods, where typical errors are 4 to 6 eV at the singles and doubles level and 0.4 to 0.8 eV at the full triples level. It remains an open question how best to generalize the Aufbau suppression approach into an even wider variety of multiconfigurational double excitations, but these early results offer strong motivation for further investigation.

Surrogate models, such as Boltzmann generators (BGs) and emulators (BEs), based on deep generative models are becoming an important tool in molecular simulation. Often, we may want to use additional external information such as sparse experimental data to refine these models. However, there is no unique way to achieve this goal. Here, we propose a method inspired by thermodynamic work from statistical mechanics to regularize the guidance of pretrained probability flow generative models (e.g., continuous normalizing flows or diffusion models) to match additional sparse information. The regularization ensures that the excess work of the guidance procedure is minimized. We developed two guiding strategies based on this method: Path Guidance, which facilitates sampling of rare transition states by concentrating probability mass on user-defined subsets, and Observable Guidance, which aligns generated distributions with experimental observables while preserving entropy. We demonstrate the framework's versatility on two coarse-grained Boltzmann emulators, showcasing its ability to sample transition configurations and to correct systematic biases using experimental data on a variety of model protein systems. Finally, we provide bounds on the distributional differences between the guided and unguided distributions. The method bridges thermodynamic principles with modern generative architectures, offering a principled, efficient, and physics-inspired alternative to standard fine-tuning in data-scarce domains. Our results highlight improved sample efficiency and bias reduction, underscoring their applicability to molecular simulations and beyond.

Ice surfaces play a central role in climate processes, astrochemistry, and materials science, yet their microscopic structure remains elusive. In particular, the degree of proton ordering at ice Ih surfaces critically influences surface reactivity, stability, and phase transitions. In this work, we employ advanced computational techniques─density functional theory to optimize equilibrium geometries, and many-body perturbation theory (GW and Bethe-Salpeter equation) to describe electronic and optical properties─to investigate ordered and partially disordered thin films of hexagonal ice (Ih). First, we analyzed six surface models featuring distinct arrangements of dangling OH bonds, quantified via an order parameter, and computed their Reflectance Anisotropy spectra, which exhibit a pronounced dependence on proton ordering. Among these, two representative models, the Ih-striped and Ih-low-ordered surfaces, emerge as the most stable. For these cases, we demonstrate that proton ordering governs the anisotropy of the optical response: the striped surface supports strongly directional excitonic states, in contrast to the nearly isotropic excitons observed in the low-ordered surface. Our results establish optical anisotropy as a robust fingerprint of proton order, providing a theoretical benchmark for polarization-resolved spectroscopic studies of ice. Furthermore, we show that excitonic effects serve as a sensitive probe of surface proton configurations, paving the way for experimental discrimination between competing models of ice surfaces under cryogenic conditions.