Rydberg excited states of molecules pose a challenge for electronic structure calculations because of their highly diffuse electron distribution. Even large and elaborate atomic basis sets tend to underrepresent the long-range tail, overly confining the Rydberg state. An approach is presented here where the molecular orbitals are variationally optimized for the excited state using a plane wave basis set in a Hartree-Fock calculation, followed by a configuration interaction calculation. The use of excited state optimized orbitals greatly enhances the convergence of the many-body calculation, as illustrated by a full configuration interaction calculation of the 2s Rydberg state of H2. A neural-network-based selective configuration interaction approach is then applied to calculations of 3s and 3p states of H2O and NH3. The obtained values of excitation energy are in close agreement with experimental measurements as well as previous many-body calculations where sufficiently diffuse atomic basis sets were used. Calculations using atomic basis sets lacking extra diffuse functions, such as aug-cc-pVTZ, give significantly higher estimates due to confinement of the Rydberg states.

Gaussian process (GP) regression provides a strategy for accelerating saddle point searches on high-dimensional energy surfaces by reducing the number of times the energy and its derivatives with respect to atomic coordinates need to be evaluated. The computational overhead in the hyperparameter optimization can, however, be large and make the approach inefficient. Failures can also occur if the search ventures too far into regions that are not represented well enough by the GP model. Here, these challenges are resolved by using geometry-aware optimal transport measures and an active pruning strategy using a summation over Wasserstein-1 distances for each atom-type in farthest-point sampling, selecting a fixed-size subset of geometrically diverse configurations to avoid rapidly increasing cost of GP updates as more observations are made. Stability is enhanced by a permutation-invariant metric that provides a reliable trust radius for early-stopping and a logarithmic barrier penalty for the growth of the signal variance. These physically motivated algorithmic changes prove their efficacy by reducing to less than a half the mean computational time on a set of 238 challenging configurations from a previously published data set of chemical reactions. With these improvements, the GP approach is established as a robust and scalable algorithm for accelerating saddle point searches when the evaluation of the energy and atomic forces requires significant computational effort.

Modeling electrocatalytic reactions at solid-liquid interfaces requires capturing both the quantum-mechanical processes at the electrode surface and the complex response of the surrounding electrochemical environment. This review examines the main theoretical frameworks and computational techniques used to describe such systems, focusing on first-principles approaches based on density functional theory (DFT). Key aspects include the treatment of reaction thermodynamics, electrode bias, solvation effects, electrolyte screening, and reaction kinetics. A broad range of methods is discussed, from thermochemical models, such as the computational hydrogen electrode, to potential-dependent formulations based on grand-canonical DFT and explicit calculation of kinetic barriers. The review also highlights recent machine-learning approaches for catalyst screening and the growing use of machine-learning-based force fields, which promise to enable efficient simulations of complex electrochemical environments over extended time and length scales with near-first-principles accuracy. The aim is not only to present the state of the art, but also to clarify the physical assumptions and approximations underlying each approach. The influence of modeling choices on reliability and computational cost is examined in detail. Alongside theoretical aspects, practical considerations are emphasized to support researchers in selecting appropriate methods and designing simulations that are both physically meaningful and computationally tractable.

