Direct path integral estimators for isotope fractionation ratios

Proton quantal delocalization and H/D translocations in (MeOH)nH+ (n = 2, 3)

In this study, we present results from path integral molecular dynamics simulations that describe the characteristics of the quantum spatial delocalizations of protons participating in OH bonds in (MeOH)2H+ and in (MeOH)3H+. The characterization was carried out by examining the overall structures of the corresponding isomorphic polymers. To introduce full flexibility in the force treatment, we have adopted a neural network fitting procedure based on second-order Møller-Plesset perturbation theory predictions. For the dimer case, we found that the spatial extent of the shared connective proton can be portrayed in terms of a prolate-like structure with typical dimensions of ∼0.1 Å. On the other hand, the dangling polymers lie confined within a thin spherical layer, spread over length scales of the order of ∼0.25 Å. In contrast, connective protons in (MeOH)3H+ exhibit larger delocalizations along the O-H bond and more localized ones along perpendicular directions, compared to their dangling counterparts. We also examined the characteristics of the relative propensities of H and D isotopes to be localized in dangling and connective positions. Physical interpretations of the different thermodynamic trends are provided in terms of the local geometrical characteristics and of the strengths of the corresponding intermolecular connectivities.

Path Integral Simulations of Condensed-Phase Vibrational Spectroscopy

Recent theoretical and algorithmic developments have improved the accuracy with which path integral dynamics methods can include nuclear quantum effects in simulations of condensed-phase vibrational spectra. Such methods are now understood to be approximations to the delocalized classical Matsubara dynamics of smooth Feynman paths, which dominate the dynamics of systems such as liquid water at room temperature. Focusing mainly on simulations of liquid water and hexagonal ice, we explain how the recently developed quasicentroid molecular dynamics (QCMD), fast-QCMD, and temperature-elevated path integral coarse-graining simulations (Te PIGS) methods generate classical dynamics on potentials of mean force obtained by averaging over quantum thermal fluctuations. These new methods give very close agreement with one another, and the Te PIGS method has recently yielded excellent agreement with experimentally measured vibrational spectra for liquid water, ice, and the liquid-air interface. We also discuss the limitations of such methods.

Routine Molecular Dynamics Simulations Including Nuclear Quantum Effects: From Force Fields to Machine Learning Potentials

We report the implementation of a multi-CPU and multi-GPU massively parallel platform dedicated to the explicit inclusion of nuclear quantum effects (NQEs) in the Tinker-HP molecular dynamics (MD) package. The platform, denoted Quantum-HP, exploits two simulation strategies: the Ring-Polymer Molecular Dynamics (RPMD) that provides exact structural properties at the cost of a MD simulation in an extended space of multiple replicas and the adaptive Quantum Thermal Bath (adQTB) that imposes the quantum distribution of energy on a classical system via a generalized Langevin thermostat and provides computationally affordable and accurate (though approximate) NQEs. We discuss some implementation details, efficient numerical schemes, and parallelization strategies and quickly review the GPU acceleration of our code. Our implementation allows an efficient inclusion of NQEs in MD simulations for very large systems, as demonstrated by scaling tests on water boxes with more than 200,000 atoms (simulated using the AMOEBA polarizable force field). We test the compatibility of the approach with Tinker-HP's recently introduced Deep-HP machine learning potentials module by computing water properties using the DeePMD potential with adQTB thermostatting. Finally, we show that the platform is also compatible with the alchemical free energy estimation capabilities of Tinker-HP and fast enough to perform simulations. Therefore, we study how NQEs affect the hydration free energy of small molecules solvated with the recently developed Q-AMOEBA water force field. Overall, the Quantum-HP platform allows users to perform routine quantum MD simulations of large condensed-phase systems and will help to shed new light on the quantum nature of important interactions in biological matter.

Quantum-mechanical exploration of the phase diagram of water

The set of known stable phases of water may not be complete, and some of the phase boundaries between them are fuzzy. Starting from liquid water and a comprehensive set of 50 ice structures, we compute the phase diagram at three hybrid density-functional-theory levels of approximation, accounting for thermal and nuclear fluctuations as well as proton disorder. Such calculations are only made tractable because we combine machine-learning methods and advanced free-energy techniques. The computed phase diagram is in qualitative agreement with experiment, particularly at pressures ≲ 8000 bar, and the discrepancy in chemical potential is comparable with the subtle uncertainties introduced by proton disorder and the spread between the three hybrid functionals. None of the hypothetical ice phases considered is thermodynamically stable in our calculations, suggesting the completeness of the experimental water phase diagram in the region considered. Our work demonstrates the feasibility of predicting the phase diagram of a polymorphic system from first principles and provides a thermodynamic way of testing the limits of quantum-mechanical calculations.

