The Gas Phase Protonation Sites of Six Naturally Occurring Nicotinoids.

The gas phase protonation sites of six naturally occurring nicotinoids, namely nicotine (NIC), nornicotine (NOR), anabasine (ANB), anatabine (ANT), cotinine (COT), and myosmine (MYO), consisting of a common Pyridine and differing non-Pyridine rings, have been determined for the first time at the physiological temperature from cryogenic ion trap infrared spectroscopy and electronic structure calculations. The protonation site on either of these two rings is related to the nicotinoid's biological activity. At room temperature, NIC is a mixture of Pyridine and Pyrrolidine (non-Pyridine) protomers, whereas NOR, ANB, ANT, and COT are pure Pyridine protomers and finally MYO is mostly a Pyroline (non-Pyridine) protomer. The nearly planar structure of MYO-H+, induced by the presence of a conjugated π system and confirmed from calculations and the UV absorption spectra, breaks from the trends observed for NIC, NOR, and ANB, since its structure is drastically different from the structures of the other nicotinoids.

Affinity of Nicotinoids to a Model Nicotinic Acetylcholine Receptor (nAChR) Binding Pocket in the Human Brain.

The binding affinity of nicotinoids to the binding residues of the α4ß2 variant of the nicotinic acetylcholine receptor (nAChR) was identified as a strong predictor of the nicotinoid's addictive character. Using ab initio calculations for model binding pockets of increasing size composed of 3, 6, and 14 amino acids (3AA, 6AA, and 14AA) that are derived from the crystal structure, the differences in binding affinity of 6 nicotinoids, namely, nicotine (NIC), nornicotine (NOR), anabasine (ANB), anatabine (ANT), myosmine (MYO), and cotinine (COT) were correlated to their previously reported doses required for increases in intracranial self-stimulation (ICSS) thresholds, a metric for their addictive function. By employing the many-body decomposition, the differences in the binding affinities of the various nicotinoids could be attributed mainly to the proton exchange energy between the pyridine and non-pyridine rings of the nicotinoids and the interactions between them and a handful of proximal amino acids, namely Trp156, Trpß57, Tyr100, and Tyr204. Interactions between the guest nicotinoid and the amino acids of the binding pocket were found to be mainly classical in nature, except for those between the nicotinoid and Trp156. The larger pockets were found to model binding structures more accurately and predicted the addictive character of all nicotinoids, while smaller models, which are more computationally feasible, would only predict the addictive character of nicotinoids that are similar to nicotine. The present study identifies the binding affinity of the guest nicotinoid to the host binding pocket as a strong descriptor of the nicotinoid's addiction potential, and as such it can be employed as a fast-screening technique for the potential addiction of nicotine analogs.

Descriptors of water aggregation.

We rely on a total of 23 (cluster size, 8 structural, and 14 connectivity) descriptors to investigate structural patterns and connectivity motifs associated with water cluster aggregation. In addition to the cluster size n (number of molecules), the 8 structural descriptors can be further categorized into (i) one-body (intramolecular): covalent OH bond length (rOH) and HOH bond angle (θHOH), (ii) two-body: OO distance (rOO), OHO angle (θOHO), and HOOX dihedral angle (ϕHOOX), where X lies on the bisector of the HOH angle, (iii) three-body: OOO angle (θOOO), and (iv) many-body: modified tetrahedral order parameter (q) to account for two-, three-, four-, five-coordinated molecules (qm, m = 2, 3, 4, 5) and radius of gyration (Rg). The 14 connectivity descriptors are all many-body in nature and consist of the AD, AAD, ADD, AADD, AAAD, AAADD adjacencies [number of hydrogen bonds accepted (A) and donated (D) by each water molecule], Wiener index, Average Shortest Path Length, hydrogen bond saturation (% HB), and number of non-short-circuited three-membered cycles, four-membered cycles, five-membered cycles, six-membered cycles, and seven-membered cycles. We mined a previously reported database of 4 948 959 water cluster minima for (H2O)n, n = 3-25 to analyze the evolution and correlation of these descriptors for the clusters within 5 kcal/mol of the putative minima. It was found that rOH and % HB correlated strongly with cluster size n, which was identified as the strongest predictor of energetic stability. Marked changes in the adjacencies and cycle count were observed, lending insight into changes in the hydrogen bond network upon aggregation. A Principal Component Analysis (PCA) was employed to identify descriptor dependencies and group clusters into specific structural patterns across different cluster sizes. The results of this study inform our understanding of how water clusters evolve in size and what appropriate descriptors of their structural and connectivity patterns are with respect to system size, stability, and similarity. The approach described in this study is general and can be easily extended to other hydrogen-bonded systems.

