The long road to calibrated prediction uncertainty in computational chemistry.

Uncertainty quantification (UQ) in computational chemistry (CC) is still in its infancy. Very few CC methods are designed to provide a confidence level on their predictions, and most users still rely improperly on the mean absolute error as an accuracy metric. The development of reliable UQ methods is essential, notably for CC to be used confidently in industrial processes. A review of the CC-UQ literature shows that there is no common standard procedure to report or validate prediction uncertainty. I consider here analysis tools using concepts (calibration and sharpness) developed in meteorology and machine learning for the validation of probabilistic forecasters. These tools are adapted to CC-UQ and applied to datasets of prediction uncertainties provided by composite methods, Bayesian ensembles methods, and machine learning and a posteriori statistical methods.

Prediction uncertainty validation for computational chemists.

Validation of prediction uncertainty (PU) is becoming an essential task for modern computational chemistry. Designed to quantify the reliability of predictions in meteorology, the calibration-sharpness (CS) framework is now widely used to optimize and validate uncertainty-aware machine learning (ML) methods. However, its application is not limited to ML and it can serve as a principled framework for any PU validation. The present article is intended as a step-by-step introduction to the concepts and techniques of PU validation in the CS framework, adapted to the specifics of computational chemistry. The presented methods range from elementary graphical checks to more sophisticated ones based on local calibration statistics. The concept of tightness, is introduced. The methods are illustrated on synthetic datasets and applied to uncertainty quantification data issued from the computational chemistry literature.

One Way Traffic: Base-to-Backbone Hole Transfer in Nucleoside Phosphorodithioate.

The directionality of the hole-transfer processes between DNA backbone and base was investigated by using phosphorodithioate [P(S- )=S] components. ESR spectroscopy in homogeneous frozen aqueous solutions and pulse radiolysis in aqueous solution at ambient temperature confirmed initial formation of G.+ -P(S- )=S. The ionization potential of G-P(S- )=S was calculated to be slightly lower than that of guanine in 5'-dGMP. Subsequent thermally activated hole transfer from G.+ to P(S- )=S led to dithiyl radical (P-2S. ) formation on the µs timescale. In parallel, ESR spectroscopy, pulse radiolysis, and density functional theory (DFT) calculations confirmed P-2S. formation in an abasic phosphorodithioate model compound. ESR investigations at low temperatures and higher G-P(S- )=S concentrations showed a bimolecular conversion of P-2S. to the σ2 -σ*1 -bonded dimer anion radical [-P-2S - . 2S-P-]- [ΔG (150 K, DFT)=-7.2 kcal mol-1 ]. However, [-P-2S - . 2S-P-]- formation was not observed by pulse radiolysis [ΔG° (298 K, DFT)=-1.4 kcal mol-1 ]. Neither P-2S. nor [-P-2S - . 2S-P-]- oxidized guanine base; only base-to-backbone hole transfer occurs in phosphorodithioate.

One Way Traffic: Base-to-Backbone Hole Transfer in Nucleoside Phosphorodithioate.

Invited for the cover of this issue are the groups of Roman Dembinski, Mehran Mostafavi, and Amitava Adhikary at the Polish Academy of Sciences, Université Paris-Saclay, and Oakland University. The image depicts a doughnut as a way of illustrating the hole transfer process. Read the full text of the article at 10.1002/chem.202000247.

Probabilistic performance estimators for computational chemistry methods: Systematic improvement probability and ranking probability matrix. I. Theory.

The comparison of benchmark error sets is an essential tool for the evaluation of theories in computational chemistry. The standard ranking of methods by their mean unsigned error is unsatisfactory for several reasons linked to the non-normality of the error distributions and the presence of underlying trends. Complementary statistics have recently been proposed to palliate such deficiencies, such as quantiles of the absolute error distribution or the mean prediction uncertainty. We introduce here a new score, the systematic improvement probability, based on the direct system-wise comparison of absolute errors. Independent of the chosen scoring rule, the uncertainty of the statistics due to the incompleteness of the benchmark datasets is also generally overlooked. However, this uncertainty is essential to appreciate the robustness of rankings. In the present article, we develop two indicators based on robust statistics to address this problem: Pinv, the inversion probability between two values of a statistic, and Pr, the ranking probability matrix. We demonstrate also the essential contribution of the correlations between error sets in these scores comparisons.

Probabilistic performance estimators for computational chemistry methods: Systematic improvement probability and ranking probability matrix. II. Applications.

In Paper I [P. Pernot and A. Savin, J. Chem. Phys. 152, 164108 (2020)], we introduced the systematic improvement probability as a tool to assess the level of improvement on absolute errors to be expected when switching between two computational chemistry methods. We also developed two indicators based on robust statistics to address the uncertainty of ranking in computational chemistry benchmarks: Pinv, the inversion probability between two values of a statistic, and Pr, the ranking probability matrix. In this second part, these indicators are applied to nine data sets extracted from the recent benchmarking literature. We also illustrate how the correlation between the error sets might contain useful information on the benchmark dataset quality, notably when experimental data are used as reference.

