From NWChem to NWChemEx: Evolving with the Computational Chemistry Landscape.

Since the advent of the first computers, chemists have been at the forefront of using computers to understand and solve complex chemical problems. As the hardware and software have evolved, so have the theoretical and computational chemistry methods and algorithms. Parallel computers clearly changed the common computing paradigm in the late 1970s and 80s, and the field has again seen a paradigm shift with the advent of graphical processing units. This review explores the challenges and some of the solutions in transforming software from the terascale to the petascale and now to the upcoming exascale computers. While discussing the field in general, NWChem and its redesign, NWChemEx, will be highlighted as one of the early codesign projects to take advantage of massively parallel computers and emerging software standards to enable large scientific challenges to be tackled.

A parallel, distributed memory implementation of the adaptive sampling configuration interaction method.

The many-body simulation of quantum systems is an active field of research that involves several different methods targeting various computing platforms. Many methods commonly employed, particularly coupled cluster methods, have been adapted to leverage the latest advances in modern high-performance computing. Selected configuration interaction (sCI) methods have seen extensive usage and development in recent years. However, the development of sCI methods targeting massively parallel resources has been explored only in a few research works. Here, we present a parallel, distributed memory implementation of the adaptive sampling configuration interaction approach (ASCI) for sCI. In particular, we will address the key concerns pertaining to the parallelization of the determinant search and selection, Hamiltonian formation, and the variational eigenvalue calculation for the ASCI method. Load balancing in the search step is achieved through the application of memory-efficient determinant constraints originally developed for the ASCI-PT2 method. The presented benchmarks demonstrate near optimal speedup for ASCI calculations of Cr2 (24e, 30o) with 106, 107, and 3 × 108 variational determinants on up to 16 384 CPUs. To the best of the authors' knowledge, this is the largest variational ASCI calculation to date.

Distributed memory, GPU accelerated Fock construction for hybrid, Gaussian basis density functional theory.

With the growing reliance of modern supercomputers on accelerator-based architecture such a graphics processing units (GPUs), the development and optimization of electronic structure methods to exploit these massively parallel resources has become a recent priority. While significant strides have been made in the development GPU accelerated, distributed memory algorithms for many modern electronic structure methods, the primary focus of GPU development for Gaussian basis atomic orbital methods has been for shared memory systems with only a handful of examples pursing massive parallelism. In the present work, we present a set of distributed memory algorithms for the evaluation of the Coulomb and exact exchange matrices for hybrid Kohn-Sham DFT with Gaussian basis sets via direct density-fitted (DF-J-Engine) and seminumerical (sn-K) methods, respectively. The absolute performance and strong scalability of the developed methods are demonstrated on systems ranging from a few hundred to over one thousand atoms using up to 128 NVIDIA A100 GPUs on the Perlmutter supercomputer.

PluginPlay: Enabling exascale scientific software one module at a time.

For many computational chemistry packages, being able to efficiently and effectively scale across an exascale cluster is a heroic feat. Collective experience from the Department of Energy's Exascale Computing Project suggests that achieving exascale performance requires far more planning, design, and optimization than scaling to petascale. In many cases, entire rewrites of software are necessary to address fundamental algorithmic bottlenecks. This in turn requires a tremendous amount of resources and development time, resources that cannot reasonably be afforded by every computational science project. It thus becomes imperative that computational science transition to a more sustainable paradigm. Key to such a paradigm is modular software. While the importance of modular software is widely recognized, what is perhaps not so widely appreciated is the effort still required to leverage modular software in a sustainable manner. The present manuscript introduces PluginPlay, https://github.com/NWChemEx-Project/PluginPlay, an inversion-of-control framework designed to facilitate developing, maintaining, and sustaining modular scientific software packages. This manuscript focuses on the design aspects of PluginPlay and how they specifically influence the performance of the resulting package. Although, PluginPlay serves as the framework for the NWChemEx package, PluginPlay is not tied to NWChemEx or even computational chemistry. We thus anticipate PluginPlay to prove to be a generally useful tool for a number of computational science packages looking to transition to the exascale.

Computing Free Energies with Fluctuation Relations on Quantum Computers.