Adsorption separation of heavier isotope gases in subnanometer carbon pores

Isotopes of heavier gases including carbon (¹³C/¹⁴C), nitrogen (¹³N), and oxygen (¹⁸O) are highly important because they can be substituted for naturally occurring atoms without significantly perturbing the biochemical properties of the radiolabelled parent molecules. These labelled molecules are employed in clinical radiopharmaceuticals, in studies of brain disease and as imaging probes for advanced medical imaging techniques such as positron-emission tomography (PET). Established distillation-based isotope gas separation methods have a separation factor (S) below 1.05 and incur very high operating costs due to high energy consumption and long processing times, highlighting the need for new separation technologies. Here, we show a rapid and highly selective adsorption-based separation of ¹⁸O₂ from ¹⁶O₂ with S above 60 using nanoporous adsorbents operating near the boiling point of methane (112 K), which is accessible through cryogenic liquefied-natural-gas technology. A collective-nuclear-quantum effect difference between the ordered ¹⁸O₂ and ¹⁶O₂ molecular assemblies confined in subnanometer pores can explain the observed equilibrium separation and is applicable to other isotopic gases.

Identification of synthesisable crystalline phases of water – a prototype for the challenges of computational materials design

We discuss the identification of experimentally realisable crystalline phases of water to outline and contextualise some of the diverse building blocks of a computational materials design process.

Accurate Description of Nuclear Quantum Effects with High-Order Perturbed Path Integrals (HOPPI)

Imaginary time path-integral (PI) simulations that account for nuclear quantum effects (NQE) beyond the harmonic approximation are increasingly employed together with modern electronic-structure calculations. Existing PI methods are applicable to molecules, liquids, and solids; however, the computational cost of such simulations increases dramatically with decreasing temperature. To address this challenge, here, we propose to combine high-order PI factorization with perturbation theory (PT). Already for conventional second-order PI simulations, the PT ansatz increases the accuracy 2-fold compared to fourth-order schemes with the same settings. In turn, applying PT to high-order path integrals (HOPI) further improves the efficiency of simulations for molecular and condensed matter systems especially at low temperatures. We present results for bulk liquid water, the aspirin molecule, and the CH₅⁺ molecule. Perturbed HOPI simulations remain both efficient and accurate down to 20 K and provide a convenient method to estimate the convergence of quantum-mechanical observables.

A combined on-the-fly/interpolation procedure for evaluating energy values needed in molecular simulations

We propose an algorithm for molecular dynamics or Monte Carlo simulations that uses an interpolation procedure to estimate potential energy values from energies and gradients evaluated previously at points of a simplicial mesh. We chose an interpolation procedure that is exact for harmonic systems and considered two possible mesh types: Delaunay triangulation and an alternative anisotropic triangulation designed to improve performance in anharmonic systems. The mesh is generated and updated on the fly during the simulation. The procedure is tested on two-dimensional quartic oscillators and on the path integral Monte Carlo evaluation of the HCN/DCN equilibrium isotope effect.

Accelerating equilibrium isotope effect calculations. II. Stochastic implementation of direct estimators

Path integral calculations of equilibrium isotope effects and isotopic fractionation are expensive due to the presence of path integral discretization errors, statistical errors, and thermodynamic integration errors. Whereas the discretization errors can be reduced by high-order factorization of the path integral and statistical errors by using centroid virial estimators, two recent papers proposed alternative ways to completely remove the thermodynamic integration errors: Cheng and Ceriotti [J. Chem. Phys. 141, 244112 (2015)] employed a variant of free-energy perturbation called "direct estimators," while Karandashev and Vaníček [J. Chem. Phys. 143, 194104 (2017)] combined the thermodynamic integration with a stochastic change of mass and piecewise-linear umbrella biasing potential. Here, we combine the former approach with the stochastic change in mass in order to decrease its statistical errors when applied to larger isotope effects and perform a thorough comparison of different methods by computing isotope effects first on a harmonic model and then on methane and methanium, where we evaluate all isotope effects of the form CH_4-xD_x/CH₄ and CH_5-xD_x ⁺/CH₅ ⁺, respectively. We discuss the reasons for a surprising behavior of the original method of direct estimators, which performed well for a much larger range of isotope effects than what had been expected previously, as well as some implications of our work for the more general problem of free energy difference calculations.

