Reactive path integral quantum simulations of molecules solvated in superfluid helium

Abstract

A hybrid ab initio path integral molecular dynamics/bosonic path integral Monte Carlo simulation method has been developed, implemented and tested, which allows for the reactive simulations of molecules, clusters or complexes solvated by superfluid ${}^4$He. The simulation takes into account “on-the-fly” the electronic structure and thus the chemical reactivity of the solutes, in conjunction with the Bose–Einstein statistics, and thus the superfluid character of this peculiar solvent. This enables investigations into cryochemical reactions taking place in helium nanodroplets, such as those used in helium nanodroplet isolation (HENDI) spectroscopy.

Introduction

Traditional cryochemistry developed two core experimental methodologies that gave most fruitful access to processes occurring at (ultra-) low temperatures. In the condensed phase, matrix isolation techniques are used in order to trap transient species, mostly radicals and highly reactive intermediates, in non-reactive noble gas environments that are kept at low temperatures in order to serve as host matrices. In the gas phase, seeded supersonic beams and various ion traps are successfully employed to perform spectroscopy on cold molecules, complexes and clusters that would not be stable otherwise. Much more recent, dating back to 1992, is the discovery of using helium nanodroplets as spectroscopically transparent carriers for atoms and molecules  [1], which stimulated the rapid development of the helium nanodroplet isolation (HENDI) spectroscopy  [2], [3], [4], [5], [6], [7], [8], [9], [10]. Helium, which undergoes a phase transition from normal to superfluid state at 2.17 K in the bulk phase, provides a gentle, ultracold environment capable of adjusting itself readily to the guest molecule, as opposed to other low-temperature media. On the other hand the pick-up process, which can be precisely controlled by the dopant vapor pressure, and the negligible friction inside the droplet enables aggregation of rare species, which overcomes limitations of other low-temperature techniques. Both advantages, together with the high resolution and other virtues of laser spectroscopy, led to numerous experimental investigations into the structure of molecules, clusters and complexes embedded in such an unusual environment at a temperature of less than one Kelvin.

At about the same time as these investigations, quantum simulation techniques of superfluid clusters were developed with the aim of understanding the properties of helium and other bosonic liquids at low temperatures. The pioneering finite-temperature quantum simulations of small ${}^4$He clusters below the $\lambda$ transition temperature were able to predict the superfluid behavior of droplets with only 64 atoms  [11], which ten years later was impressively confirmed experimentally  [12]. This led to the emergence of numerous applications, mostly based on quantum Monte Carlo techniques, which dealt with pure helium clusters at the early stage (see Refs.  [13], [14], [15], [16], [17] for review articles and references given therein for original work). Further early theoretical efforts concentrated on the solvation of molecular impurities in the helium solvent  [18], [19], [20], [21], which explained experimental results, such as the layering of the solvation shells around the solute, reduced rotational constants of molecules solvated in helium with respect to the gas phase and so on. However, these approaches investigated either single atoms and ions that have no internal degrees of freedom  [22], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37], [38], [39], or rigid molecules kept fixed at their equilibrium structure [40], [41], [42], [43], [44], [45], [46], [47], [48], [49], [50], [51], [52], [53], [54], [55], [56], [57], [58], [59], [60], [61], [62], [63], [64], [65], [66], [67], [68], [69], [70], [71], [72], [73], [74], [75], [76], [77], [78]. Clusters consisting of such rigid molecules interacting with pre-calculated “force field” potentials were also investigated along the same lines  [79], [80]. While allowing for the study of the solvation of various molecular systems by superfluid helium, these methods prevent examination of chemical transformations, such as chemical reactions, or even of large-amplitude motion. Such phenomena, although observed experimentally  [81], [82], certainly require deeper understanding, since traditionally invoked explanations seem inappropriate in this context, as discussed below.

The effort of the present work aims at pushing the simulation technology for superfluid helium clusters from non-reactive molecules to reactive ones, thus enabling the study of cryochemical reactions therein. A good example of such a reactive process in superfluid helium nanodroplets is the recently observed dissociation of a HCl molecule upon microsolvation: upon adding one water molecule after the other acid dissociation is observed upon adding the fourth water molecule  [81], [82]. But how are such processes possible at all at temperatures where the thermal energy $k_{B}T$, enabling traditional Arrhenius activation, is negligible compared to typical activation barriers? At first sight, the old concept of quantum tunneling might appear promising to explain cryochemical reactions  [83], but this seems inappropriate to explain these processes occurring in very asymmetric settings. Based on extensive ab initio simulations the novel mechanism of “aggregation-induced chemical reactions” was proposed recently to explain dissociation of HCl in microsolvation situations at the sub-Kelvin temperatures of helium nanodroplets  [81], [84]. According to this idea, successive assembly of water molecules accompanied by the stepwise formation of the aggregate molecule by molecule is mostly energetically downhill, moreover being facilitated by the kinetic energy released upon creation of hydrogen bonds as a result of aggregation, which is in turn used to overcome small energy barriers. The sequence of aggregation events eventually leads to the creation of an “activated species”, that readily undergoes a chemical reaction, after having aggregated a certain critical number of molecules. As this model was developed for the gas phase, an interesting question is if and how aggregation-induced chemical reactions are influenced by the superfluid helium environment. The coupled path integral MD/MC method being presented here aims at attacking this and other related phenomena.