As a central thermodynamic property, free energy enables the calculation of virtually any equilibrium property of a physical system, allowing for the construction of phase diagrams and predictions about transport, chemical reactions, and biological processes. Thus, methods for efficiently computing free energies, which in general is a difficult problem, are of great interest to broad areas of physics and the natural sciences. The majority of techniques for computing free energies target classical systems, leaving the computation of free energies in quantum systems less explored. Recently developed fluctuation relations enable the computation of free energy differences in quantum systems from an ensemble of dynamic simulations. While performing such simulations is exponentially hard on classical computers, quantum computers can efficiently simulate the dynamics of quantum systems. Here, we present an algorithm utilizing a fluctuation relation known as the Jarzynski equality to approximate free energy differences of quantum systems on a quantum computer. We discuss under which conditions our approximation becomes exact, and under which conditions it serves as a strict upper bound. Furthermore, we successfully demonstrate a proof of concept of our algorithm using the transverse field Ising model on a real quantum processor. As quantum hardware continues to improve, we anticipate that our algorithm will enable computation of free energy differences for a wide range of quantum systems useful across the natural sciences.

Quantum Algorithm for High Energy Physics Simulations.

Simulating quantum field theories is a flagship application of quantum computing. However, calculating experimentally relevant high energy scattering amplitudes entirely on a quantum computer is prohibitively difficult. It is well known that such high energy scattering processes can be factored into pieces that can be computed using well established perturbative techniques, and pieces which currently have to be simulated using classical Markov chain algorithms. These classical Markov chain simulation approaches work well to capture many of the salient features, but cannot capture all quantum effects. To exploit quantum resources in the most efficient way, we introduce a new paradigm for quantum algorithms in field theories. This approach uses quantum computers only for those parts of the problem which are not computable using existing techniques. In particular, we develop a polynomial time quantum final state shower that accurately models the effects of intermediate spin states similar to those present in high energy electroweak showers with a global evolution variable. The algorithm is explicitly demonstrated for a simplified quantum field theory on a quantum computer.

Mitigating Depolarizing Noise on Quantum Computers with Noise-Estimation Circuits.

A significant problem for current quantum computers is noise. While there are many distinct noise channels, the depolarizing noise model often appropriately describes average noise for large circuits involving many qubits and gates. We present a method to mitigate the depolarizing noise by first estimating its rate with a noise-estimation circuit and then correcting the output of the target circuit using the estimated rate. The method is experimentally validated on a simulation of the Heisenberg model. We find that our approach in combination with readout-error correction, randomized compiling, and zero-noise extrapolation produces close to exact results even for circuits containing hundreds of CNOT gates. We also show analytically that zero-noise extrapolation is improved when it is applied to the output of our method.

Reactions of U+ with H2, D2, and HD Studied by Guided Ion Beam Tandem Mass Spectrometry and Theory.

The kinetic energy-dependent reactions of the atomic actinide uranium cation (U+) with H2, D2, and HD were examined by guided ion beam tandem mass spectrometry. An average 0 K bond dissociation energy of D0(U+ - H) = 2.48 ± 0.06 eV is obtained by analysis of the endothermic product ion cross sections. Quantum chemistry calculations were performed for comparison with experimental thermochemistry, including high-level CASSCF-CASPT2-RASSI calculations of the spin-orbit corrections. CCSD(T) and the CASSCF levels show excellent agreement with experiment, whereas B3LYP and PBE0 slightly overestimate and the M06 approach badly underestimates the bond energy for UH+. Theory was also used to investigate the electronic structures of the reaction intermediates and potential energy surfaces. The experimental product branching ratio for the reaction of U+ with HD indicates that these reactions occur primarily via a direct reaction mechanism, despite the presence of a deep-well for UH2+ formation according to theory. The reactivity and hydride bond energy for U+ are compared with those for transition metal, lanthanide, and actinide cations, and periodic trends are discussed. These comparisons suggest that the 5f electrons on uranium are largely core and uninvolved in the reactive chemistry.

Prediction of atomization energy using graph kernel and active learning.

Data-driven prediction of molecular properties presents unique challenges to the design of machine learning methods concerning data structure/dimensionality, symmetry adaption, and confidence management. In this paper, we present a kernel-based pipeline that can learn and predict the atomization energy of molecules with high accuracy. The framework employs Gaussian process regression to perform predictions based on the similarity between molecules, which is computed using the marginalized graph kernel. To apply the marginalized graph kernel, a spatial adjacency rule is first employed to convert molecules into graphs whose vertices and edges are labeled by elements and interatomic distances, respectively. We then derive formulas for the efficient evaluation of the kernel. Specific functional components for the marginalized graph kernel are proposed, while the effects of the associated hyperparameters on accuracy and predictive confidence are examined. We show that the graph kernel is particularly suitable for predicting extensive properties because its convolutional structure coincides with that of the covariance formula between sums of random variables. Using an active learning procedure, we demonstrate that the proposed method can achieve a mean absolute error of 0.62 ± 0.01 kcal/mol using as few as 2000 training samples on the QM7 dataset.

