Analytical gradients for nuclear–electronic orbital multistate density functional theory: Geometry optimizations and reaction paths

Proton-coupled electron transfer (PCET) reactions on semiconducting metal oxide surfaces often involve charged defects in the form of electron or hole polarons. Herein, vibronically nonadiabatic PCET theory is used to model rate constants for the PCET reaction between a reduced anatase TiO2(101) surface and 4-MeO-TEMPO, where electron polarons on the TiO2 surface directly participate in the PCET reaction. This modeling strategy treats the transferring proton as well as all electrons quantum mechanically and includes the effects of excited vibronic states. The rate constant expression depends on the reorganization energy, as well as the reaction free energies and vibronic couplings for different pairs of vibronic states, and accounts for proton donor–acceptor motion. Hybrid functional periodic density functional theory (DFT) is used to calculate the parameters in the rate constant expression, and a Hubbard α-based constrained DFT approach is used to enforce charge constraints consistent with the two electronically diabatic states for the PCET reaction. This modeling strategy is applied to compute the PCET rate constants and kinetic isotope effects for reactions involving five-coordinate and six-coordinate Ti3+ defects on the TiO2(101) surface, showing that excited vibronic states contribute significantly to the rate constant for both defects, especially for deuterium. This study highlights the importance of hydrogen tunneling and excited vibronic states in interfacial PCET reactions. Such modeling strategies can be used to further understand and tailor the reactivity of metal oxide surfaces for energy conversion.

The incorporation of nuclear quantum effects and non-Born-Oppenheimer behavior into quantum chemistry calculations and molecular dynamics simulations is a longstanding challenge. The nuclear-electronic orbital (NEO) approach treats specified nuclei, typically protons, quantum mechanically on the same level as the electrons with wave function and density functional theory methods. This approach inherently includes nuclear delocalization and zero-point energy in molecular energy calculations, geometry optimizations, reaction paths, and dynamics. It can also provide accurate descriptions of excited electronic, vibrational, and vibronic states as well as nuclear tunneling and nonadiabatic dynamics. Nonequilibrium nuclear-electronic dynamics simulations beyond the Born-Oppenheimer approximation can be used to investigate a wide range of excited state processes. This Perspective provides an overview of the foundational NEO methods and enumerates the prospects for using these methods as building blocks for future developments. The conceptual simplicity and computational efficiency of the NEO approach will enhance its accessibility and applicability to diverse chemical and biological systems.

Geometries and energies of various metal aqua ions have been computed with Hartree-Fock (HF), single-reference second-order perturbation theory according to Møller and Plesset (MP2), complete active space self-consistent field (CAS-SCF) methods, multi-configurational quasi-degenerate second-order perturbation theory (MCQDPT2), and density functional theory (DFT), whereby for the latter the most widely used functionals BLYP, B3LYP, LDA ("local density approximation"), and another functional (SOP) exhibiting Slater exchange, have been investigated. The geometries as well as the energies are sensitive to the computational method, results of which can lead to different conclusions. The geometries of certain reactants, products, intermediates, or transition states cannot be obtained with all of the above-mentioned computational methods: in some cases, the structures are very poor, exhibiting large errors in the bond parameters, and in others, the targeted species cannot be calculated at all. The different computational techniques do not always predict the same coordination number of a given aqua ion or the same preferred water exchange mechanism. Inaccuracies arising from inadequate basis sets are also discussed.

Two methods are developed to enhance the stability, efficiency, and robustness of reaction path optimization using a chain of replicas. First, distances between replicas are kept equal during path optimization via holonomic constraints. Finding a reaction path is, thus, transformed into a constrained optimization problem. This approach avoids force projections for finding minimum energy paths (MEPs), and fast-converging schemes such as quasi-Newton methods can be readily applied. Second, we define a new objective function - the total Hamiltonian - for reaction path optimization, by combining the kinetic energy potential of each replica with its potential energy function. Minimizing the total Hamiltonian of a chain determines a minimum Hamiltonian path (MHP). If the distances between replicas are kept equal and a consistent force constant is used, then the kinetic energy potentials of all replicas have the same value. The MHP in this case is the most probable isokinetic path. Our results indicate that low-temperature kinetic energy potentials (<5 K) can be used to prevent the development of kinks during path optimization and can significantly reduce the required steps of minimization by 2-3 times without causing noticeable differences between a MHP and MEP. These methods are applied to three test cases, the C7eq-to-Cax isomerization of an alanine dipeptide, the (4)C1-to-(1)C4 transition of an α-d-glucopyranose, and the helix-to-sheet transition of a GNNQQNY heptapeptide. By applying the methods developed in this work, convergence of reaction path optimization can be achieved for these complex transitions, involving full atomic details and a large number of replicas (>100). For the case of helix-to-sheet transition, we identify pathways whose energy barriers are consistent with experimental measurements. Further, we develop a method based on the work energy theorem to quantify the accuracy of reaction paths and to determine whether the atoms used to define a path are enough to provide quantitative estimation of energy barriers.