Direct observation of the oxidation of DNA bases by phosphate radicals formed under radiation: a model of the backbone-to-base hole transfer.

In irradiated DNA, by the base-to-base and backbone-to-base hole transfer processes, the hole (i.e., the unpaired spin) localizes on the most electropositive base, guanine. Phosphate radicals formed via ionization events in the DNA-backbone must play an important role in the backbone-to-base hole transfer process. However, earlier studies on irradiated hydrated DNA, on irradiated DNA-models in frozen aqueous solution and in neat dimethyl phosphate showed the formation of carbon-centered radicals and not phosphate radicals. Therefore, to model the backbone-to-base hole transfer process, we report picosecond pulse radiolysis studies of the reactions between H2PO4˙ with the DNA bases - G, A, T, and C in 6 M H3PO4 at 22 °C. The time-resolved observations show that in 6 M H3PO4, H2PO4˙ causes the one-electron oxidation of adenine, guanine and thymine, by forming the cation radicals via a single electron transfer (SET) process; however, the rate constant of the reaction of H2PO4˙ with cytosine is too low (<107 L mol-1 s-1) to be measured. The rates of these reactions are influenced by the protonation states and the reorganization energies of the base radicals and of the phosphate radical in 6 M H3PO4.

Probabilistic performance estimators for computational chemistry methods: The empirical cumulative distribution function of absolute errors.

Benchmarking studies in computational chemistry use reference datasets to assess the accuracy of a method through error statistics. The commonly used error statistics, such as the mean signed and mean unsigned errors, do not inform end-users on the expected amplitude of prediction errors attached to these methods. We show that, the distributions of model errors being neither normal nor zero-centered, these error statistics cannot be used to infer prediction error probabilities. To overcome this limitation, we advocate for the use of more informative statistics, based on the empirical cumulative distribution function of unsigned errors, namely, (1) the probability for a new calculation to have an absolute error below a chosen threshold and (2) the maximal amplitude of errors one can expect with a chosen high confidence level. Those statistics are also shown to be well suited for benchmarking and ranking studies. Moreover, the standard error on all benchmarking statistics depends on the size of the reference dataset. Systematic publication of these standard errors would be very helpful to assess the statistical reliability of benchmarking conclusions.

Effect of the solvation state of electron in dissociative electron attachment reaction in aqueous solutions.

It is generally considered that the pre-solvated electron and the solvated electron reacting with a solute yield the same product. Silver cyanide complex, Ag(CN)2-, is used as a simple probe to demonstrate unambiguously the existence of a different reduction mechanism for pre-hydrated electrons. Using systematic multichannel transient absorption measurements at different solute concentrations from millimolar to decimolar, global data analysis and theoretical calculations, we present the dissociative electron attachment on Ag(CN)2-. The short-lived silver complex, Ag0(CN)22-, formed by hydrated electron with nanosecond pulse radiolysis, can be observed at room temperature. However, at higher temperatures only the free silver atom, Ag0, is detected, suggesting that Ag0(CN)22- dissociation is fast. Surprisingly, pulse radiolysis measurements on Ag(CN)2- reduction, performed by a 7 ps electron pulse at room temperature, show clearly that a new reduced form of silver complex, AgCN-, is produced within the pulse. This species, absorbing at 560 nm, is not formed by the hydrated electron but exclusively by its precursor. DFT calculations show that the different reactivity of the hydrated and pre-hydrated electrons can be due to the formation of different electronic states of Ag0(CN)22-: the prehydrated electron can form an excited state of this complex, which mainly dissociates into Ag0CN- + CN-.

The parameter uncertainty inflation fallacy.

Statistical estimation of the prediction uncertainty of physical models is typically hindered by the inadequacy of these models due to various approximations they are built upon. The prediction errors caused by model inadequacy can be handled either by correcting the model's results or by adapting the model's parameter uncertainty to generate prediction uncertainties representative, in a way to be defined, of model inadequacy errors. The main advantage of the latter approach (thereafter called PUI, for Parameter Uncertainty Inflation) is its transferability to the prediction of other quantities of interest based on the same parameters. A critical review of implementations of PUI in several areas of computational chemistry shows that it is biased, in the sense that it does not produce prediction uncertainty bands conforming to model inadequacy errors.

Energetics of Photoinduced Charge Migration within the Tryptophan Tetrad of an Animal (6-4) Photolyase.