Accurate determination of the transition states is central to an understanding of reaction kinetics. Double-endpoint methods where both the initial and final states are specified, such as the climbing image nudged elastic band (CI-NEB), identify the minimum energy path between the two and thereby the saddle point on the energy surface that is relevant for the given transition, thus providing an estimate of the transition state within the harmonic transition state theory. Such calculations can, however, incur high computational costs and may suffer stagnation on exceptionally flat or rough energy surfaces. Conversely, methods that only require the specification of an initial set of atomic coordinates, such as the minimum mode following (MMF) method, offer efficiency but can converge on saddle points that are not relevant for the transition of interest. Here, we present an adaptive hybrid algorithm that switches between the CI-NEB and the MMF methods so as to achieve faster convergence to the relevant saddle point. The method is benchmarked On for the Baker–Chan (BC) saddle point test set using the PET-MAD machine-learned potential, along with 59 transitions of a heptamer island on Pt (111) from the OptBench set. A Bayesian analysis of the performance shows a median reduction of energy and force calculations by 57% [95% CrI: 64%, −50%] relative to CI-NEB for the BC set, while a 31% reduction is found for the transitions of the heptamer island. Calculations of the BC set, where a simple switch from the CI-NEB to the MMF method is made when the magnitude of the atomic forces decreases below 0.5 eV/Å, requires 46% more force calculations than the OCI-NEB algorithm. These results show that an adaptive hybrid method mixing CI-NEB and MMF can be a highly efficient tool for high-throughput automated chemical discovery of atomic rearrangements.

Electrolyte solutions at high concentration are indispensable and yet poorly understood. In particular, the extent of speciation─the formation of complexes composed of multiple species─in concentrated ionic solutions is very challenging to obtain theoretically and experimentally, but can have a strong effect on solution properties. The literature is rife with contradictory estimates of speciation from experiments. We find that speciation affects transport properties and is therefore a prerequisite to accurately model concentrated solutions. We turn this to our advantage by showing that the viscosity can be used to determine the extent of complexation in concentrated aqueous solutions. Results of simulations as well as experimental measurements are presented. The atomistic Madrid-2019 force field is extended to model FeCl2. Solutions of FeCl2 and MgCl2 are compared, and the observed difference in viscosity is explained by more complexation in the former, a conclusion supported by recently reported X-ray absorption and neutron scattering experiments.

We present a fully variational locally scaled self-interaction corrected (SIC) energy functional using complex optimal orbitals. This represents an important milestone for fully variational SIC energy functionals, which have been shown to improve the prediction of the properties of atomic, molecular and solid state systems in general, in both ground and excited states. However, it depends on the system and property of the system whether it is beneficial to scale the SIC correction by a factor of one-half, which makes the application of SIC inconsistent. In the limit of a single electron the SIC exactly cancels the self interaction error, but overcorrects the error in regions of high density where there is large overlap between occupied orbitals. The newly implemented local scaling function, z(𝐫), which is based on an iso-orbital indicator derived from considering the kinetic energy density in the iso-electron and many electron case, and takes into account that the orbitals are complex, naturally scales the SIC correction from 0≤ z(𝐫) ≤ 1 in regions of high and low (isolated orbital) electron density. The locally scaled and fully variational SIC framework is general and applicable to atomic, molecular and solid-state systems.

Electrolyte solutions at high concentration are indispensable and yet poorly understood. In particular, the extent of speciation – the formation of complexes composed of multiple species – in concentrated ionic solutions is very challenging to obtain theoretically and experimentally, but can have a strong effect on solution properties. The literature is rife with contradictory estimates of speciation from experiments. We find that speciation affects transport properties, and is therefore, a prerequisite to accurately model concentrated solutions. We turn this to our advantage by showing that the viscosity can be used to determine the extent of complexation in concentrated aqueous solutions. Results of simulations as well as experimental measurements are presented. The atomistic Madrid-2019 force-field is extended to model FeCl_2. Solutions of FeCl_2 and MgCl_2 are compared and the observed difference in viscosity explained by more complexation in the former, a conclusion supported by recently reported X-ray absorption and neutron scattering experiments.