Acceleration of Graph Neural Network-Based Prediction Models in Chemistry via Co-Design Optimization on Intelligence Processing Units.

Atomic structure prediction and associated property calculations are the bedrock of chemical physics. Since high-fidelity ab initio modeling techniques for computing the structure and properties can be prohibitively expensive, this motivates the development of machine-learning (ML) models that make these predictions more efficiently. Training graph neural networks over large atomistic databases introduces unique computational challenges, such as the need to process millions of small graphs with variable size and support communication patterns that are distinct from learning over large graphs, such as social networks. We demonstrate a novel hardware-software codesign approach to scale up the training of atomistic graph neural networks (GNN) for structure and property prediction. First, to eliminate redundant computation and memory associated with alternative padding techniques and to improve throughput via minimizing communication, we formulate the effective coalescing of the batches of variable-size atomistic graphs as the bin packing problem and introduce a hardware-agnostic algorithm to pack these batches. In addition, we propose hardware-specific optimizations, including a planner and vectorization for the gather-scatter operations targeted for Graphcore's Intelligence Processing Unit (IPU), as well as model-specific optimizations such as merged communication collectives and optimized softplus. Putting these all together, we demonstrate the effectiveness of the proposed codesign approach by providing an implementation of a well-established atomistic GNN on the Graphcore IPUs. We evaluate the training performance on multiple atomistic graph databases with varying degrees of graph counts, sizes, and sparsity. We demonstrate that such a codesign approach can reduce the training time of atomistic GNNs and can improve their performance by up to 1.5× compared to the baseline implementation of the model on the IPUs. Additionally, we compare our IPU implementation with a Nvidia GPU-based implementation and show that our atomistic GNN implementation on the IPUs can run 1.8× faster on average compared to the execution time on the GPUs.

Switching of Protonation Sites in Hydrated Nicotine via a Grotthuss Mechanism.

The switching of the protonation sites in hydrated nicotine, probed by experimental infrared (IR) spectroscopy and theoretical ab initio calculations, is facilitated via a Grotthuss instead of a bimolecular proton transfer (vehicle) mechanism at the experimental temperature (T = 130 K) as unambiguously confirmed by experiments with deuterated water. In contrast, the bimolecular vehicle mechanism is preferred at higher temperatures (T = 300 K) as determined by theory. The Grotthuss mechanism for the concerted proton transfer results in the production of nicotine's bioactive and addictive pyrrolidine-protonated (Pyrro-H+) protomer with just 5 water molecules. Theoretical analysis suggests that the concerted proton transfer occurs via hydrogen-bonded bridges consisting of a 3 water molecule "core" that connects the pyridine protonated (Pyri-H+) with the pyrrolidine-protonated (Pyrro-H+) protomers. Additional water molecules attached as acceptors to the hydrogen-bonded "core" bridge result in lowering the reaction barrier of the concerted proton transfer down to less than 6 kcal/mol, which is consistent with the experimental observations.

The Back Door to the Surface Hydrated Electron.

We use a Mg+ metal to extend the size regime of aqueous clusters to extrapolate to the bulk limit of the vertical detachment energy (VDE) of the solvated electron to >3,200, a value between 1 to over 2 orders of magnitude larger than the one previously measured experimentally or computed theoretically. We relate the VDE to the energy difference between the Mg+(H2O)n and Mg2+(H2O)n systems and the metal's second ionization potential. The extrapolated bulk VDEs of the localized surface electron, which moves away from the metal as n increases, are 1.89 ± 0.01 eV for semiempirical (n ∼ 3,200; PM6-D3H4) and 1.73 ± 0.03 eV (n ∼ 150; HF) and 1.83 ± 0.02 eV (n ∼ 150; MP2) for ab initio, in excellent agreement with the 1.6-1.8 eV range of experimental results. The VDEs converge from above (larger values) to the bulk limit, in a manner that is qualitatively opposite from previous studies and experiments lacking a charged metal, a fact justifying the "back door" approach to the solvated electron.

Accurate Calculation of Many-Body Energies in Water Clusters Using a Classical Geometry-Dependent Induction Model.