Stochastic resolution of identity second-order Matsubara Green's function theory.

We develop a stochastic resolution of identity representation to the second-order Matsubara Green's function (sRI-GF2) theory. Using a stochastic resolution of the Coulomb integrals, the second order Born self-energy in GF2 is decoupled and reduced to matrix products/contractions, which reduces the computational cost from O(N5) to O(N3) (with N being the number of atomic orbitals). The current approach can be viewed as an extension to our previous work on stochastic resolution of identity second order Møller-Plesset perturbation theory [T. Y. Takeshita et al., J. Chem. Theory Comput. 13, 4605 (2017)] and offers an alternative to previous stochastic GF2 formulations [D. Neuhauser et al., J. Chem. Theory Comput. 13, 5396 (2017)]. We show that sRI-GF2 recovers the deterministic GF2 results for small systems, is computationally faster than deterministic GF2 for N > 80, and is a practical approach to describe weak correlations in systems with 103 electrons and more.

Revealing Disparate Chemistries of Protactinium and Uranium. Synthesis of the Molecular Uranium Tetroxide Anion, UO4.

The synthesis, reactivity, structures, and bonding in gas-phase binary and complex oxide anion molecules of protactinium and uranium have been studied by experiment and theory. The oxalate ions, AnVO2(C2O4)-, where An = Pa or U, are essentially actinyl ions, AnVO2+, coordinated by an oxalate dianion. Both react with water to yield the pentavalent hydroxides, AnVO(OH)2(C2O4)-. The chemistry of Pa and U becomes divergent for reactions that result in oxidation: whereas PaVI is inaccessible, UVI is very stable. The UVO2(C2O4)- complex exhibits a remarkable spontaneous exothermic replacement of the oxalate ligand by O2 to yield UO4- and two CO2 molecules. The structure of the uranium tetroxide anion is computed to correspond to distorted uranyl, UVIO22+, coordinated in the equatorial plane by two equivalent O atoms each having formal charges of -1.5 and U-O bond orders intermediate between single and double. The unreactive nature of PaVO2(C2O4)- toward O2 is a manifestation of the resistance toward oxidation of PaV, and clearly reveals the disparate chemistries of Pa and U. The uranium tetroxide anion, UO4-, reacts with water to yield UO5H2-. Infrared spectra obtained for UO5H2- confirm the computed lowest-energy structure, UO3(OH)2-.

Cleaving Off Uranyl Oxygens through Chelation: A Mechanistic Study in the Gas Phase.

Recent efforts to activate the strong uranium-oxygen bonds in the dioxo uranyl cation have been limited to single oxo-group activation through either uranyl reduction and functionalization in solution, or by collision induced dissociation (CID) in the gas-phase, using mass spectrometry (MS). Here, we report and investigate the surprising double activation of uranyl by an organic ligand, 3,4,3-LI(CAM), leading to the formation of a formal U6+ chelate in the gas-phase. The cleavage of both uranyl oxo bonds was experimentally evidenced by CID, using deuterium and 18O isotopic substitutions, and by infrared multiple photon dissociation (IRMPD) spectroscopy. Density functional theory (DFT) computations predict that the overall reaction requires only 132 kJ/mol, with the first oxygen activation entailing about 107 kJ/mol. Combined with analysis of similar, but unreactive ligands, these results shed light on the chelation-driven mechanism of uranyl oxo bond cleavage, demonstrating its dependence on the presence of ligand hydroxyl protons available for direct interactions with the uranyl oxygens.

Heptavalent Actinide Tetroxides NpO4- and PuO4-: Oxidation of Pu(V) to Pu(VII) by Adding an Electron to PuO4.

The highest known actinide oxidation states are Np(VII) and Pu(VII), both of which have been identified in solution and solid compounds. Recently a molecular Np(VII) complex, NpO3(NO3)2-, was prepared and characterized in the gas phase. In accord with the lower stability of heptavalent Pu, no Pu(VII) molecular species has been identified. Reported here are the gas-phase syntheses and characterizations of NpO4- and PuO4-. Reactivity studies and density functional theory computations indicate the heptavalent metal oxidation state in both. This is the first instance of Pu(VII) in the absence of stabilizing effects due to condensed phase solvation or crystal fields. The results indicate that addition of an electron to neutral PuO4, which has a computed electron affinity of 2.56 eV, counterintuitively results in oxidation of Pu(V) to Pu(VII), concomitant with superoxide reduction.

Pushing configuration-interaction to the limit: Towards massively parallel MCSCF calculations.

