Abstract The mechanisms underlying critical energy-gap fluctuations and nonadiabatic dynamics in disordered condensed phases remain elusive, and their accurate simulation poses a challenge due to the inherent complexity of these systems. In this study, we make the first attempt to leverage highly efficient coarse-grained (CG) models to investigate these fundamental mechanisms. Employing a rigorous bottom-up parameterization approach, we construct CG models for two prototypical systems: a carotenoid–porphyrin–fullerene triad in tetrahydrofuran and a perylene diimide dimer in acetonitrile. Using nonadiabatic semiclassical mapping dynamics, we evaluate a hierarchy of representations spanning fully atomistic explicit-molecule simulations, intermediate CG resolutions, and highly reduced multistate harmonic (MSH) models to assess their impact on energy-gap fluctuations and nonadiabatic properties. A distinct advantage of these explicit-molecule simulations is their ability to yield time-dependent radial distribution functions, revealing the detailed molecular picture of localized solvent rearrangement upon photoinduced charge transfer. Our results confirm that with appropriate coarse-graining of the solvent and solute, the resulting CG models faithfully preserve the structural responses, pairwise reorganization energies, and electronic population dynamics of the atomistic benchmark. Conversely, aggressive 2-heavy-atom solute mapping severely distorts local nonbonded environments. Importantly, MSH models parameterized from structurally valid CG simulations successfully reproduce explicit-molecule dynamics while discarding tens of thousands of degrees of freedom. This work provides a new, highly efficient theoretical foundation for scaling nonadiabatic dynamics simulations to larger and more complex condensed-phase systems.
Abstract Machine learning electron densities enabled fast quantum chemistry but suffered from a fidelity gap in downstream tasks due to integration errors and unbalance on uniform grids. We introduce MARGR, an adaptive framework that refines the grid only in the critical regions. We analyze its efficiency through theory and benchmark it on molecular property calculations, demonstrating that adaptive grid refinement achieves an accuracy comparable to high-resolution uniform grids at substantially reduced computational cost.
Abstract In this work, we have extended the existing “sol-3c” composite methods to meta-GGA XC functionals using the r2SCAN one. The corresponding r2SCAN-based composite methods in their pure and hybrid HF/DFT variants have been developed in combination with double-ζ (sol-def2-mSVP) and triple-ζ (pob-TZVP-rev2) quality basis sets specifically adapted to solids. The r2SCAN “sol-3c/pob-3c” methods are validated and tested on benchmark sets containing both intermolecular adducts (S66 × 7), molecular crystals (X23, DMC8) and simple inorganic solids (SS20). They are then applied to selected applications of relevance for solid-state calculations such as the study of the energetic and structural properties of layered materials, the mechanical, electronic and vibrational properties of α-quartz, relative stability of all-silica zeolite polymorphs, α-quartz vs stishovite phase transition and the adsorption of small molecules on the surface of inorganic solids and in the pore of metal–organic frameworks. The new methods show broad applicability with accurate results and an adjustable and affordable computational cost.
Abstract We present an automated path-integral generator for resonance Raman (RR) calculations supporting Herzberg–Teller effects and Duschinsky rotation, and able to treat any type of vibrational transition, including overtones, combination bands, as well as hot bands. The method constructs the required tensor invariants directly from low-rank intermediates, avoiding both transition-specific manual derivations and the explicit construction of high-order correlation functions. The generated expressions are canonicalized symbolically and evaluated through optimized contraction paths. Compared to formulations based on explicit high-rank correlation tensors, this substantially lowers the asymptotic cost for several higher-order transition classes while providing a unified workflow for automatic expression generation. We validate the approach for multiple transition families, examine the influence of finite-temperature effects through hot-band contributions, and show how it can be straightforwardly extended to resonance Raman optical activity (RROA).
Abstract Position-dependent diffusivities are central parameters in reduced stochastic descriptions of molecular transport in heterogeneous environments, but their reliable estimation from molecular dynamics simulations remains challenging. We present a residence-time approach (RTA) that extracts local diffusivities from first-exit statistics measured in biased simulations after compensation of the mean free-energy gradient. We apply the method to oxygen diffusion across a hexadecane/water slab, water permeation across a POPC lipid bilayer, and transport of water and volatile organic compounds through a model skin-barrier membrane. In the slab system, RTA diffusivities agree with independently determined bulk reference values. In the membrane systems, propagator predictions based on RTA-derived diffusivities reproduce unbiased molecular dynamics propagators over substantial lag-time ranges, while also revealing that, in some cases, no single lag-time-independent diffusivity profile captures the dynamics across all time scales. These results support residence-time statistics as a practical route for determining effective position-dependent diffusivities from biased molecular simulations.
Abstract We report the state-averaged spin–orbit driven similarity renormalization group second-order perturbation theory (SA-SO-DSRG-PT2) for efficient computations of spin–orbit coupling (SOC) effects in multireference systems. This approach adopts a set of multiconfigurational reference states obtained with a spin-free relativistic Hamiltonian and introduces SOC during the perturbative treatment of dynamical correlation. Specifically, the zeroth-order Hamiltonian is defined from the spin-free Hamiltonian, while the remaining spin-free contributions and the spin-dependent Hamiltonian enter as first-order perturbations in the DSRG transformation. This formulation generates an SOC-containing SA-DSRG-PT2 effective Hamiltonian while largely preserving the computational cost of the conventional spin-free theory. Two spin-dependent Hamiltonians are examined: the first-order Douglas–Kroll–Hess spin–orbit Hamiltonian with mean-field approximated two-electron contributions and the spin-dependent part of the one-electron X2C Hamiltonian with the screened-nuclear spin–orbit approximation. Benchmarks on main-group atoms and diatomic molecules, transition-metal elements, trivalent lanthanide cations, and actinide dioxide cations show consistent accuracy across the periodic table for both spin-dependent Hamiltonians. Applications to Co(II) single-ion magnets, including systems described with up to 1790 basis functions, further highlight the promise of SA-SO-DSRG-PT2 for efficient SOC calculations in large open-shell molecules.
Abstract X-ray absorption spectroscopy (XAS) is widely used as an element-specific probe of local electronic and geometric structures, yet its accurate simulation for correlated systems remains challenging because it requires the simultaneous treatment of large active spaces and dynamic electron correlation. In this work, we present a matrix-product-state-based multireference configuration interaction (MPS-MRCI) approach for XAS simulations in such systems. The method is assessed across closed-shell, strongly correlated, and open-shell systems. For pyrazine, systematic calculations reveal the influence of active space size and dynamic correlation on the C and N K-edge spectra. Applications to pentacene, ozone, and the allyl radical further evaluate the performance of MPS-MRCI for a large π-conjugated system, a strongly correlated biradicaloid, and an open-shell system, respectively. The calculated spectra generally reproduce the major experimental features, with natural transition orbital analysis providing insight into the corresponding core excitations through visualization of hole-particle pairs. This work establishes MPS-MRCI as a viable approach for core-level spectroscopy of correlated systems, enabling the treatment of large active spaces and the recovery of dynamic correlation.
Abstract The synthetic miniprotein chignolin features a small size, a well-defined native structure, and a rather paradigmatic two-state folding transition. Nonetheless, this deceivingly simple molecule showcases a nontrivial folding pathway, whose detailed characterization has been the subject of several studies and still presents some open questions. In this work, we investigate chignolin by employing, for the first time in a combined and integrated pipeline, transition-path theory (TPT) and the mapping entropy optimization workflow (MEOW). The former describes the folding process in terms of its natural reaction coordinate, that is, the committor; the latter pinpoints, in an unsupervised and system-agnostic manner, those residues of a biomolecule that are most informative about the structural, mechanical, and energetic organization of the configurational ensemble. Through this joint framework, we characterize in great detail the entire transition of the peptide from the unfolded to the native state. The approach allows us to identify which residues entail the largest degree of information about a conformational ensemble at a given level of progress of the folding transition; comparison with data from the literature and validation against independent observables, including the outcomes of an in silico mutation analysis, show that these residues also bear functional significance. This work analyzes the folding process of chignolin from a new perspective and showcases the integration of the MEOW protocol with TPT into a pipeline that can complement existing approaches for the investigation of proteins.
We propose an improved model, termed the gradient-corrected PCM (GCPCM), for improving the energy accuracy of the polarizable continuum model (PCM). Our previous study revealed deficiencies of PCM in describing the reaction field, i.e., the electrostatic potential generated by the solvent. These deficiencies can be partially alleviated by introducing an empirical correction to the solvent charges. As a result, solute-solvent interactions are improved at the self-consistent field level, leading to enhanced energy accuracy. The performance of GCPCM was evaluated through single-point calculations and geometry optimizations of phenol and phenolate, calculations of the free energy profile for proton transfer in glycine, and analysis of solvent responses of the HOMO and LUMO orbital energies of Brooker's merocyanine. The results demonstrate that the characteristic destabilization of charged solutes observed in conventional PCM is effectively resolved. Furthermore, despite having a computational cost comparable to that of PCM, GCPCM shows the potential to achieve an energy accuracy similar to that of 3D-RISM-SCF. The development of GCPCM enables more convenient and accurate treatment of solvation effects, which is expected to allow researchers to focus on other important challenges, such as the accurate description of electronic states.
Electrocatalytic machine-learning potentials must simultaneously describe long-range electrostatics, nonlocal charge redistribution, and electrode-potential-dependent interfacial response, which makes their physical construction and validation particularly demanding. Here, we introduce DPχ, a charge-based machine-learning potential designed for electrified metal-water interfaces. DPχ represents long-range electrostatics through Bader-basin centroids and decomposes interfacial charge into a neural-predicted chemical component and a conductor component determined self-consistently by a Siepmann-Sprik-type polarizable-electrode model under global electroneutrality. Rather than claiming broad transferability across electrocatalytic materials, we test these physical assumptions on the benchmark Pt(111)-water interface. Systematic benchmarking shows that DPχ reproduces DFT-level forces, interfacial potential drops, hydrogen-coverage-dependent electrode potentials, Volmer barriers, and interfacial vibrational signatures, while remaining robust upon system-size enlargement. These results establish DPχ as a physically consistent and reaction-ready framework for large-scale simulations of the Pt(111)-water electrochemical interface beyond AIMD spatiotemporal scales.
Redox enzymes play an essential role in nature and in biotechnological applications such as (photo)biocatalysis and biosensing. Understanding how a protein's sequence and structure tune its redox potential is very valuable for engineering proteins with tailored (photo)redox properties. Since protein redox potential measurements are laborious, reliable redox potential computations present an attractive alternative. However, redox calculations come with their own set of challenges, such as ensuring adequate sampling and accurate force fields. Dealing with charge-changing states introduces additional challenges to theory and requires special considerations in the model setup. Here, we report a protocol for computing the change in the proton-coupled one-electron redox potential associated with a D63N charge-changing mutation in a prototypal flavoprotein, Desulfovibrio vulgaris flavodoxin. An automated average protein electrostatic configuration protocol, APEC-F 2.0, was used to construct hybrid quantum mechanics/molecular mechanics (QM/MM) models. These models were used for subsequent alchemical free energy simulations in which a charged surface aspartate residue was gradually converted to an isosteric but neutral asparagine (D63N) over 40 λ windows. A thermodynamic cycle was employed to calculate the redox potential of the mutant relative to the wild-type reference. This calculation was repeated for the same D63N mutation using models prepared under slightly different conditions, focusing primarily on factors that affect the electrostatic environment in the system. Factors tested include (1) the effect of including extra salt ions in the model solution, (2) different protocols to balance the disappearing negative charge associated with the alchemical D63N mutation, and (3) the effect of accounting for the kinetic energy terms in the free energy protocol due to alchemical morphing of the atomic masses. The results indicate that such apparently minor details may have a considerable effect on the random and systematic errors obtained from the free energy simulations. The best models were shown to reproduce the experimental shift in the redox potential due to the D63N mutation with an accuracy of 0.3 kcal/mol. The associated error for this shift is 1.3 kcal/mol, calculated as the total standard deviation of triplicate simulations for each of the oxidized and reduced neutral semiquinone states of flavin. The addition of salt ions in the simulations and the proper treatment of charge-conserving coalchemical counterions in particular are found to be paramount, since less accurate models led to errors more than double the magnitude compared to the best protocol. The sources of those errors are discussed and often found to be associated with medium-range electrostatics due to missing or inaccurate second solvation/ionic shells around the mutation site.
Noncovalent interactions play a central role in chemistry, biology, and materials science, governing processes ranging from molecular recognition and protein folding to crystal packing and supramolecular assembly. The rational design of these interactions relies on principles such as directionality and complementarity, which control the organization and stability of molecular complexes. Although energetic decomposition analyses can provide insight into these features, their computational cost limits their application to large-scale systems. Here, we introduce a simple density-based descriptor,qNCIVvdW, that quantifies interaction localization and captures the balance between electrostatic and dispersion contributions. The descriptor can be evaluated globally for complete complexes or locally for individual interaction regions, requiring only structural information and avoiding computationally demanding energy partitioning schemes. This provides a scalable framework for characterizing and designing noncovalent interactions across molecular systems ranging from small complexes to large supramolecular architectures.
The high computational cost of deriving REPEAT charges via periodic density functional theory (DFT) limits the large-scale screening of metal-organic frameworks (MOFs). To address this, we developed DeGAT, a dual-expert graph attention network for the rapid prediction of partial atomic charges. By incorporating an uncertainty-driven active learning strategy on the ARC-MOF database, the model achieves a test-set R2 of 0.985, with a mean absolute error (MAE) of 0.0314 e and a species-averaged MAE (SMAE) of 0.0527 e. Subsequent grand canonical Monte Carlo and Widom insertion simulations demonstrate that CO2, N2, and water adsorption properties calculated using DeGAT charges closely match those derived from standard REPEAT charges. These results demonstrate that DeGAT enables efficient and accurate prediction of partial atomic charges in MOFs while maintaining charge neutrality and physical consistency, providing a scalable parametrization scheme for high-throughput screening and molecular simulations of porous materials.
Time-dependent density functional theory (TDDFT) has been widely used to model electronic excitations but remains unreliable for charge-transfer and doubly excited states, core excitations, and the complex topology near conical intersections (CI). We introduce a simple exciton model to arbitrary open-shell singlet excited (OSE) states by incorporating nonlocal singlet correlation energy, derived from the corresponding triplet excited state, thereby bypassing the direct optimization of singlet states and approximate spin-projection schemes. This model achieves near-quantitative accuracy for low-lying valence and Rydberg excitations and successfully reproduces the absorption spectrum of Chlorophyll a and zinc phthalocyanine. For core excitations, we found an excellent agreement across K and L-edge excitations of first and second group elements and also successfully capture the distinct experimental X-ray absorption fingerprints of imidazole-imidazolium system in water without any empirical shifting. More importantly, this model also captures near-degeneracy behavior in the vicinity of the CI for photoisomerization of cis-trans azobenzene, where standard TDDFT fails. These findings, collectively, position this simple model as a reliable and alternative ΔSCF-like approach for modeling nontrivial excited-state phenomena.
Predicting the performance of batteries using analytical and computational models plays an important role in the design of battery packs and management systems. Currently, these models rely on extensive experimental parametrization; but fundamentally, these parameters arise from atomistic interactions between components of the electrolyte and the resulting correlated motion of the ions. In this work, we demonstrate a comprehensive approach to atomistic simulations of discharging batteries, evaluating electrochemical potential gradients in the electrolyte and using Onsager mobility coefficients to relate the resulting forces to the flux of lithium ions between the electrodes. This work unifies four different sets of simulations: (1) mobility, which observes the molecular flux of species in response to constant forces; (2) thermodynamic susceptibility, which observes the response of species concentration to external potentials; (3) bulk modulus, which observes the density response to pressure; and (4) battery discharge, which generates a steady flux of cations and observes the resulting concentration and potential gradients. The resulting model self-consistently describes ion-transport battery electrolytes.
Predicting spectra and photochemical pathways relies on the computation of molecular excited states. However, the practicality of near-term variational quantum algorithms is threatened by the prohibitive growth of measurement overheads. Integrating Quantum Subspace Expansion with the Contextual Subspace (CS) method, we propose a CS-QSE framework to resolve the scaling bottlenecks in excited-state energy estimation. By confining excitation operators to a compact CS, the framework removes redundant degrees of freedom and reduces Pauli string counts, thereby easing the measurement burden. The primary strength of CS-QSE lies in its substantial reduction of the operator-pool size. Benchmarking on LiH, HF, H2O, and HCl shows that CS-QSE achieves errors within the target tolerance of 1.6 × 10-3 Ha relative to full configuration interaction benchmarks in the same basis set while mitigating the prohibitive scaling inherent in the full-space method. Numerical simulations reveal that the operator pool size is consistently pruned by over 90% across all systems, with the reduction reaching as high as 99.7% for molecules such as HCl. This CS-QSE framework establishes a resource-efficient route for molecular simulations tailored to near-term noisy intermediate-scale quantum devices.
The performance of batteries is heavily influenced by the properties of their solvents, which play a crucial role in ion solvation, conductivity, and electrochemical stability. However, traditional force field models often fall short in accurately predicting these properties. This study investigates the potential of adaptive force matching (AFM) to predict a range of solvent-related properties using only electronic structure theory. To demonstrate this approach, we developed AFM models based on B3LYP-D3(BJ) for three common molecules present in battery solvents: dimethylamine (DMA), dimethyl carbonate (DMC), and tetrahydrofuran (THF). Our results reveal that AFM models significantly outperform traditional force fields, including OPLS-AA (CM1A), GAFF2, and GROMOS 54A7, in predicting key properties such as density, heat of vaporization, viscosity, diffusion constant, boiling temperature, and free energy of vaporization. Notably, the average percentage error of AFM models is approximately 13%, substantially lower than that of traditional force fields (21% for GAFF2, the best-performing empirical model). These findings underscore the promise of AFM in predicting physical properties of battery solvents, with far-reaching implications for chemistry, materials science, and related fields where accurate predictions are essential for understanding complex phenomena and designing innovative materials and systems.
In this pilot study, structural properties and 1H NMR chemical shifts of an Ala-Glu-Pro-Phe peptide dissolved in an aqueous solution and an aqueous mixture of the 1-ethyl-3-methylimidazolium trifluoroacetate ionic liquid (IL) have been scrutinized using an integrated computational protocol based on classical molecular dynamics (MD) simulations and combined quantum mechanics/molecular mechanics (QM/MM) approaches for NMR shielding constants. Two long-lived trans and cis Pro isomers of the tetrapeptide have been considered. MD simulations as long as 400 ns were found to be too short to ensure complete sampling of the conformational phase space of the peptide in aqueous IL solution, and thus, the analysis was performed for four distinct conformations of either isomer of the peptide separately. Solvent molecules within the first solvation shell of the solute were treated quantum mechanically in the QM/MM calculations of 1H NMR shieldings. The constituent ions of the IL were seen to condense around the tetrapeptide, abundantly displacing water molecules from its first solvation shell, and the preference for the imidazolium cations to condense around the peptide in solution was identified. Prominent hydrogen bonding together with π-π stacking interactions between the benzene ring in the Phe residue and the imidazolium ring of the cations assist in maintaining the ionic cage that surrounds the peptide. Accordingly, the computed 1H NMR signals of the peptide in aqueous IL solution are shifted w.r.t. their values in aqueous solution. The computational results allow identifying several 1H NMR signals which could be potential spectral NMR markers of the conformational or the isomeric state of the tetrapeptide.
We propose a general strategy to discretize the Dyson series without applying direct numerical quadrature to high-dimensional integrals and extend this framework to open quantum systems. The resulting discretization can also be interpreted as a Strang splitting combined with a Taylor expansion. Based on this formulation, we develop a deterministic iterative method for simulating system-bath dynamics. We propose two numerical schemes, which are first-order and second-order in the time step Δt, respectively. In the second-order scheme, we can safely omit most terms arising from the Strang splitting and Taylor expansion while maintaining second-order accuracy, leading to a substantial reduction in computational complexity. For the second-order method, we achieve a time complexity of O(M322KmaxKmax2) and a space complexity of O(M222KmaxKmax), where M denotes the number of system levels and Kmax the number of time steps within the memory length. Compared with existing methods, our approach requires substantially less memory and computational effort for multilevel systems (M ≥ 3). Numerical experiments are carried out to illustrate the validity and efficiency of our method.
We compared the performance of two mass repartitioning models, HMR3 and HMR2, with tripled and doubled hydrogen masses, against the model with standard masses (SM). Due to heavier hydrogens, HMR3 and HMR2 afford longer integration steps of 4 and 3.5 fs, respectively, whereas SM used a 1 fs time step as a reference. All-atom replica exchange molecular dynamics simulations of the antimicrobial peptide PGLa binding to an anionic DMPC/DMPG bilayer were used as a case study. Across all conditions, HMR3 and HMR2 simulations each collected 36 μs of sampling, which is 50% more than our previous SM simulations. Therefore, the motivation for our work was to evaluate the HMR performance in complex biomolecular systems coupled with an advanced sampling algorithm. Our investigations led to three main conclusions. First, by affording longer simulations, HMR3 and HMR2 models provide better equilibration of PGLa peptide binding to the lipid bilayer. Specifically, they eliminated the metastable surface-bound (SB) state observed in SM. In the mass repartitioned models, PGLa exclusively sampled the inserted state (I), which also appeared as a dominant state in SM but alongside the SB state. Importantly, this better equilibration of peptide binding improved the consistency with experimental data. Second, HMR3 and HMR2 closely reproduce a broad range of structural properties previously reported for SM after correcting them for equilibration. Additionally, HMR3 and HMR2 demonstrate excellent agreement between themselves in sampling the conformational ensemble. These findings indicate that the mass repartitioning models preserve structural properties in this complex biomolecular system. Third, the actual acceleration of sampling by HMR3 is about 3.5-fold, which is lower than the theoretical 4-fold gain. This slight underperformance may be due to the slow dynamics of heavy hydrogens. HMR2 shows qualitatively similar results. An increase in the integration step in SM to 2 fs or adjustments in the computations of long-range interactions may reduce the computational gains offered by mass repartitioning, but do not negate their advantages in efficiency. HMR2, and particularly HMR3, are excellent options for accelerating conformational sampling in complex biomolecular systems.
