Liquid crystal (LC) phases formed by anisotropic particles have long attracted interest due to their unique combination of fluidity and directional order, as well as their prevalence in natural systems. Among them, chiral entities such as helices exhibit exotic LC phases, like the cholesteric and screw nematic, in addition to isotropic, smectic, and crystal phases that are observed in systems of achiral particles. Chiral particles are ubiquitous in nature across length scales. Helices, being the simplest examples of chiral particles, are important components in a variety of biological systems in the form of DNA, protein fragments, and others. In this work, we employ molecular dynamics to investigate the thermodynamics and structural characteristics of LC phases formed by a system of chiral particles compared to those formed by a system of achiral particles of the same aspect ratio. We consider a system of soft chiral rods modeled as made of fused beads and a system of soft repulsive spherocylinders (SRS) for this comparison. We evaluate the role of chirality in phase behavior with special emphasis on cholesteric and screw nematic phases. Extending the two-phase thermodynamic (2PT) model, we compute entropy, free energy, and fluidicity parameters to understand the thermodynamic stability of the similar phases exhibited by both systems. Our results reveal significant differences in phase stability and structure between the two systems, with crystals of achiral particles having a higher packing fraction than crystals of chiral particles. The study also demonstrates that translational and rotational fluidicity parameters can serve as effective phase identifiers, offering a simplified approach to characterizing complex LC phases. These findings deepen our understanding of chirality-induced phase behavior and also pave the way for applying fluidicity-based analysis to other anisotropic systems.
Proteins exhibit a diverse range of structures and dynamics that are critical to their biological function. These dynamic processes span a broad spectrum of time scales and are influenced by environmental factors, including temperature, solvent composition, and the presence of binding partners or membranes. Efficient exploration of protein conformational space is essential for understanding their functional mechanisms, but this remains challenging because of the high dimensionality of the energy landscape. Our group has previously developed the molecular dynamics with excited normal modes (MDeNM) method, which is based on the kinetic excitation of normal modes (NMs) during molecular dynamics simulations. Here, we developed an adaptive extension of the method (aMDeNM), where the motions described by preselected directions of low-frequency NMs are dynamically adjusted throughout the simulation. By coupling low-frequency NM excitation with adaptive directional adjustments, aMDeNM facilitates extensive exploration of the energy landscape, overcoming the constraints of fixed, rectilinear displacements and alleviating structural stresses and environmental resistance. The method was tested on three structurally diverse test systems: T4 lysozyme, human calmodulin, andStaphylococcus aureus monofunctional transglycosylase. Our results demonstrate improved conformational sampling compared with standard MD and other enhanced sampling methods. Additionally, spectral analysis of structural oscillations along the pathways using fast Fourier transform revealed the role of low-frequency vibrations in critical conformational changes and highlighted the influence of the surrounding environment on protein dynamics. This work provides a robust framework for studying large-scale protein motions and their functional implications within complex biological environments. Importantly, aMDeNM requires only an initial structure without the need to specify predefined target states, distinguishing it from many biased sampling techniques that rely on predefined target conformations. The aMDeNM code and usage instructions are available at https://github.com/pedro-tulio/amdenm.
Accurate prediction of drug-target binding and unbinding kinetics and thermodynamics is essential for guiding drug discovery and lead optimization. However, traditional atomistic simulations are often computationally expensive to capture rare events that govern ligand (un)binding. Several enhanced sampling methods exist to overcome these limitations, but they require extensive manual intervention and introduce variability and artifacts in free energy and kinetic estimates that limit high-throughput scalability. The present work introduces seekrflow, an automated multiscale milestoning simulation pipeline that streamlines the entire workflow from a single receptor-ligand input structure to kinetic and thermodynamic predictions in a series of steps. This integrated approach minimizes manual intervention, reduces computational overhead, and enhances the reproducibility and accuracy of kinetic and thermodynamic predictions. The accuracy and efficiency of the pipeline are demonstrated across multiple receptor-ligand complexes, including inhibitors of heat shock protein 90, threonine-tyrosine kinase, and trypsin, with predicted kinetic and thermodynamic parameters that closely match experimental estimates. seekrflow establishes a new benchmark for automated and high-throughput physics-based predictions of the kinetics and thermodynamics.
Nuclear magnetic resonance (NMR) spectroscopy is a powerful tool for characterizing the structure and electronic properties of paramagnetic coordination compounds. However, accurate computation of paramagnetic NMR (pNMR) chemical shifts for large systems remains a major challenge due to the prohibitive cost of full first-principles calculations. Herein, we extend the eXtended ONIOM (XO) method to paramagnetic systems and develop an XO-pNMR approach for the efficient and precise calculation of 13C pNMR chemical shifts in large paramagnetic molecules. We validate the XO-pNMR method by investigating a series of cobalt(II) porphyrins, including cobalt(II) tetraphenylporphyrin (CoTPP) and its substituted derivatives, all of which possess a single unpaired electron. Benchmark calculations against experimental 13C chemical shifts at 308 K reveal that the CAM-B3LYP functional yields the closest agreement with experimental data in full-system calculations, whereas the PBE0-1/3 functional exhibits optimal performance in XO-pNMR calculations─achieving a mean absolute deviation of only 0.3 ppm relative to full calculations and demonstrating robust wave function stability that outperforms alternative functionals. Application of XO-pNMR to substituted Co(II) porphyrins further demonstrates that the method reliably captures substituent-induced variations in 13C chemical shifts, matching the predictive accuracy of full calculations. Collectively, our results establish XO-pNMR with the PBE0-1/3 functional as a cost-effective and reliable approach for calculating pNMR chemical shifts in large paramagnetic molecules featuring a single paramagnetic center with one unpaired electron.
We constructed a benchmark dataset (H2X100) of 100 CCSD(T)/Complete Basis Set (CBS) binding energies for the (H2X)n, X = S, Se, Te, n = 2-4 clusters, and additional CCSD(T)-quality structures and harmonic vibrational frequencies for 12 unique homodimer minima, containing a mix of the hydrogen and chalcogen bonding motifs. Using the new H2X100 dataset, we assessed the performance of 17 density functionals across the rungs of Jacob's ladder, with and without dispersion corrections, spanning the GGA, mGGA, GH-GGA, GH-mGGA, RSH-GGA, and RSH-mGGA groupings. The r2SCAN-D4, B97M-rV, and ωB97M-V functionals provide the best agreement with the CCSD(T) benchmarks and can be used in studies of larger ensembles of these systems, including simulations of their condensed phases.
In electronic structure calculations, the transcorrelated method consists of transforming the Hamiltonian so as to remove the Coulomb cusp in its eigenfunctions. As a result, the wave function can be described more accurately without increasing the size of the basis set. However, the transcorrelated Hamiltonian is non-Hermitian and nonnormal, which makes many common quantum algorithms inapplicable. Recently, a quantum eigenvalue estimation (QEVE) algorithm was proposed for non-Hermitian Hamiltonians with real spectra [FOCS 65, 1051 (2024)]. Although the asymptotic scaling of this algorithm with the desired accuracy is shown to be optimal, the constant factor in its complexity scaling has not yet been analyzed. Here, we investigate the cost of QEVE applied to transcorrelated electronic Hamiltonians of Li, Be, B, C, and N atoms and compare it to the cost of applying standard qubitization to nontranscorrelated Hamiltonians. We find that with the xTC approximation, the T gate count of QEVE in the minimal STO-6G basis is between those of standard qubitization in the cc-pVTZ and cc-pVQZ bases. The accuracy of the transcorrelated energy differs between systems: for Li and Be, it is more accurate than the cc-pVQZ energy, while for larger atoms, the error is between those of the cc-pVDZ and cc-pVTZ energies.
High-throughput computation is essential for the rational design of heterogeneous catalysts, and the quality of the initial adsorption configurations directly determines the efficiency and success rate of subsequent geometry optimization. Existing heuristic methods are constrained by rigid-body approximations, making it difficult for them to handle complex multidentate adsorption, and the low-precision force field correction schemes on which they rely lack generality. To address this issue, this work proposes a general algorithm, High-Throughput Initialization of Multidentate Adsorption Configurations (HiMac). The algorithm reformulates configuration generation as a multiobjective optimization problem. It employs forward kinematics to model molecular flexibility and combines it with a loss function based on parent molecule similarity, thereby unifying site selection with pose adjustment. A statistical learning module is further integrated to prioritize the exploration of chemically favorable sites and reduce the computational cost. Experiments show that HiMac generates high-quality initial configurations for geometry relaxation and transition-state searches, is applicable to arbitrary adsorbate-surface systems, and provides a general and efficient method to accelerate the rational design of heterogeneous catalysts.
The human serotonin transporter (hSERT) is a member of the neurotransmitter:sodium symporter (NSS) family that mediates active reuptake of serotonin from the synapse into the presynaptic neuron. During the transport cycle, hSERT alternates between outward-facing (OF) and inward-facing (IF) states to translocate its substrate between the two sides of the membrane. During the OF-to-IF state transition, serotonin (aka, 5-HT) is inwardly symported together with Na+ and Cl- ions. The return to the OF state is facilitated by cytosolic K+ binding, a step that is also proposed to act as a kinetic decision point by frustrating the outward transport of 5-HT in the direction opposite to the physiological direction of the cycle. However, as opposed to the Na+ ions, the mechanism of K+ binding, its binding site and regulation have not been thoroughly studied. Moreover, recent studies have challenged the conventional transport stoichiometry (1 5-HTin:1 Nain+:1 Clin-:1 Kout+) in hSERT, suggesting that Cl- might remain bound to the transporter during the entire cycle. To explore the role of cytosolic K+ binding to IF hSERT, we performed an extensive set of molecular dynamics simulations. Starting from the post-release IF conformation and in the presence of cytosolic K+, we generated 50 independent trajectories, each for 200 ns to study the behavior of ions. In more than half of the simulations, spontaneous K+ binding was observed at the Na2 site, a conserved cation-binding site in NSS transporters that has been implicated in controlling conformational transitions. Markov state model analysis of coupled ion dynamics quantifies K+ binding kinetics and identifies K+ occupancy of the Na2 site, with Na+ retained at Na1, as the thermodynamically dominant post-release state. In addition, Cl- remains bound to hSERT in the majority of sampled simulations, consistent with recent experimental observations and suggesting a limited role of Cl- release during this stage of the transport cycle. Together, these results provide a kinetic and mechanistic framework for understanding cytosolic K+ binding to hSERT and its potential role in facilitating the IF-to-OF transition that resets the transport cycle.
Unitary Coupled Cluster (UCC) theory is a promising variational method for electronic structure calculations, particularly for systems that exhibit strong electronic correlation and for implementation on quantum computers. However, its practical application is limited to small chemical systems with small basis sets due to its steep computational scaling, which results from its nonterminating Baker-Campbell-Hausdorff expansion. Here, we introduce an active space UCCSD(4)/MP2 approach that leverages a fourth-order many-body perturbation theory truncation of UCCSD within a selected active space while treating external excitations at the MP2 level. We explore two variants: a composite method that sums separate internal and external contributions and an interacting method that couples the amplitudes for potentially greater accuracy. We test our approach on a range of systems, including molecules from the GW100 data set in their equilibrium geometries, a moderately correlated metaphosphate hydrolysis reaction, and the strongly correlated torsion of ethylene. Our results suggest that the interacting method with canonical orbitals is robust and stable for both weakly and moderately correlated systems and accurately reproduces the full UCCSD(4) potential energy curves, including only 15-25% of the virtual orbitals in its active space. In comparison, the composite formulation exhibits greater sensitivity to the choice of orbital basis and active space size, leading to less systematic behavior across the benchmark set. For ethylene torsion, a system dominated by strong static correlation, both the composite and interacting formulations employing canonical orbitals closely track the full UCCSD(4) reference while preserving the qualitative behavior of the parent method, including the breakdown in strongly multireference regimes. This active space framework offers a computationally tractable approach for modeling correlated molecules and reactions on classical computers and provides a viable path for scaling UCC calculations for resource-constrained quantum hardware.
We present a unified implementation ready computational framework for rotation-averaged, frame-invariant evaluation of the first- (β) and second-order (γ) molecular hyperpolarizabilities in isotropic media. Starting directly from ab initio Cartesian tensors as delivered by standard electronic-structure codes, these tensors are projected onto Frobenius-orthogonal irreducible components (β: J = 1, 3; γ: L = 0, 2, 4), defining isotropic invariants (B1,B3; N0, N2, N4). Using isotropic tensor algebra, we derive closed-form expressions for the HRS channels ⟨|βZZZ |2⟩ and ⟨|βZXX |2⟩ and, within a fully Cartesian M6/M8 invariant formulation, for the THS channels ⟨|γZZZZ |2⟩ and ⟨|γZXXX|2⟩ as exact rational combinations of N0,N2,N4 consistent with established THS tensor analysis, the mixed ZXXX channel is insensitive to the scalar (L = 0) sector. These formulas are directly useable for HRS and THS, yielding third-harmonic scattering intensities and depolarization ratios directly from invariant norms without invoking symmetry assumptions. To ensure robustness of implementation, we employ three mutually consistent routes: (i) analytic SO(3) contractions with M6/M8 tensors as the reference solution; (ii) deterministic Euler-angle quadrature; and (iii) uniform-quaternion Monte Carlo with variance-based stopping criteria used exclusively as independent verification and quality control. Beyond scalar averages, we compute directional power maps Qβ2 (n), Qγ2 (n) and introduce charge-transfer normalized effective hyperpolarizabilities βeffCT and γeffCT, which scale invariant RMS amplitudes by the S0 → S1 charge-transfer length to quantify the nonlinear response per CT length, and D-π-A chromophores confirm cross-validated β/γ metrics for HRS/THS. Windows-friendly, turnkey code is provided to facilitate reproducible use by computational practitioners.
Cysteine is a chemically versatile amino acid whose protonation state is central to catalysis, redox regulation, and covalent ligand recognition. Despite the growing availability of constant-pH molecular dynamics (CpHMD) methods, cysteine is not included among the standard titratable residues of the GROMACS-based λ-dynamics implementation developed by Aho, Buslaev, and co-workers. Here, we introduce a titratable cysteine residue, CYST, for this framework and integrate it into the phbuilder workflow. The implementation describes the thiol/thiolate equilibrium of reduced cysteines within a fixed redox topology, while disulfide-linked cysteines are retained as nontitratable covalent states defined a priori. The new residue was calibrated on an ALA-CYS-ALA (ACA) tripeptide by refining the correction potential required to remove the intrinsic force-field bias along the λ coordinate. The resulting model reproduces the expected sigmoidal titration behavior of the peptide and yields a pKa of 8.33, matching the target value adopted in the calibration procedure. Transferability to proteins was then assessed using two complementary benchmarks. First, six engineered single-cysteine mutants of acyl-coenzyme A binding protein (ACBP) were used as a set of noncatalytic cysteines spanning different local environments. The calculations reproduce the overall high-pKa regime of these residues, while also showing that quantitative accuracy depends sensitively on the local conformational ensemble sampled around the introduced cysteine. Second, the implementation was applied to the reactive Cys106 of DJ-1, placing its titration in the correct acidic regime. Overall, this work establishes a practical cysteine extension of the GROMACS CpHMD framework for biomolecular simulations.
This work introduces the Hierarchical Neighborhood of Atoms (HNA) partitioning, a recursive classification of chemically equivalent atoms that enables ultrafast computation of exact symmetry-corrected atomic correspondence and RMSD for molecular conformers. By decomposing the assignment problem hierarchically, the algorithm reduces the number of evaluated combinations from the product of branch possibilities to their sum, yielding reductions of up to 16 orders of magnitude for systems such as myoglobin. Topology-unaware approaches based on linear assignment produce chemically invalid atomic correspondences in over 89% of protein conformer pairs, while topology-aware graph isomorphism methods time out on 51-66% of biologically relevant molecules. In contrast, the proposed method achieves 100% topologically correct assignments without timeouts across all 1.4 million pairs tested, with millisecond-scale mean execution times (1.3-3.8 ms for the CCD data set, 2.7-16.5 ms for the BIRD data set) and 11-42× speedups over polynomial-time methods on protein conformer ensembles. Beyond RMSD, the HNA framework can be used for symmetry-consistent comparison of atomic properties, canonical atom labeling, and force field parametrization.
Constant-potential molecular dynamics is essential for realistic simulations of electrochemical interfaces under Operando conditions. Although various constant-potential frameworks exist, most are tightly coupled to specific electronic-structure codes or numerical architectures, limiting portability and extensibility─especially for codes constrained to integer electron numbers. Here, we present a flexible constant-potential framework implemented in the i-PI driver, interfacing with multiple density functional theory (DFT) engines and, in principle, extensible to constant-potential machine-learning potentials. The method regulates and samples the electronic chemical potential by introducing an explicit electronic degree of freedom and a dedicated potentiostat module in i-PI. To bypass the integer-electron constraint without modifying the underlying DFT code, we employ a mixed-Hamiltonian interpolation scheme: two adjacent integer-charge clients are run in parallel, and their energies, forces, and electronic chemical potentials (Fermi level/work function) are linearly interpolated to obtain an effective fractional-charge description. We validate the method on a one-dimensional asymmetric double-well model and an Al(111) surface, demonstrating stable potential control and well-behaved charge fluctuations. Finally, we couple constant-potential ab initio molecular dynamics (AIMD) with enhanced sampling to study CO2 reduction on NiN4-doped graphene, enabling efficient characterization of potential-dependent reactivity and free-energy landscapes. Overall, this framework provides a portable and scalable platform for conducting rigorous constant-potential simulations across diverse electronic-structure clients and, in principle, machine-learning potentials.
Marcus theory provides a fundamental framework for describing hopping-type charge transport, and donor-bridge-acceptor (D-B-A) dimer models are widely used to evaluate reorganization energies and electron transfer rates. In practice, constrained density functional theory (CDFT) has been extensively employed to construct charge-localized states because the direct evaluation of diabatic states using post-Hartree-Fock (post-HF) methods is computationally demanding. However, in systems containing π-conjugated bridges, maintaining charge-localized states within CDFT can become more difficult, which may affect the evaluated reorganization energies. Standard remedies, including partial geometry fixation and ONIOM-based approaches, do not sufficiently resolve this limitation. In this study, we propose a local charge state fixation (LCSF) method for calculating reorganization energies in D-B-A dimers without imposing explicit constraints on the electronic density. Charged donor and acceptor monomers are optimized independently and subsequently translated and rotated to satisfy the Eckart condition with respect to the corresponding fragments of the neutral dimer. These monomers are then grafted to construct charge-localized dimer geometries, enabling reorganization energy evaluation without charge flow. The method was applied to 1,4-diphenylbutane (PS), trans,trans-1,4-diphenyl-1,3-butadiene (PD), and 1,4-diphenylbutadiyne (PT). For the PS, reorganization energies obtained from LCSF and CDFT are comparable, whereas for PD and PT, LCSF yields consistently larger values, indicating reduced electron delocalization effects. The applicability of the method was further demonstrated using Hartree-Fock and MP2 calculations. In addition, LCSF was applied to zinc oxo cluster dimers, allowing electron transfer rates to be evaluated. Although the method introduces structural constraints, it offers a computationally efficient strategy for reorganization energy calculations in dimer systems, particularly within post-HF frameworks.
The embedded cluster reference interaction site model (EC-RISM) employs statistical-mechanical integral equation theory to predict solvent site distributions and their interaction with a solute using quantum-mechanical electronic structure methods. In contrast to apparent surface charge models such as the polarizable continuum (PCM) or the conductor-like screening (COSMO) model, EC-RISM can account for directional solvent-solute interactions due to the description of the solvent based on conventional molecular force field models. Here we present an implementation of EC-RISM combined with correlated wave function methods for ground- and excited-state energies and, for the first time, also for ground- and excited-state energy gradients. This is achieved by self-consistent equilibrating the solvent reaction field with the solute charge density in the correlated and, possibly, electronically excited state. To account for excitonic coupling and nonequilibrium contributions to electronic transition energies, EC-RISM is further combined with COSMO. We present applications to the molecular structures and the absorption and emission energies of 4-(N,N-dimethylamino)benzonitrile (DMABN), the photobase 7-amino-4-methylcoumarin, and the photoacids phenol and 3-cyanophenol in aqueous solution. As expected, for systems without strong directional solvent interactions, such as DMABN, EC-RISM yields results similar to those obtained with COSMO whereas, for molecules or ions that form strong hydrogen bonds to the solvent (particularly the deprotonated photoacids) EC-RISM provides substantial improvements.
Near-degenerate electronic structures remain a major challenge for conventional single-reference density functional theory (DFT). To address this problem, we propose time-dependent ΔSCF (TDΔSCF), a novel linear-response scheme in which a non-Aufbau ΔSCF determinant serves as the reference for a subsequent TDDFT calculation. In contrast to collinear spin-flip (SF)-TDDFT, this formulation preserves the usual Coulomb and exchange-correlation response contributions while describing the target states from an electronically excited reference. We examine the performance of TDΔSCF for several prototypical problems involving near-degeneracy, including the torsional potential of ethylene, singlet-triplet gaps of representative diradicals, geometry optimizations of benzyne isomers, and bond-dissociation curves of hydrogen fluoride and F2. Across these tests, TDΔSCF shows markedly weaker functional dependence than SF-TDDFT and often yields a more balanced description of challenging singlet states. In particular, it provides smooth torsional potentials, improved singlet-triplet gaps, a consistent monocyclic structure for singlet m-benzyne, and a more satisfactory description of bond dissociation without the spurious low-lying states found in SF-TDDFT. At the same time, the method exhibits a systematic tendency to overestimate singlet energies and can lose accuracy when the underlying ΔSCF reference is not well suited to the final state. We also identify a numerical instability that can arise in non-Aufbau calculations and trace its origin to the exchange-correlation potential near uncompensated nodal regions. These results highlight both the promise and the practical limitations of TDΔSCF as a low-cost method for singlet states with near-degenerate electronic structures.
Predicting the impact of mutations on protein-ligand binding affinity is crucial in drug discovery, particularly in addressing drug resistance and repurposing existing drugs. Conventional structure-based methods are often limited by their reliance on static cocrystal structures. To address this, we integrate AlphaFold 2 (AF2) subsampling with a Siamese neural network to predict mutation-induced changes in the relative binding affinity. By leveraging AF2 subsampling, we generated conformational ensembles for Abelson tyrosine kinase (ABL) mutants, shifting the paradigm from single-point predictions to an ensemble-based approach that accounts for intrinsic structural flexibility. Furthermore, we augmented the data set by pairing the generated conformations with reference states, followed by the identification of structurally relevant states via a most probable distribution analysis. To facilitate relative affinity prediction, we developed SIGMA-Net (Siamese structure and graph-aware multistructural affinity prediction network), which was employed to discern features between wild-type and mutant states, enabling free-energy predictions with chemically meaningful accuracy. Benchmarking on the tyrosine kinase inhibitors (TKI) data set and the refined set of PDBbind, our proposed approach achieves higher correlation coefficients for five of six TKI molecules across 31 ABL mutants, outperforming molecular docking and trichannel graph network (TriG-Net). By integrating conformational sampling with Siamese learning, our method enhances both the predictive accuracy and robustness. It achieves absolute binding free energy (ABFE) prediction performance comparable to that of state-of-the-art models such as Boltz-2, whereas Boltz-2 demonstrates better performance in relative binding free energy (RBFE) prediction in the evaluated systems. This framework effectively transcends the limitations of static structure dependence, providing a transferable solution for modeling protein-ligand interactions in highly flexible drug targets.
Following our previous work [Cattin, C. J. Phys. Chem. Lett. 2026, 17(5), 1288-1295], we propose the DMTS-NC approach, a distilled multi-time-step (DMTS) strategy using nonconservative (NC) forces, to further accelerate atomistic molecular dynamics (MD) simulations using foundation neural network models such as FeNNix-Bio1. Therein, a dual-level reversible reference system propagator algorithm (RESPA) formalism couples a target accurate conservative potential to a simplified distilled representation optimized for the production of nonconservative forces. Despite being nonconservative, the distilled architecture is designed to enforce key physical priors, such as equivariance under rotation and cancellation of atomic force components. These choices facilitate the distillation process and therefore improve drastically the robustness of simulation, significantly limiting abnormal discrepancies between the two models, thus achieving excellent agreement with the forces data. Overall, the DMTS-NC scheme is found to be more stable and efficient than its conservative counterpart with additional speedups reaching 15-30% over DMTS. Requiring no fine-tuning steps, it is easier to implement and can be pushed to the limit of the system's physical resonances to maintain accuracy while providing maximum efficiency. We obtain additional speedup by combining hydrogen mass repartitioning (HMR) and High Hydrogen Friction (HHF) to further extend the largest time step up to 10 fs of our schemes while conserving stability and accuracy. As for DMTS, DMTS-NC is applicable to any neural network potential (NNPs) and can be applied to approaches that are computationally heavier than FeNNix-Bio1. We show a proof of principle applying the approach to the distillation of MACE-OFF23 with consequent speedups ranging from 3.66 to 5.64 compared to a single time step.
The level of theory required to predict ultrafast electron diffraction signals is investigated in photoexcited γ-butyrolactone. The total isotropic diffraction signal is calculated both with the independent atom model (IAM) and directly from ab initio electronic wave functions. The results are benchmarked at the equilibrium geometry, along a representative photochemical reaction coordinate (for both ground and excited electronic states), and, most importantly, for a full fewest-switches surface-hopping simulation of the nonadiabatic dynamics. The IAM qualitatively captures the time-resolved UED signal, but there are deviations between the IAM and the ab initio results, irrespective of the electronic structure method used. The errors of the IAM do not get washed out even when the nuclear wavepacket is considered. The errors primarily manifest in the intensity rather than position of scattering features, particularly at small to medium values of momentum transfer, and are sufficient to affect lifetimes and other parameters determined from the diffraction signal. Our results indicate that the IAM can lead to errors when scattering signals are inverted. Crucially, calculations beyond the IAM are necessary to account for the effects of chemical bonding, charge transfer, electronic excitation and ionization.
The Perdew-Zunger self-interaction correction (PZSIC) makes density functional approximations (DFAs) exact for all one-electron densities. However, it overcorrects in many-electron regions, introducing errors for the uniform-density limit, where uncorrected DFAs are exact. The locally scaled PZSIC (LSIC), based on the iso-orbital indicator zσ [which distinguishes single-orbital and slowly varying density regions and is used with the local spin density approximation (LSDA)], restores the uniform-density limit and significantly improves results for many properties, including chemical reaction barrier heights, atomization energies, and ionization potentials. Yet, LSIC performs poorly for weakly bonded systems, leaving many unbound, due to limitations of its iso-orbital indicator. To correct this, in this work we propose a new local scaling, LSIC-α, based on the iso-orbital indicator ασ (which additionally identifies regions of overlapping density tails). A two-parameter scaling function of ασ is fitted to a subset of the nonbonded appropriate norms for the SCAN and r2SCAN meta-GGAs, and tested on many properties of main-group atoms, molecules, and molecular complexes. LSIC-α greatly improves the interaction energies of weakly bonded systems in the S22 data set while retaining LSIC's accuracy for other properties. This work shows that the errors of LSDA (and presumably of higher-level DFAs) can be largely but not entirely repaired by a proper "do no harm" self-interaction correction.