Exploring the potential energy surface of the ethyl cation by new procedures

Hydrogen atom transfer reactions between the substrate and coenzyme are key mechanistic features of all adenosylcobalamin-dependent enzymes. For one of these enzymes, glutamate mutase, we have investigated whether hydrogen tunneling makes a significant contribution to the mechanism by examining the temperature dependence of the deuterium kinetic isotope effect associated with the transfer of a hydrogen atom from methylaspartate to the coenzyme. To do this, we designed a novel intramolecular competition experiment that allowed us to measure the intrinsic kinetic isotope effect, even though hydrogen transfer may not be rate-determining. From the Arrhenius plot of the kinetic isotope effect, the ratio of the pre-exponential factors (A(H)/A(D)) was 0.17 +/- 0.04 and the isotope effect on the activation energy [DeltaE(a(D-H))] was 1.94 +/- 0.13 kcal/mol. The results imply that a significant degree of hydrogen tunneling occurs in glutamate mutase, even though the intrinsic kinetic isotope effects are well within the semiclassical limit and are much smaller than those measured for other AdoCbl enzymes and model reactions for which hydrogen tunneling has been implicated.

Electronic structure calculations at the correlated molecular orbital theory and density functional theory levels have been used to generate a reliable set of clustering energies for up to three water molecules in carbon dioxide clusters up to n = 12. The structures and energetics are dominated by Lewis acid-base interactions with hydrogen-bonding interactions playing a lesser energetic role. The actual binding energies are somewhat larger than might be expected. The correlated molecular orbital MP2 method and density functional theory with the ωB97X exchange-correlation functional provide good results for the energetics of the clusters, but the B3LYP and ωB97X-D functionals do not. Seven CO(2) molecules form the first solvent shell about a single H(2)O with four CO(2) molecules interacting with the H(2)O via Lewis acid-base interactions, two CO(2) interacting with the H(2)O by hydrogen bonds, and the seventh CO(2) completing the shell. The Lewis acid-base and weak hydrogen bond interactions between the water molecules and the CO(2) molecules are strong enough to disrupt the trimer ring configuration for as few as seven CO(2) molecules. Calculated (13)C NMR chemical shifts for mCO(2)·(H(2)O)(n) show little change with respect to the number of H(2)O or CO(2) molecules in the cluster. The O-H stretching frequencies do exhibit shifts that can provide information about the interactions between water and CO(2) molecules.

The gas phase electron resonance spectra of NCO in its ground 2Π3/2 vibronic state and in two excited vibronic states are described. Theoretical analysis of the spectra yields effective g values for the three states. In additon the 14N magnetic hyperfine and electric quadrupole coupling constants and the electric dipole moment are determined. The theory of the Renner coupling of electronic and vibrational motion is extended, and shown to account for anomalous contributions to the g values. The theory also shows that these contributions are closely related to the Renner coupling constant.

Recently (Photochem Photobiol. 2023;100:1277-1289. doi:10.1111/php.13898), we described the anti-Kasha effect in tribenzo-6H-1,4-diazepinoporphyrazins with C2v symmetry, where the ultrafast spin changes successfully compete with the internal conversion. In this study, we show the presence of this effect in 2 (3),9 (10),16(17),23(24)-tetra-tert-butyl-29H,31H-phthalocyanine (1) and 1,4-di-[2-(2-methoxyethoxy)ethoxy]-29H,31H-phthalocyanine (2), which also possess reduced molecular symmetry and do not bear 6H-1,4-diazepine fragments. The anti-Kasha effect in 1 and 2 supplemented by Mg(II) tribenzo-6H-1,4-diazepinoporphyrazinates 3 and 4 exhibits a close-to-linear dependence on energy gap value between the zero vibrational levels of two lowest singlet excited states S1 0 and S2 0 (these states are degenerate in D4h symmetry) and enhances with increase. The theoretical kinetic model of excited state dynamics, which takes into account the observed effects and follows Fermi's golden rule, predicts the presence of an additional excited state with enhanced spin-orbit coupling compared to S1 0, S2 0 and the corresponding triplet states, which is not predicted by TDDFT calculations in the Born-Oppenheimer approximation. The combination of the above indicates that the key role in the observed anti-Kasha effect and the mechanism of dissipation of the excited state in porphyrazines and their analogs is played by vibronic excited states, which requires theoretical research methods beyond the Born-Oppenheimer approximation.