The parameterization of simulation-based models is a central yet laborious task in computational chemistry and physics, often driven by human intuition and manual iteration. Automating this task necessitates the definition of suitable objective functions, which tend to be expensive to evaluate, noisy, non-differentiable, or composed of heterogeneous contributions originating from separate sets of simulations. Gradient-free and black-box optimization algorithms are powerful tools which are particularly well-suited to minimizing such objective functions. Here, we introduce ChemFit, a flexible Python framework for the definition, composition, and massively concurrent evaluation of simulation-based objective functions, which is designed to operate in conjunction with these algorithms. We demonstrate the broad applicability of this approach by using ChemFit for three representative examples of increasing complexity and real-world relevance. First, we obtain the parameters of the Lennard-Jones potential for liquid argon from experimental measurements of the density. Second, we parameterize a polarizable and flexible potential energy function to reproduce the structure of small H_2O clusters obtained from density functional theory calculations. Finally, we tune a small subset of the parameters of a residue-level coarse-grained protein force-field, with the goal to reproduce the experimental critical solution temperature of the low complexity domain of the wild-type hnRNPA1 sequence and an arginine-enriched mutant of this protein. hnRNPA1 is an RNA-binding protein linked to amyotrophic lateral sclerosis. Together, these examples illustrate how ChemFit enables scalable, reproducible, and optimizer-agnostic parameter fitting for broadly applicable multiscale models.

Polaron-mediated charge transport in α-Fe2O3 plays a central role in its performance as a gas-sensing material, yet the atomistic interaction between surface adsorbates and polarons remains insufficiently understood. Here, density functional theory with Hubbard-U correction (DFT+U) combined with nudged elastic band calculations is used to investigate polaron formation, migration, and quenching at the Fe-terminated α-Fe2O3 (0001) surface. The calculated activation energy for small-polaron hopping in bulk α-Fe2O3 is found to be 0.12 eV, in excellent agreement with experimental measurements, confirming the validity of the computational approach. Slab calculations show that migration of the polaron from bulk to the surface lowers the energy by 0.12 eV, indicating preferential localization of charge carriers at the gas-solid interface. Adsorption of NO2 induces substantial electron transfer (0.72 e-) from the oxide to the molecule, eliminating the localized Fe2+ polaron state and thereby suppressing polaronic conductivity. These results provide a direct microscopic explanation for the resistance increase of hematite-based sensors upon exposure to oxidizing gases. More broadly, the study establishes how surface adsorption can modulate charge transport α-Fe2O3 through control of polaron populations, offering design principles for improved iron oxide gas sensors.

Calculations of the lowest valence π* as well as the 3s and higher energy 3pσ Rydberg excited states of the CO2 molecule are carried out using density functionals with variational optimization of the orbitals, an approach involving relatively little computational effort. Five functionals with varying degree of exchange are used in combination with real or complex-valued orbitals that are optimized by finding saddle points on the electronic energy surface corresponding to the excited states. When the PBE functional is used in combination with complex orbitals, the calculated excitation energy is found to be within 0.3 eV of multireference configuration interaction reference values, and the results are further improved with hybrid functionals. In contrast, linear-response time-dependent density functional theory calculations give errors up to 1.9 eV for the most diffuse 3pσ excitation and exhibit stronger dependence on both the excitation character and the functional used. Calculated C-O dissociation curves using the PBE functional and the orbital-optimized approach compare remarkably well with the reported multireference configuration interaction and equation-of-motion coupled-cluster singles and doubles calculations. Thanks to the low computational cost, these results demonstrate that orbital-optimized density functional calculations can be a promising route for modelling photorelaxation in condensed-phase CO2, for example in the context of interstellar cosmic-ray radiation driven process involving high-energy Rydberg states.

A general polarizable embedded (PE) quantum mechanics/molecular mechanics scheme for periodic systems is presented, describing mutual polarization of the two subsystems. The QM system, described with density functional theory (DFT), is coupled to a single center multipole expansion (SCME) model, characterising H_2O molecules in the MM region. In SCME the H_2O molecules are ascribed anisotropic dipole and quadrupole polarizabilities and permanent multipoles up to and including the hexadecapole. Our embedding scheme illustrates a smooth and efficient convergence pattern of the periodic interaction potential by introducing a single and clustered multipole expansion points in the far-field. By choosing the near- and far-field expansion of the potential carefully the PE-QM/MM calculation matches the level of accuracy of a the QM calculation. In the short range, the electrostatic interaction between the QM and MM subsystems is damped with a real-space and pair-wise isotropic damping functions - resulting in a screened interaction and preventing over-polarization. In molecular dynamics simulations the two subsystems are separated with the elastic scattering assisted flexible inner region [Kirchhoff et. al. JCTC, 2021, 17, 9, 5863] - ensuring a smooth transition in the radial distribution at the boundary between the two subsystems.