A new large-scale parallel multiconfigurational self-consistent field (MCSCF) implementation in the open-source NWChem computational chemistry code is presented. The generalized active space approach is used to partition large configuration interaction (CI) vectors and generate a sufficient number of batches that can be distributed to the available cores. Massively parallel CI calculations with large active spaces can be performed. The new parallel MCSCF implementation is tested for the chromium trimer and for an active space of 20 electrons in 20 orbitals, which can now routinely be performed. Unprecedented CI calculations with an active space of 22 electrons in 22 orbitals for the pentacene systems were performed and a single CI iteration calculation with an active space of 24 electrons in 24 orbitals for the chromium tetramer was possible. The chromium tetramer corresponds to a CI expansion of one trillion Slater determinants (914 058 513 424) and is the largest conventional CI calculation attempted up to date.

Strengthening of the Coordination Shell by Counter Ions in Aqueous Th4+ Solutions.

The presence of counterions in solutions containing highly charged metal cations can trigger processes such as ion-pair formation, hydrogen bond breakages and subsequent re-formation, and ligand exchanges. In this work, it is shown how halide (Cl-, Br-) and perchlorate (ClO4-) anions affect the strength of the primary solvent coordination shells around Th4+ using explicit-solvent and finite-temperature ab initio molecular dynamics modeling methods. The 9-fold solvent geometry was found to be the most stable hydration structure in each aqueous solution. Relative to the dilute aqueous solution, the presence of the counterions did not significantly alter the geometry of the primary hydration shell. However, the free energy analyses indicated that the 10-fold hydrated states were thermodynamically accessible in dilute and bromide aqueous solutions within 1 kcal/mol. Analysis of the results showed that the hydrogen bond lifetimes were longer and solvent exchange energy barriers were larger in solutions with counterions in comparison with the solution with no counterions. This implies that the presence of the counterions induces a strengthening of the Th4+ hydration shell.

Gas phase uranyl activation: formation of a uranium nitrosyl complex from uranyl azide.

Activation of the oxo bond of uranyl, UO2(2+), was achieved by collision induced dissociation (CID) of UO2(N3)Cl2(-) in a quadrupole ion trap mass spectrometer. The gas phase complex UO2(N3)Cl2(-) was produced by electrospray ionization of solutions of UO2Cl2 and NaN3. CID of UO2(N3)Cl2(-) resulted in the loss of N2 to form UO(NO)Cl2(-), in which the "inert" uranyl oxo bond has been activated. Formation of UO2Cl2(-) via N3 loss was also observed. Density functional theory computations predict that the UO(NO)Cl2(-) complex has nonplanar Cs symmetry and a singlet ground state. Analysis of the bonding of the UO(NO)Cl2(-) complex shows that the side-on bonded NO moiety can be considered as NO(3-), suggesting a formal oxidation state of U(VI). Activation of the uranyl oxo bond in UO2(N3)Cl2(-) to form UO(NO)Cl2(-) and N2 was computed to be endothermic by 169 kJ/mol, which is energetically more favorable than formation of NUOCl2(-) and UO2Cl2(-). The observation of UO2Cl2(-) during CID is most likely due to the absence of an energy barrier for neutral ligand loss.

Activation of CH4 by Th(+) as studied by guided ion beam mass spectrometry and quantum chemistry.

The reaction of atomic thorium cations with CH4 (CD4) and the collision-induced dissociation (CID) of ThCH4(+) with Xe are studied using guided ion beam tandem mass spectrometry. In the methane reactions at low energies, ThCH2(+) (ThCD2(+)) is the only product; however, the energy dependence of the cross-section is inconsistent with a barrierless exothermic reaction as previously assumed on the basis of ion cyclotron resonance mass spectrometry results. The dominant product at higher energies is ThH(+) (ThD(+)), with ThCH3(+) (ThCD3(+)) having a similar threshold energy. The latter product subsequently decomposes at still higher energies to ThCH(+) (ThCD(+)). CID of ThCH4(+) yields atomic Th(+) as the exclusive product. The cross-sections of all product ions are modeled to provide 0 K bond dissociation energies (in eV) of D0(Th(+)-H) ≥ 2.25 ± 0.18, D0(Th(+)-CH) = 6.19 ± 0.16, D0(Th(+)-CH2) ≥ 4.54 ± 0.09, D0(Th(+)-CH3) = 2.60 ± 0.30, and D0(Th(+)-CH4) = 0.47 ± 0.05. Quantum chemical calculations at several levels of theory are used to explore the potential energy surfaces for activation of methane by Th(+), and the effects of spin-orbit coupling are carefully considered. When spin-orbit coupling is explicitly considered, a barrier for C-H bond activation that is consistent with the threshold measured for ThCH2(+) formation (0.17 ± 0.02 eV) is found at all levels of theory, whereas this barrier is observed only at the BHLYP and CCSD(T) levels otherwise. The observation that the CID of the ThCH4(+) complex produces Th(+) as the only product with a threshold of 0.47 eV indicates that this species has a Th(+)(CH4) structure, which is also consistent with a barrier for C-H bond activation. This barrier is thought to exist as a result of the mixed ((4)F,(2)D) electronic character of the Th(+) J = (3)/2 ground level combined with extensive spin-orbit effects.