Cryptochromes and photolyases are flavoproteins that undergo cascades of electron/hole transfers after excitation of the flavin cofactor. It was recently discovered that animal (6-4) photolyases, as well as animal cryptochromes, feature a chain of four tryptophan residues, while other members of the family contain merely a tryptophan triad. Transient absorption spectroscopy measurements on Xenopus laevis (6-4) photolyase have shown that the fourth residue is effectively involved in photoreduction but at the same time could not unequivocally ascertain the final redox state of this residue. In this article, polarizable molecular dynamics simulations and constrained density functional theory calculations are carried out to reveal the energetics of charge migration along the tryptophan tetrad. Migration toward the fourth tryptophan is found to be thermodynamically favorable. Electron transfer mechanisms are sought either through an incoherent hopping mechanism or through a multiple sites tunneling process. The Jortner-Bixon formulation of electron transfer (ET) theory is employed to characterize the hopping mechanism. The interplay between electron transfer and relaxation of protein and solvent is analyzed in detail. Our simulations confirm that ET in (6-4) photolyase proceeds out of equilibrium. Multiple site tunneling is modeled with the recently proposed flickering resonance mechanism. Given the position of energy levels and the distribution of electronic coupling values, tunneling over three tryptophan residues may become competitive in some cases, although a hopping mechanism is likely to be the dominant channel. For both reactive channels, computed rates are very sensitive to the starting protein configuration, suggesting that both can take place and eventually be mixed, depending on the state of the system when photoexcitation takes place.

Molecular Isomer Identification of Titan's Tholins Organic Aerosols by Photoelectron/Photoion Coincidence Spectroscopy Coupled to VUV Synchrotron Radiation.

The chemical composition of Titan organic haze is poorly known. To address this issue, laboratory analogues named tholins are synthesized and analyzed by methods often requiring an extraction process in a carrier solvent. These methods exclude the analysis of the insoluble tholins' fraction and assume a hypothetical chemical equivalence between soluble and insoluble fractions. In this work, we present a powerful complementary analysis method recently developed on the DESIRS VUV synchrotron beamline at SOLEIL. It involves soft pyrolysis of tholins at ∼230 °C and electron/ion coincidence analysis of the emitted volatile compounds photoionized by tunable synchrotron radiation. By comparison with reference photoelectron spectra (PES), the spectral information collected on the detected molecules yields their isomeric structure. The method is more readily applied to light species (m/z ≤ 69), while for heavier ones, the number of possibilities and the lack of PES reference spectra in the literature limit its analysis. A notable pattern in the analyzed tholins is the presence of species containing adjacent doubly bonded N atoms, which might be a signature of heterogeneous incorporation of N2 in tholins.

Prediction uncertainty of density functional approximations for properties of crystals with cubic symmetry.

The performance of a method is generally measured by an assessment of the errors between the method's results and a set of reference data. The prediction uncertainty is a measure of the confidence that can be attached to a method's prediction. Its estimation is based on the random part of the errors not explained by reference data uncertainty, which implies an evaluation of the systematic component(s) of the errors. As the predictions of most density functional approximations (DFA) present systematic errors, the standard performance statistics, such as the mean of the absolute errors (MAE or MUE), cannot be directly used to infer prediction uncertainty. We investigate here an a posteriori calibration method to estimate the prediction uncertainty of DFAs for properties of solids. A linear model is shown to be adequate to address the systematic trend in the errors. The applicability of this approach to modest-size reference sets (28 systems) is evaluated for the prediction of band gaps, bulk moduli, and lattice constants with a wide panel of DFAs.

Calibration of forcefields for molecular simulation: sequential design of computer experiments for building cost-efficient kriging metamodels.

We present a global strategy for molecular simulation forcefield optimization, using recent advances in Efficient Global Optimization algorithms. During the course of the optimization process, probabilistic kriging metamodels are used, that predict molecular simulation results for a given set of forcefield parameter values. This enables a thorough investigation of parameter space, and a global search for the minimum of a score function by properly integrating relevant uncertainty sources. Additional information about the forcefield parameters are obtained that are inaccessible with standard optimization strategies. In particular, uncertainty on the optimal forcefield parameters can be estimated, and transferred to simulation predictions. This global optimization strategy is benchmarked on the TIP4P water model.

Oxidation of bromide ions by hydroxyl radicals: spectral characterization of the intermediate BrOH•-.

The reaction of (•)OH with Br(-) has been reinvestigated by picosecond pulse radiolysis combined with streak camera absorption detection and the obtained spectro-kinetics data have been globally analyzed using Bayesian data analysis. For the first time, the absorption spectrum of the intermediate species BrOH(•-) has been determined. This species absorbs in the same spectral domain as Br(2)(•-): the band maximum is roughly at the same wavelength (λ(max) = 352 nm instead of 354 nm) but the extinction coefficient is smaller (ε(max) = 7800 ± 400 dm(3) mol(-1) cm(-1) compared with 9600 ± 300 dm(3) mol(-1) cm(-1)) and the band is broader (88 nm versus 76 nm). Quantum chemical calculations have also been performed and corroborate the experimental results. In contrast to Br(2)(•-), the existence of several water-BrOH(•-) configurations leading to different transition energies may account for the broadening of the absorption spectrum in addition to the higher number of degrees of freedom.