The excited electronic states involved in the optical cycle preparation of a pure spin state of the negatively charged NV-defect in diamond are calculated using the HSE06 hybrid density functional and variational optimization of the orbitals. This includes the energy of the excited triplet as well as the two lowest singlet states with respect to the ground triplet state. In addition to the vertical excitation, the effect of structural relaxation is also estimated using analytical atomic forces. The lowering of the energy in the triplet excited state and the resulting zero-phonon line triplet excitation energy are both within 0.1 eV of the experimental estimates. An analogous relaxation in the lower energy singlet state using spin purified atomic forces is estimated to be 0.06 eV. These results, obtained with a hybrid density functional, improve on previously published results using local and semi-local functionals, which are known to underestimate the band gap. The good agreement with experimental estimates demonstrates how time-independent variational calculations of excited states using density functionals can give accurate results and, thereby, provide a powerful screening tool for identifying other defect systems as candidates for quantum technologies.

Calculations based on density functional theory and statistical mechanics reveal how Fe site preference in olivine is altered at interfaces. Although the M1 site is favored in bulk olivine, surface metal sites provide greater stabilisation for high-spin Fe2+. This enrichment accounts for the enhanced reactivity of olivine interfaces towards dissolution, carbonation and catalysis.

Copper‐based catalysts are of particular interest for electrochemical reduction of (CO2RR) as products beyond CO can form. To improve activity and selectivity, several studies have focused on the addition of other elements as substitutional impurities. Although the adsorption of a single CO molecule has often been used as a descriptor for CO2RR activity, our recent calculations using the RPBE functional showed that multiple CO molecules can bind to first‐row transition metal impurities. Here, we extend the study to second‐row transition metals and also to a functional that explicitly includes dispersion interaction, BEEF‐vdW. The binding energy of the first CO molecule on the impurity atom is found to be significantly larger than on the clean Cu(111) and Cu(100) surfaces, but the differential binding energy generally drops as more CO molecules adsorb. The dispersion interaction is found to make a significant contribution to the binding energy, in particular for the last and weakest bound CO molecule, the one that is most likely to participate in CO2RR. In some cases, four CO admolecules can bind more strongly on the impurity atom than on the clean copper surface. The adsorption of CO causes the position of the impurity atom to shift outwards and in some cases, even escape from the surface layer. The C─O stretch frequencies are calculated in order to identify possible experimental signatures of multiple CO adsorption.

A method for time-reversible numerical integration of the deterministic Landau-Lifshitz Gilbert equation by means of a second order Suzuki-Trotter decomposition is presented and tested against commonly used second order predictor-corrector methods. We find the time-reversibility of the Suzuki-Trotter integrator to be superior by several orders of magnitude while the computational effort is similar. Calculations of trajectories backwards in time are useful, for example, when evaluating dynamical corrections to transition state theory.

