Superfluid helium has not only fascinated scientists for centuries but is also the ideal matrix for the investigation of chemical systems under ultra-cold conditions in helium nanodroplet isolation experiments. Together with related experimental techniques such as helium tagging photodissociation spectroscopy, these methods have provided unique insights into many interesting systems. Complemented by theoretical work, they were additionally able to greatly expand our general understanding of manifestations of superfluid behavior in finite sized clusters and their response to molecular impurities. However, most theoretical studies up to now have not included the reactivity and flexibility of molecular systems embedded in helium. In this perspective, the theoretical foundation of simulating fluxional molecules and reactive complexes in superfluid helium is presented in detail. Special emphasis is put on recent developments for the converged description of both the molecular interactions and the quantum nature of the nuclei at ultra-low temperatures. As a first step, our hybrid path integral molecular dynamics/bosonic path integral Monte Carlo method is reviewed. Subsequently, methods for efficient path integral sampling tailored for this hybrid coupling scheme are discussed while also introducing new developments to enhance the accurate incorporation of the solute⋯solvent coupling. Finally, highly accurate descriptions of the interactions in solute⋯helium systems using machine learning techniques are addressed. Our current automated and adaptive fitting procedures to parameterize high-dimensional neural network potentials for both the full-dimensional potential energy surface of solutes and the solute⋯solvent interaction potentials are concisely presented. They are demonstrated to faithfully represent many-body potential functions able to describe chemically complex and reactive solutes in helium environments seamlessly from one He atom up to bulk helium at the accuracy level of coupled cluster electronic structure calculations. Together, these advances allow for converged quantum simulations of fluxional and reactive solutes in superfluid helium under cryogenic conditions.
Experimental techniques such as vibrational photodissociation spectroscopy using messenger tagging1–4 or helium nanodroplet isolation (HENDI) experiments5–10 have become increasingly popular in recent years for the investigation of complex molecules, clusters, and chemical reactions. Both classes of experimental techniques use helium either as a very gentle tagging agent or to fully embed the chemical system of interest in a superfluid helium environment. They have led to great insights into the properties of a wealth of different systems such as neutral11–17 and protonated water clusters,18,19 HCl⋯water aggregates,20,21 reactive intermediates22 such as or metal clusters,23–28 and amino acid peptides or even small proteins,29–32 to name but a few. Moreover, such experiments turned out to be fundamentally important for understanding the bosonic nature of few to many interacting 4He atoms (as well as of para-H2 molecules) and the so-called molecular superfluidity33,34 via the superfluid response of bosonic helium around molecular impurities.5,35–43
While spanning the full range from small sized helium aggregates to bulk like superfluid helium nanodroplets, the experimental techniques operate under unique temperature conditions on the order of a few kelvin in the case of photodissociation tagging spectroscopy or even at sub-kelvin temperatures35 as established in HENDI experiments. This grants access to the system of interest close to its rovibrational ground state essentially in the absence of thermal fluctuations. For a detailed understanding of these systems under such ultra-cold conditions, quantum simulation techniques that take nuclear quantum effects (NQEs) into account are key to provide insights into the underlying physical and chemical processes from a fully atomistic point of view. They can therefore assist, complement, and guide experimental studies to shed light on many fundamental questions. Pioneering such quantum simulations were, for example, able to predict44 that finite clusters consisting of only 64 4He atoms already feature “manifestations of superfluid behavior” below about 2 K, which has been confirmed experimentally ten years later.37 Many of these simulation techniques are based on the Feynman–Kac formulation of quantum statistical mechanics via discretized imaginary time path integrals.45–48 This formalism is able to provide, in principle, the exact quantum structural and thermodynamic properties of chemically complex and even reactive many-body systems at finite temperatures, notably also including Bose–Einstein statistics49 and thus superfluid response.48
Indeed, atomistic quantum simulation approaches mainly based on finite-temperature path integral Monte Carlo (PIMC) techniques, but also using various quantum Monte Carlo approaches to access the quantum ground state, have greatly expanded our knowledge of single atoms and ions50–61 as well as of chemically inert non-reactive molecules.38,41,42,62–86 These vary from systems interacting with a few up to many helium atoms with a particular focus on the quantum (micro-)solvation shell structure around these solute species and on the molecular superfluidity of that bosonic quantum solvent. Most of the pioneering studies have used force fields for the description of the molecular impurities in superfluid helium, where vibrational degrees of freedom were usually neglected, while rigid-body translations and rigid-body rotations were increasingly taken into account. There also exist dedicated hybrid schemes where exclusively inter-molecular degrees of freedom within molecular complexes or a single coordinate describing large-amplitude motion within a molecule are explicitly considered in the sampling together with helium,65,73,86 while other degrees of freedom such as the intra-molecular vibrations are kept frozen. Yet, a large body of literature exists where the molecular solutes are even fully clamped in space in the sense of space-fixed classical point particles, contributing only a fixed external potential that acts on the helium atoms.
Although these PIMC approaches provide deep insights into the quantum solvation of chemically complex molecular systems at ultra-low temperatures, they allow one to study neither chemical reactions nor intra-molecular vibrations in the presence of bosonic helium, the latter being of overriding importance for flexible, floppy, or fluxional molecules. In contrast, quantum simulation techniques based on path integral molecular dynamics (PIMD)87–89 and, in particular, ab initio path integral (AIPI) simulations90–94 are ideally suited to study NQEs on fully flexible molecules, clusters, or molecular complexes including even covalent chemical reactions via the latter approach.95–99 Such molecular dynamics based sampling, unfortunately, is not readily applicable to superfluid solutions since bosonic exchange cannot be efficiently included in PIMD for systems consisting of many bosonic particles despite promising recent progress to be addressed later.
The recently developed hybrid path integral molecular dynamics/bosonic path integral Monte Carlo (HPIMD/MC) method100 combines above-mentioned advantages by efficiently sampling the bosonic solvent using path integral Monte Carlo while treating the fully flexible and reactive solute using ab initio path integral molecular dynamics. This has, for the first time, opened the door to investigate quantum solutes including their full reactivity inside superfluid helium environments. The pioneering simulations100,101 of reactive HCl⋯water clusters with this hybrid technique relied on “on-the-fly” evaluations of the interactions within the solute species based on concurrent electronic structure calculations in the spirit of AIPI simulations.90–94 However, this approach imposes rather high computational costs and is usually limited to a density functional theory based description of the electronic structure. In the end, this severely limits the convergence of these simulations in many respects, namely, in terms of (i) the statistical error of the computed properties, (ii) the systematic error due to the finite path integral discretization, (iii) the accuracy of the underlying electronic structure method for describing the solute itself, and (iv) the solute⋯helium interactions. In certain cases, computationally very efficient but specialized force fields can also be applied in the case of molecular fluxionality, such as developed for the highly flexible protonated methane molecule,102 for which, indeed, intriguing couplings between the large-amplitude motion and the bosonic exchange in the helium surrounding were disclosed only recently.103 However, this approach has only poor capability for generalization to other interesting and larger chemical systems and, in addition, is limited in terms of accuracy by the underlying functional form, in particular when it comes to describing complex chemical reactions. An additional complication of essentially all simulation approaches for solute⋯helium complexes is that carefully developed, system-specific interaction potentials for the coupling of the solute and the He solvent are essential.104–106 This has the disadvantage that for each system of interest, new interaction potentials with helium have to be computed and analytically fitted, which turns out to be a daunting task if molecular flexibility and chemical reactivity need to be taken into account.105,106
In the following, we review in detail how essentially fully converged quantum simulations of fluxional and reactive solutes in superfluid helium can be performed due to most recent methodological advances while also presenting some unpublished developments that contribute to solving the above-mentioned limitations (i)–(iv). For this purpose, we first introduce our HPIMD/MC method100 and also allude to various other simulation approaches. Subsequently, we address the issue of converging path integral sampling at finite but ultra-low temperatures within our hybrid approach by showing how the challenges that arise from ultra-cold conditions107,108 in conjunction with the bosonic quantum nature of helium can be met. Thereafter, we outline a general strategy that provides computationally efficient and largely automated access to accurate descriptions of both inter- and intra-molecular interactions. Our particular focus is on recent developments based on machine learning techniques that rely on artificial neural network representations of both solute⋯helium109 and intra-solute110 full-dimensional potential energy surfaces (PES) at the level of essentially converged coupled cluster electronic structure calculations. We finally conclude with an outlook on the simulation of reactive impurities in bosonic quantum solvents and spell out remaining challenges.
II. HYBRID PATH INTEGRAL MOLECULAR DYNAMICS/BOSONIC PATH INTEGRAL MONTE CARLO
In the following, we first briefly discuss different approaches to simulate molecular impurities in superfluid helium before reviewing the main aspects of our general hybrid strategy for reactive solutes in quantum solvents. Most finite temperature path integral approaches for the simulation of molecules in superfluid helium rely on path integral Monte Carlo (PIMC) sampling for which various algorithms have been specifically designed to describe the rotation of rigid molecular impurities in 4He, e.g., established in Refs. 34, 71, 77, 81, and 111–113. However, as already mentioned in the Introduction, PIMC and QMC approaches, in general, are not ideally suited for the efficient simulation of chemically complex molecular systems if molecular vibrations or even chemical reactions shall be included since a lot of the standard sampling moves would usually be rejected. Thus, algorithms introducing specific MC moves are needed to ensure efficient sampling in some of these cases, which have been devised repeatedly such as rigid-body DMC to sample inter-molecular interactions, while the interacting molecules themselves are treated as frozen entities.114–118 Rigorous ground state QMC methods have also been tailored for the specific needs of simulating molecules in helium environments.119–122 Although granting access to the correct bosonic ground state, these approaches prevent the analysis of finite temperature effects on properties such as the superfluid density.
Therefore, in order to sample the quantum partition function of a molecular system solvated by a bosonic solvent, it is advantageous to split the system into two coupled subsystems, which can be sampled individually using the most efficient method. This is the main idea behind the HPIMD/MC method introduced in Ref. 100, where the interested reader can find all details, while we focus here on a short summary of the fundamental concepts. Each subsystem is described with a sampling algorithm that is best suited for its individual challenges. Molecular impurities usually feature high flexibility and possibly reactivity and, thus, are best described by PIMD propagation. This offers great sampling advantages for chemically complex molecular systems that are governed by highly anisotropic and strong intra-molecular interactions, since the forces required in molecular dynamics propagation guide best the sampling of the partition function, while MC moves relying only on energy differences would be rejected most of the time (except in the limit of using tiny displacements where, however, sampling of configuration space becomes prohibitively slow or upon introducing system-specific tailored MC moves). Quantum solvents that are subject to relatively weak intermolecular interactions, such as 4He but also para-H2, are efficiently sampled by PIMC, which additionally enables one to include their bosonic nature. These two methods, which will be described in detail in Sec. III, are combined using the rigorous so-called MDMC hybrid sampling scheme introduced originally by LaBerge and Tully for a distinctly different purpose, namely, for simulating classical fluids.123
In our HPIMD/MC approach, the total system is evolved in two different phases, as schematically depicted in Fig. 1(a). During the PIMD phase of the sampling, the solute is propagated according to its equations of motion including the forces due to the intra-molecular interactions, Umol, and all those due to the helium environment, Uint,
which is kept fixed during this phase. The helium environment acts essentially as an external potential in this phase of the sampling algorithm as in the classical MDMC method.123 Note that the forces need only to be evaluated for the last (Monte Carlo) accepted configuration of the solvent according to the hybrid scheme introduced in Ref. 123. Afterward, the solvent is sampled as a PIMC walk through helium configuration and “permutation space” (see Sec. III for details) while now keeping all solute positions fixed according to MDMC sampling.123 The interaction potential is evaluated for each helium trial move and contributes to the change in potential energy,
New helium configurations are then accepted according to the usual Metropolis criterion124 applied to the Boltzmann-weighted change in the thermal path integral action. As discussed in Sec. III, it is important to perform sufficiently many MC steps during this phase to efficiently move through configuration and especially “permutation space” of the solute. These two phases alternate and together allow for the correct sampling of molecular systems solvated by superfluid helium. In Sec. III, details of the path integral sampling for the two subsystems are presented, which allows for a converged and ergodic description of reactive systems in superfluid helium.
III. PATH INTEGRAL SAMPLING
A. Path integral Monte Carlo sampling of bosonic helium
For the simulation of bosonic helium, one has to account for quantum-mechanical particle exchange, which is at the very heart of many quantum phenomena such as the Lambda transition from normal to superfluid helium that occurs at 2.17 K in bulk helium at ambient pressure. For 4He, Bose–Einstein statistics has to be obeyed, which means that the wavefunction and the corresponding density matrix need to be symmetric under exchange of any two atoms. This implies that in addition to the configuration space, the “permutation space” has to be sampled if quantum exchange statistics shall be enforced. In the path integral isomorphism of a classical system of P beads connected by harmonic springs, this means that originally closed paths (also called “ring polymers” in some literature) of different originally distinguishable atoms can now be connected to each other,49 thus lifting the distinguishable nature of identical particles according to Maxwell–Boltzmann statistics and thereby establishing Bose–Einstein statistics. Such global (topological) changes are achieved in a non-continuous manner so that sampling via molecular dynamics is significantly hampered, whereas suitable global PIMC moves can be introduced to rigorously sample symmetric exchange permutations, thus generating superfluid helium at finite temperatures. For further background and details on the theory and implementation of bosonic path integrals, we refer to the classic review by Ceperley.48
Two highly efficient sampling algorithms based on PIMC have been developed for the task of simultaneously sampling permutation and configuration space in an ergodic manner. We note only in passing that in the realm of PIMD sampling, there have been several attempts to also generalize PIMD to bosonic systems.125–128 This also includes largely129 unpublished efforts by two of us, which in the end were not useful for simulating superfluid helium in view of the strongly repulsive He⋯He interactions that result in a radial distribution function that is essentially zero up to a distance of about 2 Å, thereby strongly suppressing the particle overlap and thus exchange permutations. Very recently,130 a PIMD method has been proposed that avoids explicit enumeration of permutations of bosons subject to a very weakly repulsive model interaction potential resulting in a pair correlation function that is still sizable even for zero distance. In stark contrast to PIMD, successful PIMC methods to sample symmetric permutations are known since long. The first bosonic PIMC algorithm, called the bisection algorithm,48,131 relies on resampling sections of a ring polymer by a Lévy construction that guarantees acceptance of the move for a free particle. For weakly interacting particles, this provides a very effective move with rather high acceptance ratios. Permutation space sampling is incorporated by resampling more than one section of different ring polymers at the same time, where the interconnection is swapped. Note that this includes parts of imaginary time paths of different atoms, thus making these moves rather unfavorable in the general case. In order to increase the probability of accepting such a change of the permutation state, the free particle transition probability of pair permutations between neighboring atoms can be used as a bias in the selection of appropriate permutation changes.48 This algorithm has been implemented quite some time ago by HF in the CP2k program package132,133 and underlies our original ab initio HPIMD/MC method.100,101
Alternatively, the so-called worm algorithm has been introduced for efficient sampling of bosonic systems.134,135 This technique opens the cyclic ring polymers, which transforms the sampling of the diagonal density matrix into off-diagonal elements [cf. open move in Fig. 1(d)]. Afterward, the opened ring polymer, called worm, can connect with that of other atoms via similar MC moves as the ones used for the bisection algorithm, thus effectively changing the permutation state of the system [cf. swap move in Fig. 1(g)]. However, in contrast to the bisection algorithm, only individual imaginary time sections need to be sampled and accepted, thus increasing the overall acceptance probability of permutation changes. After closing [cf. Fig. 1(d)], which moves the sampling back to the diagonal density matrix elements, the system ends up in a different permutation state. Originally introduced in the grand canonical ensemble,134,135 the algorithm can be reduced to a canonical version by a suitable combination of the original moves and omission of some others.136 Our canonical version of the worm algorithm (as introduced in the Appendix and Fig. 1) has recently been implemented by us in the CP2k program package.132,133 As usual in MC sampling, specialized moves can be designed for the sake of improving efficiency, such as the centroid move [cf. Fig. 1(b)] and other moves to resample different parts of the ring polymer [cf. Figs. 1(c), 1(e), and 1(f)] in order to enhance permutation sampling also in the limit of small finite systems where a molecule is microsolvated by only a few helium atoms, which, however, fairly strongly interact with the solute.103,137 The different moves used in our CP2k implementation of the canonical worm algorithm are described in detail in the Appendix.
Due to the simple additive two-body nature of the He⋯He potential, the pair density approximation has been developed to efficiently converge the path integral discretization.48,138 We note that the pair potential is not only a good approximation in the bulk phase of helium, but also in the helium microsolvation shell around even fairly strongly interacting molecules as demonstrated in Sec. III C of Ref. 106, where the three-body contribution to the interaction energy has been shown to not exceed one percent as estimated by explicit coupled cluster electronic structure calculations. The pair density approximation allows one to efficiently perform PIMC simulations of bosonic helium even at temperatures on the order of 1 K or less. In this approach, the many-body density matrix is expanded in one-body and two-body contributions to the many-body action, where the nontrivial two-body contribution to the total action is approximated by the sum over the two-body contributions of the exact pair action for all pairs of helium atoms. These two-body contributions can be efficiently calculated using the matrix squaring method combined with the partial wave expansion.48 In practice, they are precomputed and tabulated on a regular grid, and further fitted to a few one-dimensional spline tables for quick evaluation during the simulation as implemented long ago by HF in the CP2k program package.132,133 Overall, the number of replicas P required to converge the path integral describing the system can be significantly reduced by this approach (see Ref. 48 for more details).
Having access to the bosonic density matrix at sufficiently low, yet finite temperatures allows one to compute the superfluid response.48 The superfluid fraction fs, which is a central property for the study of superfluidity, can be computed based on different estimators in PIMC simulations. The estimator commonly applied for finite sized systems comes from linear response theory and relies on the mean-squared area of the paths (i.e., the area spanned by the ring polymer) sampled for a container at rest divided by the classical moment of inertia,44
where is the projected area,
along some direction and Ic is the corresponding classical moment of inertia, which can be expressed as
The results obtained from this and related estimators should be interpreted with care, especially in inhomogeneous systems containing molecular impurities that can impact both the projected area and the moment of inertia. In addition, the estimator will still result in non-zero superfluid fractions even for very small clusters—notably even for a single helium atom that per definition cannot exchange. More recent work on molecules in small bosonic para-H2 clusters therefore proposed to renormalize the superfluid density (yielding the so-called exchange superfluid fraction) based on additional simulations that exclude bosonic exchange139 essentially to account approximately for artificial contributions to the projected area that is unrelated to bosonic exchange paths.34
Other work has been devoted to defining locally resolved superfluid response,63,67,140–142 where pragmatic proxies for superfluidity based on simple cutoffs for sufficiently long exchange cycle lengths were also proposed.63,142,143 Clearly, superfluidity is a genuinely macroscopic concept, and any local perspective or analysis of finite bosonic systems introduces necessarily specific viewpoints beyond what is a physical observable in the strict sense. It is for these reasons wise to speak about “manifestations of superfluid behavior”44 or to use the notion of “molecular superfluidity”33,34 instead of simply calling significant bosonic exchange phenomena “superfluidity” in such a context.
In periodic boundary conditions, another estimator has been derived that relies on the so-called winding number of paths wrapping through the periodic box,48,144 whereas the condensate fraction of superfluid helium has been computed for the first time in Ref. 131. This was actually a breakthrough in the mid-1980s toward quantifying superfluidity in quantum simulations. In such a case, the classical moment of inertia becomes and the projected area can be written as W · rbox/2, where the integer W denotes the winding number. In this approach, all non-winding paths are ignored, even if they participate in permutations because their contribution is negligible at large values of rbox (i.e., in the thermodynamic limit). Based on these considerations, the superfluid fraction can be defined as48
where the winding number W counts how often a path wraps around the periodic boundary conditions along a given direction times the length of the simulation cell in that direction. At variance with the above-mentioned formula for finite bosonic systems, this estimator is rigorous for computing the bulk superfluid fraction within periodic supercells once corrected for finite-size effects.145 We note in passing that the topological concept of winding numbers can also be used to numerically compute path integrals on restricted spaces such as encountered when quantizing rigid-body rotational motion.146–149
B. Path integral molecular dynamics for solute molecules
As already explained, PIMD is used in our hybrid HPIMD/MC approach in order to sample the configuration space of the solute molecule that is embedded in a bosonic helium environment because MD is well suited for the solute subsystem and has several advantages over MC. In particular, the configuration space is sampled according to Newton’s equations of motion due to the classical isomorphism, and thus, no specialized moves need to be designed. In addition, thermostats to sample the canonical (NVT) ensemble, such as massive Nosé–Hoover chain thermostats,92,150 or stochastic friction-based techniques, such as the Path Integral Langevin Equation (PILE)151 thermostat, further improve the sampling efficiency, as explicitly shown for doped helium clusters at low temperatures.152 Moreover, time dependent properties can, in principle, be approximately computed in the framework of PIMD using variants such as Centroid MD (CMD)153 or Ring Polymer MD (RPMD)154 if the temperature is not too low. Moreover, both methods are plagued with serious deficiencies (namely, the curvature problem in CMD and spurious resonances for RPMD) when it comes to computing vibrational densities of states or infrared spectra,155,156 which can be partially cured in the case of RPMD.157–160
However, performing path integral simulations at very low, but finite, temperatures is a computational challenge as the number of beads P required to converge the discretized path integral numerically becomes prohibitively high. For this reason, most simulations of clusters or molecules embedded in 4He have been performed using either a completely frozen solute molecule or a rigid-rotor approximation for the description of the impurity. These simulations are thus completely neglecting the vibrational degrees of freedom, which requires the highest number of beads to converge. While this enormous simplification is usually a valid approximation for quasi-rigid molecules, such as SF6, OCS, or CH4, to name but a few such cases, it breaks down for more flexible molecules, such as the hydronium (H3O+) cation subject to pronounced umbrella tunneling161 or the protonated methane () molecule103,137 in view of its hydrogen scrambling. In such cases, one needs to perform converged path integral simulations while accounting for all degrees of freedom, which is a nearly impossible task at ultra-low temperatures on the order of 1 K, if standard techniques based on the so-called primitive (lowest-order Trotter) approximation of the thermal path integral action are used. As we explicitly demonstrated by converged benchmark calculations using the PILE thermostat151 for some representative showcase molecules,107,108 namely, , CH4 and , Trotter numbers of about P ≈ 10 000 are required to converge intra-molecular observables at about 1 K.
Various methods have been designed in order to accelerate the convergence of path integrals with respect to the number of discretization replica P. This includes the mentioned pair density approximation used to describe the He⋯He interactions, which is based on approximating the high temperature N-body density matrix as a product of pair density matrices.48 While this approach is perfectly suited for the pairwise interaction of point-like particles used to describe the He⋯He interactions based on the tabulated or splined pair functions as explained earlier, it is not easily applicable for more complex cases, such as flexible or even reactive solute species and the corresponding solute⋯He interactions due to the curse of dimensionality. Fortunately, other techniques based on different approaches such as higher-order Trotter discretization schemes to go beyond the standard primitive approximation,162–167 perturbation theory,168 or colored noise thermostats169,170 can also be used to effectively enhance the convergence of the path integral as a function of increasing discretization number P. We note in passing that in PIMD (or PIMC), fully converged values of different quantities, such as energies or structural properties, can, in general, be estimated by extrapolation to the P → ∞ limit71,152,171 by making use of the asymptotic nature of the Trotter discretization.172–174 This approach however requires us to perform calculations for different numbers of replica in order to extrapolate.
Among these methods, the use of colored noise thermostats has been shown to significantly accelerate the convergence of PIMD simulations175–178 even down to ultra-low temperatures and for chemically complex and highly anharmonic molecular systems.107,108 The main idea behind these colored noise thermostats, such as PIGLET175,176 and PIQTB,177 is to achieve a frequency dependent thermalization so that the energy distribution among the different normal modes of angular frequency ω is given by the quantum harmonic expression,
In practice, the target distribution of Eq. (7) is imposed through the use of a Langevin thermostat with a random force, which is not a white noise as in the classical case, but a colored noise. The power spectral density of this colored noise is directly proportional to the quantum energy distribution of Eq. (7). These thermostats are thus constructed in the spirit of a thermal bath consisting of quantum harmonic oscillators and are able to approximately include part of the NQEs, such as the zero point energy and quantum delocalization, in standard MD simulations of classical point particles.169,170,179–188
It has also been recognized that these colored noise thermostats can be used in the context of PIMD simulations as a way to accelerate the convergence with respect to the number of beads.175–177 In this case, the quantum-mechanical thermal bath is applied to the whole ring polymer, and the target energy distribution thus needs to be modified. In the path integral framework, one can write the average energy of a harmonic normal mode of angular frequency ω as
with being the average energy of the kth ring polymer normal mode of angular frequency and ωP = P/βℏ being the angular frequency of the quantum springs. As before, the colored noise is used to enforce the proper quantum energy distribution in the harmonic approximation but ensures in the path integral framework that the average energy [Eq. (8)] is given by the correct quantum expression of Eq. (7) for any value of P. Thus, the energy distribution that needs to be enforced on the ring polymer normal modes must be solution of the following equation:
with m being the mass and γ being the friction coefficient of the Langevin thermostat. Within this framework, part of the NQEs is then already approximately included through the thermostat so that fewer replicas are necessary to converge toward the exact quantum limit P → ∞.
The striking computational advantage when using colored noise thermostats is illustrated with the help of Fig. 2, where the convergence with respect to the number of beads of the distribution functions of H–H distances for the quasi-rigid methane (CH4) molecule and its highly fluxional cousin protonated methane () is presented. Clearly, the use of colored noise thermostats, here PIQTB,177 allows for a significant acceleration of the convergence with respect to P compared to standard PIMD simulations (using PILE151 in the present case). In particular, it should be noted that the distributions obtained using a small number of beads (such as P = 32) are already very close to the fully converged ones even in the difficult case of the highly fluxional protonated methane molecule. We note in passing that the PIQTB177 and PIGLET175,176 approaches yield quite similar convergence acceleration at ultra-low temperatures, as demonstrated recently in Ref. 108.
Having shown that colored noise thermostatting greatly improves the path integral convergence of structural properties by adding artificial fluctuations according to quantum harmonic oscillators, one might wonder if these do not distort the partly anisotropic quantum delocalization properties of nuclei within molecules. As demonstrated in Fig. 3 of Ref. 108, it is reassuring to find that this is not the case, although colored noise thermostats tend to overestimate the particle delocalization if too few beads are used, yet they convincingly reproduce anisotropic NQEs.
Thus, the use of colored noise thermostats allows one to perform converged PIMD simulations of molecular systems even down to ultra-low temperatures on the order of one K while accounting for the full flexibility of the molecule including notably large-amplitude motion and fluxionality.
C. Tackling the helium⋯solute coupling
As introduced in Secs. III A and III B, ultra-low temperatures on the order of 1 K require many thousands of replicas to yield a converged representation of the path integral in the primitive approximation. At the same time, the pair density approximation,48,138 utilized for PIMC simulations of helium, provides the converged results using one to two orders of magnitude fewer beads, while Trotter convergence in PIMD simulations can be accelerated via colored noise thermostats by roughly the same factor of 10–100 as we have just seen. However, it is still possible that the solute part, which features usually very high frequency vibrations that are harder to converge, requires a higher discretization of the path integral than the helium environment when using the pair density approximation. If these distinctly different acceleration approaches shall be combined, as required in our HPIMD/MC scheme, one possibly has to couple different numbers of path integral replica for the helium environment (using the pair density approximation in PIMC) and the hosted impurity (using some colored noise thermostatting in PIMD). In recent years, different approaches for such a task have been developed, ranging from ring polymer contraction190–192 to multiple-time stepping in imaginary time.193 In the following, such ideas are transferred to our coupled HPIMD/MC approach in order to achieve a converged quantum description of the solute while coupling it efficiently and most accurately to the quantum solvent at ultra-low temperatures. Our approach relies mainly on the multiple-time stepping technique introduced in Ref. 193.
In order to derive a meaningful solute⋯solvent bead coupling, let us start by focusing on the Hamiltonian of the MC and MD subsystem individually. During the PIMD propagation of the solute in the fixed environment of the solvent, the solute is fully described by the solute Hamiltonian Ĥmol and its coupling via the interaction potential with the solvent Ûint so that ĤPIMD = Ĥmol + Ûint. Here, it is important to realize that the He⋯solute interaction is relatively weak compared to the intramolecular solute interactions. Such energetically low-lying modes of a system exhibit less zero point energy than higher modes. As a consequence, it is expected that the He⋯solute interaction can be converged with fewer path integral beads than the solute. Having this in mind, one can derive a two-step Trotter splitting of the density matrix as shown in the following.
A first Trotter splitting with a low discretization Plow is introduced to split off the He⋯solute interactions (given that the He⋯He interactions are represented using the pair density matrix approximation where rather small Trotter numbers ensure convergence). In the next step, the remaining matrix element involving only the solute can be further discretized using a second Trotter splitting resulting in a larger Trotter number P = Plow · Pstride, which is a multiple of Plow since Pstride > 1. One thus arrives at an effective potential of the PIMD subsystem of the form
where Pstride defines the so-called striding length.
Consequently, the He⋯solute interaction is only applied to every Pstrideth bead of the molecular solute; a scheme of the resulting coupling approach is presented in Fig. 3. As apparent from Eq. (11), these beads, which are directly coupled to the quantum solvent, experience this interaction Pstride-times stronger and distribute it via the harmonic couplings within the solute’s ring polymer to the adjacent beads. This also indicates the limitations of this approach: Once the instantaneous forces due to the He–solute interaction times Pstride are too large, the Trotter splitting will break down. As outlined above, the approximation is therefore only valid in the limit of sufficiently weak interactions and thus small Uint.
Let us next focus on the advantages of this coupling via striding compared to other approaches, especially at very low temperatures. In particular, ring polymer contraction and an averaging scheme are considered as other approaches to couple two systems with different numbers of beads. These methods are not derived here, but details for the contraction approach can, for example, be found in Ref. 190. The main idea of the first approach is to reduce the number of replicas of the path integral ring polymer by discarding the highest frequencies of the system in the normal mode representation. Averaging achieves a reduction in the number of beads by a simple geometric average of adjacent beads. In particular, the normal mode contraction approach has been shown to be a powerful tool to boost the performance of path integral simulations under ambient conditions191,192 where rather low bead numbers suffice to ensure convergence.
However, for the purpose of studying systems at ultra-low temperatures using the HPIMD/MC method, multiple-time stepping in imaginary time has some advantages. To highlight this, it is useful to consider the PIMC subsystem in more detail. During the MC phase of our hybrid approach, the system is described by the Hamiltonian ĤPIMC that is comprised of the solvent Hamiltonian ĤHe and the coupling via the interaction potential with the solute Ûint so that ĤPIMC = ĤHe + Ûint. In this case, it is important that the solvent interacts with the solute in representative configurations of its quantum ensemble since it is kept frozen during MC sampling. However, the averaging and contraction based methods significantly alter the structural and delocalization properties of the solute. At the same time, multiple-time stepping in imaginary time does not modify the structural properties at all because it comes down to essentially skipping a certain amount of replica. In the path integral formalism, all replica are equivalent and thus provide the same structural information. The structural properties are of major interest for the correct coupling of the PIMC and PIMD subsystems within our hybrid sampling approach, but the delocalization of the PIMD systems also needs to be transferred correctly to the helium environment. Explicit tests by us have shown that multiple-time stepping in imaginary time is able to conserve the correct quantum delocalization even for large striding lengths Pstride, while both the averaging and contraction methods provide always smaller particle delocalization, in particular, at ultra-low temperatures. This clearly highlights that imaginary time striding provides the best representation of a reduced ring polymer for the purpose of coupling it with the helium environment under ultra-cold conditions.
After having confirmed that multiple-time stepping in imaginary time is a valid option for the coupling of two differently discretized subsystems, we will focus next on quantifying the remaining uncertainties of this approach. For this purpose, conventional PIMD techniques can be used to describe solute⋯He complexes with only a single He atom since this treatment is known to straightforwardly converge to the numerically exact result for P → ∞, although a very high number of beads are needed to provide such fully converged benchmarks at ultra-low temperatures. In direct comparison to this path integral sampling of the solute⋯He complex, our HPIMD/MC method is used to describe the same complex using a combination of all above-introduced techniques in terms of solute and solvent sampling, thus using a substantially lower number of replicas compared to the conventional simulations. For the sake of one-to-one comparison of the PIMD benchmark to the HPIMD/MC scheme, the identical additive intra-molecular solute and intermolecular solute⋯He coupled cluster potentials as introduced in Secs. IV B and IV C, respectively, are employed, although a global surface for the entire solute ·HeN complex could be easily generated for N = 1. This outlined validation procedure has been performed for the CH5+ · He and CH4 · He complexes at 1.67 K in order to assess fluxional and quasi-rigid solute molecules, respectively.
In Fig. 4, the H⋯He and O⋯He distance distribution functions of these two representative solute⋯He complexes are compared for the two simulation approaches. In order to investigate if the hybrid coupling between the PIMC and PIMD regimes is valid, these two distance distributions were chosen as the most important features. For both structural properties and both systems, our hybrid HPIMD/MC approach, using all introduced methods for efficient sampling at ultra-low temperatures, reproduces the conventional results with many more replica in the frame of standard PIMD simulations extremely well (still, some minor differences can be detected). The shoulder and main peak in the H⋯He distance distribution functions at smaller and larger distances, respectively, result from the different equivalent hydrogen atoms in the two systems with respect to the single He atom, as resolved in terms of the two sets of dashed lines (see the caption). These subtle features are nicely reproduced by the hybrid approach with only slight shifts in the distribution. Note that the corresponding average values differ only by about 1%, as visualized by the vertical lines.
In conclusion, this analysis demonstrates that our HPIMD/MC method for efficient treatment of the quantum convergence is able to provide an accurate description of solutes in helium for the given interaction potentials. While for a single He atom, it is, in principle, possible to resort to the global potential energy surface of the complex and thus to conventional PIMD techniques (such as PILE or massive Nosé–Hoover chain thermostatting), albeit at the price of employing very large Trotter numbers, this approach is impossible, in practice, for the simulation of larger solute⋯helium clusters including bosonic exchange involving many 4He atoms. The reviewed and introduced acceleration techniques therefore open the door to efficiently describe complexes solvated by superfluid helium and thus shed new light on these intriguing systems at ultra-low temperatures.
IV. INTRA- AND INTERMOLECULAR INTERACTIONS
A. Helium⋯helium interactions
The van der Waals interactions between helium atoms, which are among the weakest intermolecular interactions known, lead to the truly unique properties of small helium clusters194–197 and are deeply connected to the formation of the superfluid phase in liquid helium.48 For the simulation of most of these properties, a pair potential fitted to highly accurate ab initio calculations198 is sufficient to provide reliable results. This also enables the usage of the pair density matrix approximation,48 described in Sec. III A, and therefore provides a computationally efficient and accurate approach for the simulation of many particle helium systems. Note that a lot of work has been devoted to the development of more sophisticated helium⋯helium potentials based on ab initio methods,199–209 which include adiabatic, relativistic, or higher order terms in the many body expansion. However, for the most properties of interest in the present context, a sufficient agreement between computer simulation and experiment is already provided when using the long established pair potentials due to the work of Aziz and co-workers including the original empirical 1979 variant.210
In order to illustrate the quality of the results obtained with the 1995 pair potential198 that we usually employ in our simulations (combined with the pair density matrix approach using a Trotter discretization P such that P × T = 80 K and the bosonic PIMC sampling using our implementation of these techniques in the CP2k code, in particular the canonical worm sampling as defined in the Appendix), different properties of simulated bulk helium are compared to experiment in Fig. 5. First, the specific heat CV, computed by numerically differentiating the total energy with respect to the temperature, is compared to the experimental values as a function of temperature. The well-known divergence of the specific heat capacity of bulk 4He, called Lambda transition,212 when moving from the classical liquid phase to the superfluid phase of helium under ambient pressure conditions, is recovered at around 2.1–2.2 K, although finite size effects broaden the peak considerably compared to experiment even when using 512 4He atoms (see Ref. 131 for the pioneering work and Ref. 145 for a finite size scaling study of the Lambda transition). Last but not least, the superfluid fraction fs is also shown as a function of temperature in the bottom panel of Fig. 5. Again, the experimental behavior is well described in the simulations, and the phase transition close to the experimental value, 2.17 K, is indeed recovered, if the number of atoms is sufficiently large. We stress that we use a truncated octahedral supercell in our implementation in order to minimize the number of atoms located beyond the largest non-periodic spherical cutoff and also compute the winding number and thus the superfluid density accordingly, as detailed in Sec. II E of Ref. 100.
As demonstrated based on the physical observables collected in Fig. 5, the weak intermolecular interactions are faithfully described by the standard He⋯He pair potential used198 and the bosonic quantum nature of 4He is also nicely reproduced by our PIMC simulations in agreement with the pioneering work.131,144
What remains challenging is the description of the intra-molecular solute interactions, i.e., the full-dimensional global PES of the solvated species, as well as the solute–helium interactions, where small errors are sufficient to tip the balance with respect to the delicate helium–helium interactions. This demands both very accurate correlated electronic structure calculations and their accurate representation in terms of computationally efficient potentials that provide analytical forces in addition. In order to address these remaining challenges, our most recent developments based on machine learning potentials are presented in Secs. IV B and IV C.
B. Intra-solute interactions: Full-dimensional global potential energy surfaces
In most theoretical studies of helium (micro-)solvation, the molecular solute species has been treated either as a classical point-like particle system62–64,67,83–85,213 or as a quantum rigid rotor,41,42,66,68–72,74,75,77–82 thus excluding vibrational degrees of freedom. These studies have led to a wealth of pioneering insights into quantum solvation as well as the general understanding of superfluidity and its manifestations at the microscopic scale. In particular, for quasi-rigid molecules, such as SF6, OCS, or CO2, to name but a few, it is expected to be an excellent approximation to restrict the description to rotational and translational degrees of freedom by freezing the molecular skeleton, thus neglecting any small-amplitude intra-molecular motion. At the same time, these approaches prevent a more detailed analysis of the influence of helium solvation on the solute molecule beyond rotational dynamics and are severe simplifications in the case of more flexible species, in particular so for those undergoing large-amplitude motion such as fluxional molecules. Other cases are those with soft inter-molecular degrees of freedom such as hydrogen bonded dimers where an interesting impact of helium microsolvation on tunneling splittings has been discovered,65,73 albeit at 0 K and while keeping the interacting monomers frozen.
Including intra-solute motion and thus the internal flexibility of molecules, complexes, or clusters forces one to deal with full-dimensional potential energy surfaces of these species (which finally must also be coupled to helium via solute⋯He interactions for all possible solute configurations and, thus, at each point of the respective PES, as detailed in Sec. IV C). The situation gets even more complicated if the solute species can undergo chemical reactions. For small gas phase molecules, vibrational degrees of freedom have usually been included by means of specialized parameterizations of physical or purely mathematical functional forms to represent the PES.214–218 This task can quickly become cumbersome and difficult, which easily hampers any routine exploration of chemically complex species including their reactions. In the realm of quantum simulations at finite temperatures, this problem has been solved long ago by introducing ab initio path integrals (AIPI)90–94 where the electronic structure problem is solved on-the-fly while sampling the path integral using molecular dynamics techniques. This avoids the pre-parameterization of typically high-dimensional potential energy surfaces or interaction potentials, yet the price to pay is using computationally efficient electronic structure theory, which means resorting to density functional theory in practice. As is well known, the quality of these density-based methods very much depends on the specific case and on the choice of the functional, although excellent descriptions can be achieved as demonstrated repeatedly, for instance, for protonated methane219–222 based on extensive benchmarking.
In the last few years, new approaches based on machine learning and thus on general functional forms have been established as an alternative to the more traditional approaches, as reviewed, for instance, in Refs. 223–227. Among many distinct approaches,228–239 high-dimensional neural network potentials (NNPs)228 combined with atom centered symmetry functions240 have been shown to reliably and accurately reproduce electronic structure based reference calculations. Moreover, they allow for a high level of automation of the fitting process, as reviewed in detail in Refs. 224, 225, 241, and 242. This enables the rapid development of highly accurate potential energy surfaces for inter- and intra-molecular interactions.243–248 Recently, such approaches have been demonstrated to even grant efficient access to custom-made high-dimensional interaction potentials at a quality level108–110,249,250 that is broadly considered to be the “gold standard” in quantum chemistry, namely, coupled-cluster theory using singles, doubles and perturbative triple excitations [CCSD(T)].251
In order to illustrate in the present perspective the power of these modern machine-learned approaches as applied to quantum simulations, the first machine-learned CCSD(T) PIMD simulation108 being published in 2018, our recently established adaptive automated fitting procedure110 (that has been employed earlier to generate the NNP of bare H5O2+, which has been used previously in Ref. 108) is applied here to four different systems. All four species are interesting targets for helium solvation, namely, the methane molecule, its fluxional cousin protonated methane, as well as the hydronium and Zundel cations since they represent a most challenging mix of quasi-rigidity, low-dimensional large-amplitude motion, and high-dimensional fluxionality. Their global PES are described by high-dimensional NNPs,224,225,241,242 which are trained to energies obtained from explicitly correlated coupled cluster calculations. These include single, double, and perturbative triple excitations in conjunction with the F12a correlation factor252 and scaled triple contributions252 in the aug-cc-pVTZ basis set [i.e., CCSD(T*)-F12a/AVTZ], as implemented in the Molpro package.253 This reference method has been demonstrated to yield results very close to the complete basis set (CBS) limit252 due to using this explicitly correlated triple-zeta basis.
The generation of our CCSD(T*)-F12a/AVTZ quality NNPs is performed using an automated fitting procedure that identifies by itself the most representative training points110 by drawing on the established properties of neural network potentials.224 Key to our procedure is to perform the least number of computationally demanding CCSD(T*)-F12a/AVTZ electronic structure calculations by an adaptive sampling strategy [while not using forces that are increasingly demanding to compute numerically or analytically for large N-atom systems when using consistently the CCSD(T*)-F12a/AVTZ methodology]. Our strategy targets specifically these parts of the configuration space of the specific system, which are relevant for the envisaged usage of the NNP in computer simulations. We stress that this subspace is certainly different for low-temperature quantum simulations close to the quantum ground state compared to thermally activated dissociation or reactive scattering calculations, which must be carefully considered when using a certain NNP. In other words, no NNP should ever be used for purposes for which it has not been explicitly trained.
The main steps of our automated iterative procedure110 are the following:
A large number of configurations are generated by classical MD or PIMD simulations at different temperatures, which are initially performed on-the-fly using general purpose, efficient, and easy-to-handle DFT-based methods. Once a dynamically stable NNP becomes available, this NNP is used to continue these MD or PIMD simulations, thus greatly improving sampling efficiency.
Relevant configurations are selected specifically either by comparing two (or more) different NNPs, which feature large differences in regions that are still under-represented in the MD or PIMD sampling procedure, or by detecting situations where the configuration space of the existing training set is left. Note that for the very first iteration, a small number of randomly chosen configurations on the order of as few as 10–20 configurations are used to seed the procedure.
Reference energies (at the desired coupled cluster level) are computed exclusively for those specifically selected configurations (rather than using random or otherwise sampled configurations from the above-mentioned MD or PIMD runs), which are then added to the dataset.
New NNPs are fitted to the improved dataset, and the cycle starts over at the stage of using NNP MD or PIMD simulations.
These four steps are repeated until a sufficiently accurate neural network potential is obtained. Our procedure is described in more detail in Ref. 110 and grows the training set for which the coupled cluster calculations are performed in an optimal fashion, thus allowing one to reduce as much as possible the number of these expensive electronic structure reference calculations. In practice, the training sets used for all systems considered in this perspective are composed of only on the order of 10 000 points and the associated root mean square errors (RMSEs) are typically well below about 0.1 kJ mol−1.
In order to explicitly validate the accuracy of the resulting NNPs, short quantum PIMD trajectories at 1.67 K have been generated using the NNPs. Afterward, single-point coupled cluster calculations have been carried out along these trajectories at each and every individual PIMD sampling step. The resulting comparison of the energies along the trajectories is depicted in Fig. 6 for our four representative solute species, namely, CH4, , H3O+, and from top to bottom. For all four systems, almost perfect agreement between the NNP and the explicit coupled cluster reference is observed, thus highlighting the quality of our approach (note that only a few red circles are added to visualize the CC lines that perfectly overlay the NNP lines on the scale set by the respective energy fluctuations). At the same time, the evaluation of the NNPs is many orders of magnitude faster than the coupled cluster reference, granting access to the converged electronic structure results at negligible cost.
Moreover, the individual NNPs can be developed very efficiently due to the automated procedure (so far successfully applied up to the protonated water tetramer in Ref. 110). This procedure iteratively improves the training set while keeping the demanding reference calculations to a minimum. This opens the door for routine exploration of medium sized molecules and molecular complexes under ultra-cold conditions including all their degrees of freedom as well as the quantum nature of the nuclei while employing an essentially converged description of the interactions.
Transcending this specific class of applications, we are confident that our adaptive automated NNP fitting procedure110 can also be used to generate global potential energy surfaces, including analytical forces and Hessians, for distinctly different purposes than (PI)MD, such as semiclassical wavepacket propagation or full quantum dynamics.
C. Solute⋯helium interactions
In order to investigate the impact of molecular motion on helium solvation and superfluidity (and vice versa), an accurate description of the intermolecular interactions is essential. Indeed, these solute⋯helium interactions are usually rather weak and can be comparable to the He⋯He interactions, and the competition between these two different interactions gives rise to complex solvation patterns that can change drastically with the number of 4He atoms. While wavefunction based ab initio methods such as coupled cluster theory are able to provide the required level of accuracy, they usually are computationally too demanding to be used “on-the-fly” in path integral simulations with notable exceptions.254 Indeed, most studies of doped helium clusters have thus used physically motivated analytical interaction potentials fitted to reproduce reference data usually obtained by accurate ab initio methods.41,62,64,66–68,70,72–76,78–85,103,137,213 However, this requires the design of specific interaction potentials for each system under study, which can quickly become a time consuming and complex task. This is especially complicated if one wants to rigorously reproduce the fine structure of helium solvation from small clusters composed of a few 4He atoms up to much larger ones, including possibly periodic bulk while still accounting for the full large-amplitude flexibility of the solute, including chemical reactions.105,106,255 For this reason, previous studies have been focused on a relatively small set of dopant molecules for which accurate potentials are available, such as SF6, OCS, or CO2, and among others.104
As for the above-mentioned intra-solute interactions, machine learning techniques such as neural networks can also be used to accurately describe the solute⋯helium interactions in an efficient way. Here, we use the same family of high dimensional NNPs224,228 combined with atom centered symmetry functions240 to describe the atomic environment, as introduced before for the description of the solute potential energy surface. Because the typical values of the solute⋯helium interaction energies are usually completely different from the energy range associated with the inter-molecular interactions of the solute, it is beneficial to separate between inter- and intra-molecular interactions. Consequently, we use two different NNPs to describe these two types of interactions. Moreover, this separation is ideally suited for our hybrid HPIMD/MC sampling approach presented earlier, which relies on a different sampling strategy for the solute and the solvent. Now, the NNP is trained to reproduce the interaction energy between one helium atom and the solute molecule, thus accounting for two-body terms in a many-body expansion. The total interaction energy is then obtained in a pairwise fashion by adding all solute⋯helium two-body terms in the system. Thus, three-body and higher order terms are neglected in our approach. This has been shown to be an excellent approximation in the case of the He interactions for which three-body terms have been shown by extrapolated coupled cluster calculations to account for less than 1% of the total interaction energy.106
The interaction NNPs are obtained in a very similar fashion as just described for the intra-solute potential energy surfaces using the adaptive automated fitting procedure.110 However, small modifications of the procedure have been introduced109 for the present task as follows. The initialization of the procedure is extended,109 since it is not possible to perform reasonably accurate on-the-fly electronic structure based PIMD simulations of solute⋯helium complexes in view of the significant dispersion (also known as London or van der Waals) contributions to these interactions, which are difficult to faithfully include in DFT based AIPI simulations. First of all, the solute configurations are not sampled from previous MD or PIMD trajectories, but randomly extracted from the full dataset of reference coupled cluster calculations, which has been generated previously to fit the solute NNP. Second, a large number of positions for the helium atom are generated on randomly perturbed grids around these solute configurations as elaborated in the original article,109 where an unbiased simple sampling MC procedure is applied on the grid to stochastically select a computationally manageable subset of points. At this stage, a short range cutoff is used in addition to avoid generating helium positions that are too close to the solute molecule, which would result in very high (repulsive) energies that could distort the NNP fit with no gains. Similarly, a long range cutoff is used to avoid generating the helium position too far away from the solute species for which the interaction energy approaches zero.
Akin to what is done in the case of the intra-solute interactions, as explained in Sec. IV B, two different NNPs are used to predict the interaction energies associated with the generated configurations and only the most relevant ones are selected based on the difference of predicted energies. However, the absolute values of the interaction energies in the repulsive (short range) region are significantly higher than the ones for helium positions far away from the molecule. Consequently, the absolute difference between the two NNP predictions, ΔE, is also significantly higher in the repulsive region, and these short range configurations are over-represented if only the configurations with the highest ΔE are selected. To circumvent this problem, a large number of configurations associated with the highest ΔE are extracted and then a Metropolis importance sampling MC procedure (instead of simple sampling MC) is used to select configurations with a probability distribution p(E) ∝ exp(−βE). The value of the β coefficient is not the inverse thermal energy 1/kBT as usual but is also randomly chosen within meaningful bounds at each new iteration to automatically ensure a satisfactory sampling of helium configurations, which includes both short range and long range configurations. The β parameter bounds are adaptively chosen up to four times some average interaction strength of the solute species and the helium atom. Only after new configurations have been selected by this pre-screening technique, the associated reference interaction energies are computed using the standard “supermolecule” approach based on coupled cluster calculations. Here, the same electronic structure method as before for the solute molecule is used, namely, CCSD(T*)-F12a/AVTZ, but with the counterpoise correction applied to correct for the basis set superposition error256 that remains even when computing interaction energies using explicitly correlated (F12) wavefunctions with a triple-zeta basis set.
This iterative fitting procedure, which is critically based on random sampling protocols to maximize ignorance as well as automation, has been shown to yield very accurate NNPs to describe the helium solvation of small protonated water clusters using datasets composed of only about 50 000 configurations for which the coupled cluster calculations need to be performed.109 In the meantime, we have fitted several NNPs using this procedure for different dopants, namely, methane (CH4), protonated methane (), hydronium (H3O+), and the Zundel complex () with training and test RMSEs that are lower than 0.3 kJ mol−1 even for the complex case of the highly flexible molecule (while always using datasets containing about 50 000 configurations, which already includes the sampling of sufficiently many intra-molecular configurations of the solute species due to its vibrational displacements in view of the choice of the solute configurations as described above).
The ability of these NNPs to describe the solute⋯helium interactions is illustrated in Fig. 7, which shows helium densities in the three-dimensional Cartesian space as obtained from path integral simulations at T = 1.67 K for different doped helium clusters, namely, CH4 · He15, , H3O+ · He6, and from top to bottom. For these particular simulations designed to benchmark the accuracy of the solute⋯helium interactions as given by the corresponding interaction NNP compared to the coupled cluster reference, the molecule is kept frozen in space (specified in the caption) and the solute⋯helium interaction energies needed to perform PIMC in the Cartesian continuum are interpolated from values that have been precomputed on rather dense grids using either the NNP or the reference method (being the counterpoise-corrected CCSD(T*)-F12a/AVTZ calculations for these interactions). We note in passing that these grids contain typically 250 000 points at which the CC calculations have been carried out (which is one order of magnitude more than what is required to parameterize the interaction NNP), but solely for the purpose of carefully benchmarking the solute⋯helium NNPs “at work.” This procedure allows for visual comparison of helium densities obtained using the NNP and the coupled cluster method at the level of the relevant orientationally resolved spatial modulations due to the building up of quantum solvation shells. From the densities shown in Fig. 7, it is clear that the NNPs are able to very accurately reproduce the helium densities in all cases. Overall, this technique permits one to conduct path integral simulations of helium solvation at coupled cluster accuracy and without restriction on the number of 4He atoms.
V. CONCLUSIONS AND OUTLOOK
In summary, the fully atomistic quantum simulation of molecules as well as molecular complexes and clusters solvated by superfluid helium has made significant steps forward due to many methodological improvements in the last few years. One of these developments is the hybrid (ab initio) path integral molecular dynamics/bosonic path integral Monte Carlo approach (HPIMD/MC)100 that has been reviewed in detail in this perspective. This approach enables the description of fluxional and even reactive systems in bosonic helium while accounting for their quantum nature.101 The combination of (i) colored noise path integral thermostats175–177 for the simulation of quantum solutes at ultra-low temperatures107,108 with (ii) the pair density approximation and the sampling of the exchange permutations131,134,135,138 for the bosonic quantum solvent in the framework of (iii) the HPIMD/MC technique100 provides an essentially converged quantum statistical description of dopants in superfluid helium. A striding approach to couple the solute with the solvent more accurately is introduced herein and shown to allow for the convergence of the molecular dynamics and Monte Carlo processes using different path integral methods and thus different levels of Trotter discretization to sample the entire path integral that represents the full solute–solvent system in the context of hybrid PIMD/MC sampling. Moreover, a canonical version of the so-called worm algorithm to efficiently sample bosonic permutations even for strongly interacting solute⋯helium systems containing only a few 4He atoms is presented herein and validated using as many as 512 bosonic atoms.
As recently demonstrated for the protonated cousin of methane, , the HPIMD/MC technique can provide unique insights into the delicate couplings between large-amplitude motion leading to molecular fluxionality137 and bosonic exchange statistics in the helium environment.103 However, these recent studies were restricted to only a few helium atoms in the microsolvation shell of the solute due to the delicate nature of the solute⋯helium interactions that render the fitting of traditional pairwise interaction potentials difficult at best.105,106
Recent developments based on neural network potentials (NNPs)109,110 allow us to systematically overcome these previous shortcomings. By means of adaptive automated improvement, these NNPs provide a straightforward way for the investigation of complex molecular systems including clusters in superfluid helium. This approach gives access to interactions described at the essentially converged coupled cluster level, namely, CCSD(T*)-F12a/aug-cc-pVTZ including the counterpoise correction when computing interaction energies. With these techniques now at hand, it will be possible to study the entire progression from small helium clusters to helium droplets to periodic bulk superfluid solutions that host fluxional and even reactive molecular species using a quantum simulation approach that is essentially converged at all levels. We hope that this will stimulate new experimental work toward applying helium messenger and nanodroplet isolation spectroscopies to those molecular species that are subject to large-amplitude motion, such as fluxional or floppy molecules, and to address cryochemical reactions beyond what has been achieved in the last decade.20,21
In the context of rotating impurities in superfluid helium, as amply referred to in the previous text, it will be especially interesting to understand the salient couplings between the superfluid response and the dynamics of fluctuating charged molecules that feature much stronger solute⋯helium interactions compared to neutral species. Another possible direction for future work is the study of other bosonic quantum solvents that can also potentially feature “manifestations of superfluid behavior,” such as para-H2 for which both pure140–142,257–259 and doped clusters33,139,260–263 have already been the subject of many insightful theoretical studies as comprehensively reviewed not long ago.34
Let us now address some of the major computational challenges awaiting progress in the future. Accounting for the bosonic quantum nature of helium is already a very challenging task, but fermionic exchange can also become important at ultra-low temperatures, if light particles such as hydrogen atoms are present in the solute—not to mention the use of 3He as a fermionic quantum solvent.37 Although it has recently been demonstrated to be of minor importance for the quantum structural properties specifically of protonated methane (see the detailed analysis in the supplementary material of Ref. 103), it remains unclear more generally speaking how the fermionic nature of the nuclei would influence different properties such as the entangled264 rovibrational dynamics being intricately coupled to the nuclear spin statistics within the solute species.264–266 However, accounting for fermionic exchange in finite-temperature path integral simulations of sufficiently large and strongly interacting many-body systems is notoriously difficult and is even today one of the unsolved challenges in the field of computer simulations, although many attempts have been made for decades trying to overcome the difficulty. More importantly, for a more direct comparison with spectroscopic experiments, it would be beneficial to compute time dependent observables, which unfortunately is not directly accessible in imaginary time path integral simulations. Specifically in the context of rigid-body rotation of molecules in bosonic 4He or para-H2 environments, several dedicated methods have been designed to extract spectroscopic information from imaginary time correlation functions.33,34,66,68,70,71,77,261,263 More direct approaches have also been designed in the framework of path integral molecular dynamics to generate approximate quantum dynamics at finite temperatures, thus providing access to time-correlation functions as noted above in passing. Yet, they seem to break down when reaching cryogenic temperatures due to different intrinsic problems. Clearly, computing the accurate quantum dynamical properties of fluxional and possibly reactive molecular solutes in bosonic solvents at cryogenic temperatures, in particular rotation-resolved infrared spectra as commonly recorded in helium tagging photodissociation and HENDI action spectroscopy, remains another grand challenge.
F.B. and C.S. contributed equally to this work.
The data that support the findings of this study are available from the corresponding author upon reasonable request.
This research was partially supported by Grant No. MA 1547/19 to D.M. and is also part of the Cluster of Excellence “RESOLV”: funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy (Grant No. EXC 2033 – 390677874). In addition, we would like to thank Professor Joachim Werner for providing the auxiliary basis set functions of He needed to evaluate the F12a correlation factor within the Molpro package. C.S. acknowledges partial financial support from the Studienstiftung des Deutschen Volkes and the Verband der Chemischen Industrie. The computational resources were provided by HPC@ZEMOS, HPC-RESOLV, and BOVILAB@RUB.
There are no conflicts to declare.
APPENDIX: CANONICAL WORM ALGORITHM TO SAMPLE BOSONIC EXCHANGE PERMUTATIONS IN PERIODIC AND FINITE SYSTEMS
In the following, the details of our implementation of the canonical worm algorithm are presented (see the main text for general background and panels B to G of Fig. 1 for its pictorial presentation). In order to formulate a canonical version of the worm algorithm that is well integrated into our hybrid HPIMD/MC simulation approach, some restrictions on the set of allowed MC moves are required. First of all, the moves obviously need to conserve the number of particles. Therefore, the “remove” and “insert” moves of the original grand canonical implementation134,135 are entirely excluded. Second, as a result thereof, the chemical potential μ has to be removed from all acceptance probabilities of the different moves. To achieve this, the sets of moves from Refs. 134 and 135 are combined in such a way that the number of ring polymer springs is conserved when proposing a new configuration, thereby removing the chemical potential from the acceptance probability. The precise nature of the combined moves used in our canonical worm algorithm is described later in detail, and we refer to Ref. 136 for a previous version of canonical worm sampling.
The grand canonical formulation propagates through the configuration space solely by performing updates on the worm itself and, thus, rarely resides in a closed state. In order to utilize the worm algorithm in our HPIMD/MC scheme, it is necessary to let every MC cycle end in a closed configuration (since, otherwise, no PIMD propagation step involving that path configuration would be possible). If this is not the case, the MC cycle is automatically extended by a fixed number of steps. However, such simple prolongation of the MC sweeps is not efficient by itself. To increase the probability of ending in a closed configuration, two additional moves are introduced, which can sample and thus generate closed states. Therefore, it is not required to immediately open the worm after it has been closed as done in the grand canonical version. These steps can also be applied during the open state and sample parts of the system not containing the head or tail of the worm to increase configuration space sampling. In general, all moves are accepted with a probability of
except for the open/close and swap moves for which the correct acceptance probabilities are explicitly given at a later stage. Note that ΔS corresponds to the change of action associated with the MC move but not including the change in the associated free particle action, since this is taken care of by the way the MC moves are proposed and that the position of a bead j of atom i is thereby denoted as . We treat all beads to be identical, and therefore, an updated imaginary time section can be formally distributed among two atoms, but for the sake of simplicity, our description here considers those sections that do not extend beyond the Pth bead.
1. Centroid move
The traditional centroid move [Fig. 1(b)] is the first update, which is additionally introduced. To perform a centroid update, an atom is selected randomly and the exchange cycle it belongs to (in open configurations also the head/tail bead) is translated by a random vector, where the components are drawn independently from an uniform distribution [−Δmax, Δmax]. As shown in Fig. 1(b), the change in action depends on all beads of all atoms from the selected permutation cycle, and this update is therefore computationally demanding. Note that the change in the kinetic energy contribution to the action is zero because the kinetic and potential energies can be evaluated independently.
2. Staging move
The well-known staging move [Fig. 1(c)] is the second update, which is not part of the original grand canonical formulation of the worm algorithm. It performs a local update on any imaginary time section as long as it does not contain the head or tail bead. Thereby, parts of the system are propagated that are not influenced by the worm. An atom and a bead index (s) are selected randomly and are used as the first anchor bead. A second anchor bead (e) is chosen by moving from the first anchor l + 1 beads along the imaginary time axis. Note that the second anchor can be located on a different atom, depending on the start anchor bead index and the permutation state of the system. The l beads between s and e are recursively resampled according to
where is a vector of random numbers with its elements drawn from a Gaussian distribution with zero mean and a variance of
and τ = β/P here and in what follows. This update is local, and the action difference for the acceptance criterion is therefore only needed for the l resampled beads as well as the associated harmonic springs. This move is immediately rejected if this update is performed in the open state, and one of the resampled beads (not the anchor beads) happens to be the head/tail bead of the worm. Note that the staging move alone enables the sampling of distinguishable particles, but, in practice, the centroid move from Sec. 1 of the Appendix can additionally enhance sampling efficiency especially for very weakly bound systems in the limit of small cluster sizes.
3. Open/close move
The open move [Fig. 1(d)] is adapted from the grand canonical worm algorithm and modified to fulfill the requirements of the canonical version. This move is the only way to switch from sampling the diagonal elements of the density matrix to the off-diagonal components. An atom (i) and one of its beads (s) are selected randomly to serve as an anchor bead. After that, a second bead is chosen by moving along imaginary time by l beads, which will be the tail bead (t) after the move. Similarly as in the grand canonical version, the beads between s and t are removed, but this imposes a problem in the canonical version, since it requires the inclusion of the chemical potential into the acceptance probability. Therefore, the move is combined with an add move of the original algorithm such that the chemical potential μ cancels. A set of l new beads is constructed starting at the bead with index s + 1 by drawing the new coordinates from a Gaussian distribution centered around bead s with a variance of
The bead number s + 2 is drawn from a Gaussian distribution with the same variance but centered around the newly created bead s + 1. This continues until bead s + l is reached, which shares the same formal time slice as bead t but is treated independently as the head bead. The proposed open configuration is accepted with a probability of
where C0 is a dimensionless parameter. It should be on the order of unity for a bulk system at its saturated vapor pressure. Then, the probability of finding the system in a closed system is roughly 50%, whereas smaller values of C0 result in a higher fraction of closed configurations. The free density matrix appearing in the denominator is given by
with λ = ℏ2/(2mHe). Note that the density matrix appearing above and below is evaluated at l τ instead of τ.
The open move is complemented by the close move, which brings the system back to the diagonal elements of the density matrix. Starting from and including the head bead, l beads are removed such that the bead labeled s is the head bead. Next, a series of l − 1 beads are constructed between bead s and t according to the staging move. The new closed configuration is accepted with a probability of
where C0 is exactly the same dimensionless parameter as for the open move.
4. Head/tail move
The staging update is able to resample any imaginary time section but relies on two anchor beads that stay fixed during the update. As a consequence, it is not possible to move the head or tail bead with the staging update since these beads can exclusively be anchor beads. Therefore, specialized moves to resample the ends of the worm are introduced as follows. The move head and move tail update [Fig. 1(e)] are identical but act on the different ends of the worm. From the head or the tail end of the worm, l beads are removed like the first step in the close move and then resampled equivalently to the second step in the open move. This effectively moves the ends of the worm, which is accepted with a probability according to Eq. (A1).
5. Crawl move
The crawl move [Fig. 1(f)] is similar to the movement of the head or tail. However, instead of reconstructing the imaginary time section that has just been deleted, the other end of the worm is extended by the same number of beads. This shifts the worm in imaginary time by ±l depending on the direction of the crawl and effectively leads to the worm crawling through the configuration space. Although this move is not strictly required, it greatly helps in preventing situations where the worm’s head and tail are trapped in separate local minima, thus making it difficult to close the worm. These are situations that are easily encountered when using a rather small number of 4He atoms around solutes that have strongly anisotropic interactions with deep potential minima that are well separated in space.
6. Swap move
The swap move [Fig. 1(g)] is the only considered MC move that is able to change the permutation state of the bosonic system by opening a closed exchange cycle and attaching it to one end of the worm. In order to enhance the acceptance probability, a swap partner is selected based on the free particle bias. Starting with a worm opening at atom w, a swap partner atom k ≠ w is selected with probability,
The corresponding reverse selection probability, which is required according to the Metropolis criterion, is given by
Next, a new path is constructed according to the staging move with and , providing the anchor beads, and accepted with a probability
This makes the bead the new head bead.