Indirect dynamics in a highly exoergic substitution reaction.

The highly exoergic nucleophilic substitution reaction F(-) + CH3I shows reaction dynamics strikingly different from that of substitution reactions of larger halogen anions. Over a wide range of collision energies, a large fraction of indirect scattering via a long-lived hydrogen-bonded complex is found both in crossed-beam imaging experiments and in direct chemical dynamics simulations. Our measured differential scattering cross sections show large-angle scattering and low product velocities for all collision energies, resulting from efficient transfer of the collision energy to internal energy of the CH3F reaction product. Both findings are in strong contrast to the previously studied substitution reaction of Cl(-) + CH3I [Science 2008, 319, 183-186] at all but the lowest collision energies, a discrepancy that was not captured in a subsequent study at only a low collision energy [J. Phys. Chem. Lett. 2010, 1, 2747-2752]. Our direct chemical dynamics simulations at the DFT/B97-1 level of theory show that the reaction is dominated by three atomic-level mechanisms, an indirect reaction proceeding via an F(-)-HCH2I hydrogen-bonded complex, a direct rebound, and a direct stripping reaction. The indirect mechanism is found to contribute about one-half of the overall substitution reaction rate at both low and high collision energies. This large fraction of indirect scattering at high collision energy is particularly surprising, because the barrier for the F(-)-HCH2I complex to form products is only 0.10 eV. Overall, experiment and simulation agree very favorably in both the scattering angle and the product internal energy distributions.

Cation-cation interactions in [(UO2)2(OH)n](4-n) complexes.

The structures and bonding of gas-phase [(UO2)2(OH)n](4-n) (n = 2-6) complexes have been studied using density functional theory (DFT), MP2, and CCSD(T) methods with particular emphasis on ground state structures featuring cation-cation interactions (CCIs) between the uranyl groups. An interesting trend is observed in the stabilities of members of this series of complexes. The structures of [(UO2)2(OH)2](2+), [(UO2)2(OH)4], and [(UO2)2(OH)6](2-) featuring CCIs are found at higher energies (by 3-27 kcal/mol) in comparison to their conventional µ2-dihydroxo structures. In contrast, the CCI structures of [(UO2)2(OH)3](+) and [(UO2)2(OH)5](-) are respectively degenerate with and lower in energy than the structures with the µ2-dihydroxo format. The origin of this trend lies in the symmetry-based need to balance the coordination numbers and effective atomic charges of each uranium center. The calculated IR vibrational frequencies provide signature probes that can be used in differentiating the low-energy structures and in experimentally confirming the existence of the structures featuring CCIs.

Coordination and hydrolysis of plutonium ions in aqueous solution using Car-Parrinello molecular dynamics free energy simulations.

Car-Parrinello molecular dynamics (CPMD) simulations have been used to examine the hydration structures, coordination energetics, and the first hydrolysis constants of Pu(3+), Pu(4+), PuO2(+), and PuO2(2+) ions in aqueous solution at 300 K. The coordination numbers and structural properties of the first shell of these ions are in good agreement with available experimental estimates. The hexavalent PuO2(2+) species is coordinated to five aquo ligands while the pentavalent PuO2(+) complex is coordinated to four aquo ligands. The Pu(3+) and Pu(4+) ions are both coordinated to eight water molecules. The first hydrolysis constants obtained for Pu(3+) and PuO2(2+) are 6.65 and 5.70, respectively, all within 0.3 pH unit of the experimental values (6.90 and 5.50, respectively). The hydrolysis constant of Pu(4+), 0.17, disagrees with the value of -0.60 in the most recent update of the Nuclear Energy Agency Thermochemical Database (NEA-TDB) but supports recent experimental findings. The hydrolysis constant of PuO2(+), 9.51, supports the experimental results of Bennett et al. [Radiochim. Acta 1992, 56, 15]. A correlation between the pKa of the first hydrolysis reaction and the effective charge of the plutonium center was found.