We incorporate geometry-dependent distributed multipole and polarizability surfaces into an induction model that is used to describe the 3- and 4-body terms of the interaction between water molecules. The moment expansion is carried out up to the hexadecapole with the multipoles distributed on the atom sites. Dipole-dipole, dipole-quadrupole, and quadrupole-quadrupole distributed polarizabilities are used to represent the response of the multipoles to an electric field. We compare the model against two large databases consisting of 43,844 3-body terms and 3,603 4-body terms obtained from high level ab initio calculations previously used to fit the MB-pol and q-AQUA classical interaction potentials for water. The classical induction model with no adjustable parameters reproduces the ab initio 3-/4-body terms contained in these two databases with a root-mean-square error (RMSE) of 0.104/0.058 and a mean-absolute error (MAE) of 0.054/0.026 kcal/mol, respectively. These results are on par with the ones obtained by fitting the same data using over 14,000 (for the 3-body) and 200 (for the 4-body) parameters via Permutationally Invariant Polynomials (PIPs). This demonstrates the accuracy of this physically motivated model in describing the 3- and 4-body terms in the interactions between water molecules with no adjustable parameters. The triple-dipole-dispersion energy, included in the calculation of the 3-body energy, was found to be small but not quite negligible. The model represents a practical, efficient, and transferable approach for obtaining accurate nonadditive interactions for multicomponent systems without the need to perform tens of thousands of high level electronic structure calculations and fitting them with PIPs.

A Perspective on Sustainable Computational Chemistry Software Development and Integration.

The power of quantum chemistry to predict the ground and excited state properties of complex chemical systems has driven the development of computational quantum chemistry software, integrating advances in theory, applied mathematics, and computer science. The emergence of new computational paradigms associated with exascale technologies also poses significant challenges that require a flexible forward strategy to take full advantage of existing and forthcoming computational resources. In this context, the sustainability and interoperability of computational chemistry software development are among the most pressing issues. In this perspective, we discuss software infrastructure needs and investments with an eye to fully utilize exascale resources and provide unique computational tools for next-generation science problems and scientific discoveries.

Isotope Effects in the Zundel-Eigen Isomerization of H+(H2O)6.

The isomerization pathway between the energetically low-lying Zundel and Eigen isomers of the protonated water hexamer was investigated using high-level ab initio calculations including a treatment of zero-point corrections. On the basis of these calculations, the Zundel-Eigen isomerization was found to proceed through a stable intermediate isomer, which consists of a four-membered ring with two single acceptor water molecules. The inclusion of vibrational zero-point energy is shown to be important for accurately establishing the relative energies of the three relevant isomers involved in the Zundel-Eigen isomerization. Diffusion Monte Carlo calculations including anharmonic vibrational effects show that all three isomers of H+(H2O)6 and D+(D2O)6 have well-defined structures. The energetic ordering of the three isomers changes upon deuteration. The implications of these results for the vibrational spectra of H+(H2O)6 and D+(D2O)6 are also discussed.

Helical Organic and Inorganic Polymers.