Review: 8163 refs.].

We have firstly developed the ab initio linear combination of gaussian type orbital (LCGTO) program of generalized spin orbital (GSO) density functional theory (DFT) and applied it for multicenter metal clusters with noncollinear spin structures. It was found that two- and three-dimensional spin states becomes ground-states for Cr 3 (D 3h ) and CR 4 (T d ) respectively, indicating importance of GSO treatment of DFT for these clusters. Generalized spin orbital density functional theory multicenter metal clusters noncollinear molecular magnetism

Model calculations of kinetic isotope effects for the reactions: C2H6+·CH3 → ·C2H5+CH4 CH4 + ·CF3 → ·CH3+CHF3 are reported. Transition state geometries were those calculated by Dewar and coworkers using the MNDO semiempirical method. Transition state force fields were formulated from empirical expressions for stretching, bending, linear bending and interaction valence force constants by using bond orders as disposable parameters. Although it proved impossible to assign a best force field to either reaction, the calculated isotope effects generally were in satisfactory agreement with experiment, and were better than those calculated from the MNDO potential energy surface. Hydrogen tunnelling is apparently implicated in the CH4+ CF3 reaction.

The aim of this paper is to prove that safe success in finding reaction paths (RPs) can only be expected from global path-determining methods. Some extensions of the mathematical arguments leading to the introduction of the DDRP (dynamically defined reaction path) method have been sketched. Four cases involving relaxation of analyticity, variability of the gradient field, minimum energy (reaction) paths (MEPs) and “golf pocket holes” on the potential energy surface (PES), and the rather strange consequences of the main theorem of the DDRP method giving a rigorous mathematical basis to chemical intuition in reaction kinetics have been discussed. The discussions show that the DDRP method - when changing the conditions and parameters - may, in essence, involve all other global methods. It has been shown that the DDRP method works in a stable way even for non-analytic though smooth energy functions; moreover, the gradient field can be replaced by other vector fields resulting in better convergence to the reaction path. As a by-product, the question of the existence of MEPs can safely be handled and golf pocket holes are constructed on the PES in order to prove that local methods have chance to search faithfully the RPs in complicated systems only if the energy function can be restored from its arbitrarily small pieces.

The elementary vapor-phase reaction between Cl atoms and HCHO has been studied by ab initio methods. Calculations at the MP2, MP3, MP4(SDTQ), CCSD, CCSD(T), and MRD-CI levels of theory show that the reaction is characterized by a low electronic barrier; excluding the effects of spin-orbit splitting in Cl, our best estimate at the MRD-CI/aug-cc-pVTZ//RHF-RCCSD(T)/aug-cc-pVTZ level of theory predicts a Born-Oppenheimer barrier height of 0.7 kJ mol-1. The energies of the lowest six electronic states as resulting from MRD-CI calculations are presented at discrete points along the reaction path, and two avoided crossings are found in the transition state region. The spin-orbit splitting in Cl is also calculated along the reaction path; it is not negligible in the transition state region and is found to increase the barrier by only 1.4 kJ mol-1 at the RCCSD(T)/aug-cc-pVTZ transition state geometry. The minimum energy path of the reaction connects an energetically weakly stabilized adduct on the flat potential surface on the reactant side and an energetically strongly stabilized postreaction adduct. The reaction rate coefficient and the kinetic isotope effects were calculated using improved canonical variational theory with small curvature tunneling (ICVT/SCT), and the results were compared to experimental data. The experimental reaction rate coefficient is reproduced within its uncertainty limits by variational transition state theory with interpolated single-point energy corrections (ISPE) at the MP4(SDTQ) level of theory and by conventional transition state theory with interpolated optimized energies (IOE) at the MRD-CI//RCCSD(T) level of theory and interpolated optimized geometries at the RCCSD(T) level of theory on an MP2/aug-cc-pVTZ potential energy surface when employing scaled vibrational frequencies.

AtomREM: Non-empirical seeker of the minimum energy escape paths on many-dimensional potential landscapes without coarse graining