Reduction of earth alkaline metal salts in THF solution studied by picosecond pulse radiolysis.

Picosecond pulse radiolysis of tetrahydrofuran (THF) solutions containing earth alkaline metal salt, M(II)(ClO4)2, at different concentrations are performed using two different supercontinua as probe pulse, one covering the visible and another the near-infrared (NIR) down to the visible. Two types of line scan detectors are used to record the absorption spectra in the range from 400 to 1500 nm. Because of the strong overlap between the spectra of the absorbing species in the present wavelength range, global matrices were built for each M(II) system, by delay-wise binding the matrix for pure THF with the available matrices for this cation. The number of absorbers was assessed by Singular Value Decomposition of the global matrix, and a MCR-ALS analysis with the corresponding number of species was performed. The analysis of the results show clearly that solvated electron reacts with the earth alkaline metal molecule and the product has an optical absorption band very different than that of solvated electron in pure THF. So, contrarily to the case of solution containing free Na(+), in the presence of Mg(II), Ca(II) and Sr(II) the observed absorption band is not only blueshifted, but its shape is also drastically changed. In fact with Na(+) solvated electron forms a tight-contact pair but with earth alkaline metal cation solvated electron is scavenged by the undissociated molecule M(II)(ClO4)2. In order to determine the structure of the absorbing species observed after the electron pulse, Monte Carlo/DFT simulations were performed in the case of Mg(II), based on a classical Monte Carlo code and DFT/PCM calculation of the solute. The UV-visible spectrum of the solute is calculated with the help of the TDDFT method. The calculated spectrum is close to the experimental one. It is due to two species, a contact pair and an anion.

A new method for in vivo assessment of corneal transparency using spectral-domain OCT.

Corneal transparency is essential to provide a clear view into and out of the eye, yet clinical means to assess such transparency are extremely limited and usually involve a subjective grading of visible opacities by means of slit-lamp biomicroscopy. Here, we describe an automated algorithm allowing extraction of quantitative corneal transparency parameters with standard clinical spectral-domain optical coherence tomography (SD-OCT). Our algorithm employs a novel pre-processing procedure to standardize SD-OCT image analysis and to numerically correct common instrumental artifacts before extracting mean intensity stromal-depth (z) profiles over a 6-mm-wide corneal area. The z-profiles are analyzed using our previously developed objective method that derives quantitative transparency parameters directly related to the physics of light propagation in tissues. Tissular heterogeneity is quantified by the Birge ratio Br and the photon mean-free path (ls) is determined for homogeneous tissues (i.e., Br~1). SD-OCT images of 83 normal corneas (ages 22-50 years) from a standard SD-OCT device (RTVue-XR Avanti, Optovue Inc.) were processed to establish a normative dataset of transparency values. After confirming stromal homogeneity (Br <10), we measured a median ls of 570 µm (interdecile range: 270-2400 µm). By also considering corneal thicknesses, this may be translated into a median fraction of transmitted (coherent) light Tcoh(stroma) of 51% (interdecile range: 22-83%). Excluding images with central saturation artifact raised our median Tcoh(stroma) to 73% (interdecile range: 34-84%). These transparency values are slightly lower than those previously reported, which we attribute to the detection configuration of SD-OCT with a relatively small and selective acceptance angle. No statistically significant correlation between transparency and age or thickness was found. In conclusion, our algorithm provides robust and quantitative measurements of corneal transparency from standard SD-OCT images with sufficient quality (such as 'Line' and 'CrossLine' B-scan modes without central saturation artifact) and addresses the demand for such an objective means in the clinical setting.

A Novel Pump-Pump-Probe Resonance Raman Approach Featuring Light-Induced Charge Accumulation on a Model Photosystem.

Light-induced charge accumulation is at the heart of biomimetic systems aiming at solar fuel production in the realm of artificial photosynthesis. Understanding the mechanisms upon which these processes operate is a necessary condition to drive down the rational catalyst design road. We have built a nanosecond pump-pump-probe resonance Raman setup to witness the sequential charge accumulation process while probing vibrational features of different charge-separated states. By employing a reversible model system featuring methyl viologen (MV) as a dual electron acceptor, we have been able to watch the photosensitized production of its neutral form, MV0, resulting from two sequential electron transfer reactions. We have found that, upon double excitation, a fingerprint vibrational mode corresponding to the doubly reduced species appears at 992 cm-1 and peaks at 30 µs after the second excitation. This has been further confirmed by simulated resonance Raman spectra which fully support our experimental findings in this unprecedented buildup of charge seen by a resonance Raman probe.