An efficient and scalable implementation of a method for locating first-order saddle points on the energy surface of a magnetic system is presented, along with several applications in which the mechanisms of various magnetic transitions are identified. The starting point for the iterative search algorithm can be anywhere, even close to a local energy minimum representing an initial state of the system, and the final state need not be specified. Convergence on a saddle point is obtained by inverting the component of the gradient along the minimum mode, thereby effectively transforming the neighborhood of the saddle point to that of a local minimum. The method requires only the lowest two eigenvalues and corresponding eigenvectors of the Hessian of the system's energy and they are found using a quasi-Newton limited-memory Broyden-Fletcher-Goldfarb-Shanno solver for the minimization of the Rayleigh quotient without explicit evaluation of the Hessian. The method is applicable to large systems, as it does not introduce additional scaling overhead to the computational complexity determined by the interactions present in the system. Applications are presented to transitions in systems that reveal significant complexity of coexisting magnetic states, such as skyrmions, skyrmion bags, skyrmion tubes, chiral bobbers, and globules. The identification of new metastable three-dimensional (3D) textures, such as magnetic bobbers with extended equilibrium distance between the base and the terminating Bloch point, and magnetic globules appearing as isolated states in 3D due to magnetostatic interactions, demonstrates the usefulness of the method for the characterization of complex energy surfaces of magnetic systems. When combined with rate theory within the harmonic approximation, the method can be used for simulations of the long timescale dynamics of complex magnetic systems characterized by multiple metastable states.

The effect of substituting a hydrogen atom by a chlorine atom or a methyl group by a trichloromethyl (CCl3) group at the stereogenic center of light-driven second genera- tion molecular motors is calculated in order to assess the effect on rotational speed and the separation of the absorption peaks of the isomers. While experimental and theoret- ical studies have previously been carried out for fluorine substitution, this is the first study of chlorine substitution. Five well-characterized base molecules are studied and the trends are compared with the effect of fluorine substitution. The trichloromethyl substitution is found to accelerate the rotation more than a trifluoromethyl (CF3) sub- stitution by reducing the life-time of the metastable state, due to larger steric hindrance in the metastable state than in the transition state for the thermal helix inversion (THI). A larger increase in the separation of the absorption peaks of the two isomers is also ob- tained. The Cl atom substitution, however, changes the energy landscape significantly, making the M isomer lower in energy than the P isomer, and raising the energy barrier for THI beyond that of the back transition, thus quenching the rotation.

The task of locating first order saddle points on high-dimensional surfaces describing the variation of energy as a function of atomic coordinates is an essential step for identifying the mechanism and estimating the rate of thermally activated events within the harmonic approximation of transition state theory. When combined directly with electronic structure calculations, the number of energy and atomic force evaluations needed for convergence is a primary issue. Here, we describe an efficient implementation of Gaussian process regression (GPR) acceleration of the minimum mode following method where a dimer is used to estimate the lowest eigenmode of the Hessian. A surrogate energy surface is constructed and updated after each electronic structure calculation. The method is applied to a test set of 500 molecular reactions previously generated by Hermes, E. D. [ J. Chem. Theory Comput. 2022, 18, 6974-6988]. An order of magnitude reduction in the number of electronic structure calculations needed to reach the saddle point configurations is obtained by using the GPR. Despite the wide range in stiffness of the molecular degrees of freedom, the calculations are carried out using Cartesian coordinates and are found to require similar number of electronic structure calculations as an elaborate internal coordinate method implemented in the Sella software package. The present implementation of the GPR surrogate model in C++ is efficient enough for the wall time of the saddle point searches to be reduced in 3 out of 4 cases even though the calculations are carried out at a low Hartree-Fock level.

The switching mechanisms in artificial spin ice systems are investigated with focus on shakti and modified shakti lattices. Minimum energy paths are calculated using the geodesic nudged elastic band (GNEB) method implemented with a micromagnetic description of the system, including the internal magnetic structure of the islands and edge modulations. Two switching mechanisms, uniform magnetization rotation and domain wall formation, are found to have comparable activation energy. The preference for one over the other depends strongly on the saturation magnetization and the magnetic ordering of neighboring islands. Surprisingly, these mechanisms can coexist, leading to an enhanced probability of magnetization reversal. These results provide valuable insight that can help control internal magnetization switching processes in spin ice systems and help predict their thermodynamic properties.