Despite being a staple of synthetic plastics and biomolecules, helical polymers are scarcely studied with Gaussian-basis-set ab initio electron-correlated methods on an equal footing with molecules. This article introduces an ab initio second-order many-body Green's function [MBGF(2)] method with nondiagonal, frequency-dependent Dyson self-energy for infinite helical polymers using screw-axis-symmetry-adapted Gaussian-spherical-harmonics basis functions. Together with the Gaussian-basis-set density-functional theory for energies, analytical atomic forces, translational-period force, and helical-angle force, it can compute correlated energy, quasiparticle energy bands, structures, and vibrational frequencies of an infinite helical polymer, which smoothly converge at the corresponding oligomer results. These methods can handle incommensurable structures, which have an infinite translational period and are hard to characterize by any other method, just as efficiently as commensurable structures. We apply them to polyethylene (2/1 helix), polyacetylene (Peierls' system) and polytetrafluoroethylene (13/6 helix) to establish the quantitative accuracy of MBGF(2)/cc-pVDZ in simulating their (angle-resolved) ultraviolet photoelectron spectra and of B3LYP/cc-pVDZ or 6-31G** in reproducing their structures, infrared and Raman band positions, phonon dispersions, and (coherent and incoherent) inelastic neutron scattering spectra. We then predict the same properties for infinitely catenated chains of nitrogen or oxygen and discuss their possible metastable existence under ambient conditions. They include planar zigzag polyazene (N2)x (Peierls' system), 11/3-helical isotactic polyazane (NH)x, 9/4-helical isotactic polyfluoroazane (NF)x, and 7/2-helical polyoxane (O)x as potential high-energy-density materials.

Many-Body Effects in Aqueous Systems: Synergies Between Interaction Analysis Techniques and Force Field Development.

Interaction analysis techniques, including the many-body expansion (MBE), symmetry-adapted perturbation theory, and energy decomposition analysis, allow for an intuitive understanding of complex molecular interactions. We review these methods by first providing a historical context for the study of many-body interactions and discussing how nonadditivities emerge from Hamiltonians containing strictly pairwise-additive interactions. We then elaborate on the synergy between these interaction analysis techniques and the development of advanced force fields aimed at accurately reproducing the Born-Oppenheimer potential energy surface. In particular, we focus on ab initio-based force fields that aim to explicitly reproduce many-body terms and are fitted to high-level electronic structure results. These force fields generally incorporate many-body effects through (a) parameterization of distributed multipoles, (b) explicit fitting of the MBE, (c) inclusion of many-atom features in a neural network, and (d) coarse-graining of many-body terms into an effective two-body term. We also discuss the emerging use of the MBE to improve the accuracy and speed of ab initio molecular dynamics.

The primary gas phase hydration shell of hydroxide.

The number of water molecules in hydroxide's primary hydration shell has been long debated to be three from the interpretation of experimental data and four from theoretical studies. Here, we provide direct evidence for the presence of a fourth water molecule in hydroxide's primary hydration shell from a combined study based on high-resolution cryogenic experimental photoelectron spectroscopy and high-level quantum chemical computations. Well-defined spectra of OH-(H2O)n clusters (n = 2 to 5) yield accurate electron binding energies, which are, in turn, used as key signatures of the underlying molecular conformations. Although the smaller OH-(H2O)3 and OH-(H2O)4 clusters adopt close-lying conformations with similar electron binding energies that are hard to distinguish, the OH-(H2O)5 cluster clearly has a predominant conformation with a four-coordinated hydroxide binding motif, a finding that unambiguously determines the gas phase coordination number of hydroxide to be four.

An extensive assessment of the performance of pairwise and many-body interaction potentials in reproducing ab initio benchmark binding energies for water clusters n = 2-25.

We assess the performance of 7 pairwise additive (TIP3P, TIP4P, TIP4P-ice, TIP5P, OPC, SPC, SPC/E) and 8 families of many-body potentials (q-AQUA, HIPPO, AMOEBA, EFP, TTM, WHBB, MB-pol, MB-UCB) in reproducing high-level ab initio benchmark values, CCSD(T) or MP2 at the complete basis set (CBS) limit for the binding energy and the many-body expansion (MBE) of water clusters n = 2-11, 16-17, 20, 25. By including a large range of cluster sizes having dissimilar hydrogen bonding networks, we obtain an understanding of how these potentials perform for different hydrogen bonding arrangements that are mostly outside of their parameterization range. While it is appropriate to compare the results of ab initio based many-body potentials directly to the electronic binding energies (De's), the pairwise additive ones are compared to the enthalpies at T = 298 K, ΔH(298 K), as the latter class of force fields are parametrized to reproduce enthalpies (implicitly accounting for zero-point energy corrections) rather than binding energies. We find that all pairwise additive potentials considered overestimate the reference ΔH values for the n = 2-25 clusters by >13%. For the water dimer (n = 2) in particular, the errors are in the range 83-119% for the pairwise additive potentials studied since these are based on an effective rather than the true 2-body interaction specifically designed as a means of partially accounting for the missing many-body terms. This stronger 2-body interaction is achieved by an enhanced monomer dipole moment that mimics its increase from the gas phase monomer to the condensed phase value. Indeed, for cluster sizes n ≥ 4 the percent deviations become slightly smaller (albeit all exceeding 13%). In contrast, we find that the many-body potentials perform more accurately in reproducing the electronic binding energies (De's) throughout the entire cluster range (n = 2-25), all reproducing the ab initio benchmark binding energies within ±7% of the respective CBS values. We further assess the ability of a subset of the many-body potentials (MB-UCB, q-AQUA, MB-pol, and TTM2.1-F) to also reproduce the magnitude of the ab initio many-body energy terms for water cluster sizes n = 7, 10, 16 and 17. The potentials show an overall good agreement with the available benchmark values. However, we identify characteristic differences upon comparing the many-body terms at both the ab initio-optimized geometries and the respective potential-optimized geometries to the reference ab initio values. Additionally, by applying this analysis to a wide range of cluster sizes, trends in the MBE of the potentials with increasing cluster size can be identified. Finally, in an attempt to draw a parallel between the pairwise additive and many-body potentials, we report the analysis of the individual molecular dipole moments for water clusters with 1 to ∼4 solvation shells with the TTM2.1-F potential. We find that the internally solvated water molecules have in general a larger molecular dipole moment ranging from 2.6-3.0 D. This justifies the use of an enhanced, with respect to the gas-phase value, molecular dipole moment for the pairwise additive potentials, which is intended to fold in the many body terms into an effective (enhanced) pairwise interaction through the choice of the charges. These results have important implications for the development of future generations of efficient, transferable, and highly accurate classical interaction potentials in both the pairwise additive and many-body categories.

A Formulation of the Many-Body Expansion (MBE) for Periodic Systems: Application to Several Ice Phases.

We introduce a new formulation of the many-body expansion (MBE) for periodic systems and apply it to 7 ice polymorphs (Ih, II, VIII, IX, XIII, XIV, and XV). This new formulation is built via a hierarchical procedure that connects gas-phase clusters that mimic unit cells over finite supercells to infinite solids. For periodic systems, the method is validated by showing that the lattice energies computed up to the 4-body in the MBE reproduce the lattice energies obtained using periodic boundary conditions with an Ewald summation for the 7 ice polymorphs. This development makes it possible to quantify, for the first time, the many-body contributions to the lattice energy of various ice polymorphs, which vary significantly among the 7 ice phases, amounting to between 7 and 24% of the total lattice energies. This development opens the door for obtaining insights into solid-state properties, while leveraging the computational benefits of the MBE.

A critical comparison of CH⋯π versus π⋯π interactions in the benzene dimer: obtaining benchmarks at the CCSD(T) level and assessing the accuracy of lower scaling methods.

We have established CCSD(T)/CBS (Complete Basis Set) limits for 3 stationary points on the benzene dimer potential energy surface, corresponding to the π⋯π (parallel displaced or PD(C2h), minimum) and CH⋯π (T-shaped or T(C2v), transition state) and tilted T-shaped (or TT(Cs), minimum) bonding scenarios considering both the structure and binding energy. The CCSD(T)/CBS binding energies are -2.65 ± 0.02 (PD), -2.74 ± 0.03 (T), and -2.83 ± 0.01 kcal mol-1 (TT). To this end, the CH⋯π is ∼0.2 kcal mol-1 stronger than the π⋯π interaction, whereas the tilting of the CH donating benzene molecule with respect to the other benzene is worth 0.1 kcal mol-1. As previously discussed in the literature, the MP2 level of theory does not provide a close match for either the energy or structure, yet the SCS-MP2 yields structures in excellent agreement with respect to the CCSD(T) result. It is found that the SCS-MI-MP2 also gives optimized structures very close to SCS-MP2 (within ∼0.01 Å of the benchmark). Despite the closer match in structure, the spin-biased MP2 methods (SCS-, SCS-MI-, and SOS-MP2) incorrectly predict the relative stabilities of the isomers. That said, none of the spin biased MP2 methods offers a good compromise between energy and structure for the systems examined. Finally, the CCSD(T)/CBS benchmarks were used to assess the performance of 13 DFT functionals selected from different rungs of Jacob's ladder. Several functionals such as TPSS-D3, B3LYP-D3, B97-D, B97-D3, and B2PLYP-D3 provided a good description of the binding energies for both CH⋯π and π⋯π interactions, yielding values within 6% of the CCSD(T)/CBS benchmark values. Unlike the MP2 methods, these functionals correctly predict the relative stability of the PD(C2h) and T(C2v) dimers. Further, we find that there is no systematic improvement as Jacob's ladder is ascended (increased complexity of functional). The best functionals that result in a good compromise between structure and energy accuracy are B97-D3 and B2PLYP-D3 for both the CH⋯π and π⋯π interaction. Despite the impressive performance of these functionals, a challenge that remains is ensuring the transferability of these density functionals in accurately describing the interaction between dimers of larger aromatic molecules, the latter requiring high-level benchmarks for these systems.