Another class of questions that motivated the development of the present method refer to fluxional (or floppy) molecules or complexes which are characterized by multi-dimensional large-amplitude motion. A prominent example, inspired by pioneering messenger vibrational spectroscopy experiments [85], [86], is protonated methane, CH${}_5^{+}$, which undergoes vivid pseudo-rotational transitions between equivalent $C_s$ structures even at very low temperatures where nuclear quantum effects are key to the process [87], [88]. This fast, essentially barrierless interconversion, known as “full hydrogen scrambling”, can be restricted by attaching solvent molecules, e.g. molecular hydrogen  [89], which makes the problem more tractable to spectroscopic measurements. Although much has been learned about this system, including the assignment of its broadband vibrational spectrum  [90] as well as those of all its H/D isotopologues  [91], the detailed effect of microsolvation on its fluxionality remains controversial. Similarly, in the case of hydronium and Zundel complexes (i.e. the protonated water monomer H₃O⁺ and dimer H₅O₂⁺), the vibrational spectra are known to depend on the number and chemical nature of the particular messenger species, thus leaving their fingerprints in the spectra [92], [93]. An interesting question is if the bulk-like solvation inside helium affects large-amplitude motion of floppy molecules in a way analogous to messenger species in the gas phase? The superfluid helium environment would not only be a perfect matrix to investigate the impact of attaching messenger molecules to “fragile species” in the limit of very weak coupling, but also to probe if the molecule–helium interactions might be able to affect the fluxionality of bare protonated methane or the properties of bare protonated water clusters in the first place.

The method to be outlined in detail in the following sections rests on using ab initio path integral molecular dynamics (PIMD) [94], [95], [96], [97], [98] to treat the embedded molecule or cluster, including explicitly its electronic structure “on-the-fly”, which is the chemically reactive solute. Here, we use the cost-efficient Generalized Gradient Approximation (GGA) to Kohn–Sham density functional theory in order to generate the ground-state intra-solute interactions along the PIMD trajectory within the Born–Oppenheimer approximation. Of course more elaborate methods such as wavefunction-based approaches could be used—if practical. This reactive solute species is coupled to the explicit but chemically inert helium solvation environment via a parameterized interaction potential  [99] (much like in well-established hybrid quantum-mechanics/molecular-mechanics QM/MM approaches) fitted to coupled cluster electronic structure reference data. The superfluidity of the helium solvent is established by using bosonic path integral Monte Carlo (PIMC) for indistinguishable ⁴He atoms that interact pairwise via a parameterized interatomic potential following Ceperley and Pollock  [14], [100], [101]. The coupling of the quantum solute to the quantum solvent, in the sense of a hybrid PIMD/PIMC sampling scheme, is established based on ideas of LaBerge and Tully  [102]. In summary, the first “hybrid” aspect of our technique concerns the treatment of the solute/solvent interactions, whereas the second “hybrid” feature is related to the statistical sampling of the coupled solute/solvent density matrix. The first one, “QM/MM”, is an approximation since the full electronic structure of the helium solvation environment is replaced by a force field, and thus its interaction with the reactive solute, which precludes any changes of the solvent’s electronic structure being consistent with considering the solvent to be chemically inert. In our particular case of a reactive solute in an inert solvent, it should be noted that no covalent bonds pass through the QM–MM interface but only noncovalent solute–solvent interactions. Next, given that He–He interactions are well known, e.g. in terms of the effective He–He pair potential due to Aziz and coworkers, and given that the He–solute interactions have been computed at the coupled cluster level and subsequently carefully parameterized for the specific reactive species treated herein  [99], it can be anticipated that the QM/MM Ansatz should be an excellent approximation to study these solute–solvent systems. The second feature, “hybrid MD/MC”, is rigorous since, the algorithm to couple PIMD propagation of the reactive solute to PIMC sampling of the surrounding chemically inert solvent assures –by construction –the correct canonical sampling of the equilibrium density matrix of the full system. It is noted in passing that, of course, no dynamical information is provided by the molecular dynamics propagation of the ab initio path integral describing the solute subsystem, where “time” is used as a convenient way to parameterize statistical sampling instead of counting configurations or moves as in PIMC.

The method is implemented in the CP2K simulation package  [103], drawing on the Quickstep module for density functional based electronic structure calculations  [104], and thus enables full finite-temperature quantum simulations of reactive solute molecules or clusters in a superfluid helium environment, which can be a few ${}^4$He atoms, a finite droplet or periodic bulk. The ability of the algorithm to enable the investigation of chemical reactions taking place in superfluid solvent is demonstrated, in the sense of a proof of concept, by simulating the elementary chemical reaction of the dissociation of an acid molecule inside a small water cluster, which in turn is embedded in superfluid helium, thus mimicking the situation met in typical HENDI experiments.