Ab initio thermodynamics of liquid and solid water

Thermodynamic properties of liquid water as well as hexagonal (Ih) and cubic (Ic) ice are predicted based on density functional theory at the hybrid-functional level, rigorously taking into account quantum nuclear motion, anharmonic fluctuations, and proton disorder. This is made possible by combining advanced free-energy methods and state-of-the-art machine-learning techniques. The ab initio description leads to structural properties in excellent agreement with experiments and reliable estimates of the melting points of light and heavy water. We observe that nuclear-quantum effects contribute a crucial [Formula: see text] to the stability of ice Ih, making it more stable than ice Ic. Our computational approach is general and transferable, providing a comprehensive framework for quantitative predictions of ab initio thermodynamic properties using machine-learning potentials as an intermediate step.

Dual QM and MM Approach for Computing Equilibrium Isotope Fractionation Factor of Organic Species in Solution

A dual QM and MM approach for computing equilibrium isotope effects has been described. In the first partition, the potential energy surface is represented by a combined quantum mechanical and molecular mechanical (QM/MM) method, in which a solute molecule is treated quantum mechanically, and the remaining solvent molecules are approximated classically by molecular mechanics. In the second QM/MM partition, differential nuclear quantum effects responsible for the isotope effect are determined by a statistical mechanical double-averaging formalism, in which the nuclear centroid distribution is sampled classically by Newtonian molecular dynamics and the quantum mechanical spread of quantized particles about the centroid positions is treated using the path integral (PI) method. These partitions allow the potential energy surface to be properly represented such that the solute part is free of nuclear quantum effects for nuclear quantum mechanical simulations, and the double-averaging approach has the advantage of sampling efficiency for solvent configuration and for path integral convergence. Importantly, computational precision is achieved through free energy perturbation (FEP) theory to alchemically mutate one isotope into another. The PI-FEP approach is applied to model systems for the 18O enrichment found in cellulose of trees to determine the isotope enrichment factor of carbonyl compounds in water. The present method may be useful as a general tool for studying isotope fractionation in biological and geochemical systems.

Can Path Integral Molecular Dynamics Make a Good Approximation for Vapor Pressure Isotope Effects Prediction for Organic Solvents? A Comparison to ONIOM QM/MM and QM Cluster Calculation

Isotopic fractionation of volatile organic compounds (VOCs), which are under strict measures of control because of their potential harm to the environment and humans, has an important ecological aspect, as the isotopic composition of compounds may depend on the conditions in which such compounds are distributed in Nature. Therefore, detailed knowledge on isotopic fractionation, not only experimental but also based on theoretical models, is crucial to follow conditions and pathways within which these contaminants are spread throughout the ecosystems. In this work, we present carbon and, for the first time, bromine vapor pressure isotope effect (VPIE) on the evaporation process from pure-phase systems-dibromomethane and bromobenzene, the representatives of aliphatic and aromatic brominated VOCs. We combine isotope effects measurements with their theoretical prediction using three computational techniques, namely path integral molecular dynamics, QM cluster, and hybrid ONIOM models. While evaporation of both compounds resulted in normal bromine VPIEs, the difference in the direction of carbon isotopic fractionation is observed for the aliphatic and aromatic compounds, where VPIEs are inverse and normal, respectively. Even though theoretical models tested here turned out to be insufficient for quantitative agreement with the experimental values, cluster electronic structure calculations, as well as two-layer ONIOM computations, provided better reproduction of experimental trends.

Nuclear Quantum Effects in Water at the Triple Point: Using Theory as a Link Between Experiments

One of the most prominent consequences of the quantum nature of light atomic nuclei is that their kinetic energy does not follow a Maxwell-Boltzmann distribution. Deep inelastic neutron scattering (DINS) experiments can measure this effect. Thus, the nuclear quantum kinetic energy can be probed directly in both ordered and disordered samples. However, the relation between the quantum kinetic energy and the atomic environment is a very indirect one, and cross-validation with theoretical modeling is therefore urgently needed. Here, we use state of the art path integral molecular dynamics techniques to compute the kinetic energy of hydrogen and oxygen nuclei in liquid, solid, and gas-phase water close to the triple point, comparing three different interatomic potentials and validating our results against equilibrium isotope fractionation measurements. We will then show how accurate simulations can draw a link between extremely precise fractionation experiments and DINS, therefore establishing a reliable benchmark for future measurements and providing key insights to increase further the accuracy of interatomic potentials for water.