Inclusion of trial functions in the Langevin equation path integral ground state method: Application to parahydrogen clusters and their isotopologues

On the accuracy and efficiency of different methods to calculate Raman vibrational shifts of parahydrogen clusters

The Raman vibrational frequency shifts of pure parahydrogen and orthodeuterium clusters of sizes N = 4–9 are calculated using the Langevin equation path integral ground state method. The shifts are calculated using three different methods; the results obtained from each are compared to experiment and variance properties are assessed. The first method requires the direct calculation of energies from two simulations: one when the cluster is in the v = 0 vibrational state and one when the cluster has v = 1 total quantum of vibration. The shift is directly calculated from the difference in those two energies. The second method requires only a v = 0 simulation to be performed. The ground state energy is calculated as usual and the excited state energy is calculated by using the distribution of the v = 0 simulation and the ratio of the density matrices between the v = 1 state and the v = 0 state. The shift is calculated from the difference in those two energies. These first two are both exact methods. The final method is based on perturbation theory where the shift is calculated by averaging the pairwise difference potential over the pair distribution function. However, this is an approximate approach. It is found that for large enough system sizes, despite the approximations, the perturbation theory method has the strongest balance between accuracy and precision when weighing against computational cost.

Ground state chemical potential of parahydrogen clusters of size N = 21-40

We report the ground state chemical potential of parahydrogen clusters between N = 21-40 calculated using the Langevin equation Path Integral Ground State method. There has been much debate in the past whether the chemical potential size evolution in this region is jagged (indicating magic number cluster sizes) or if it is smooth (indicating some quantum melting below 1 K). We compare to previous diffusion Monte Carlo and Path Integral Ground State (PIGS) results, including very recent Variational Path Integral Molecular Dynamics (VPIMD) calculations [S. Miura, J. Chem. Phys. 148, 102333 (2018)]. We find that the ground state chemical potential is not a smooth curve and that magic number clusters are present, consistent with VPIMD and PIGS Monte Carlo results.

Path integral simulations of confined parahydrogen molecules within clathrate hydrates: Merging low temperature dynamics with the zero-temperature limit

Clathrate hydrates, or cages comprised solely of water molecules, have long been investigated as a clean storage facility for hydrogen molecules. A breakthrough occurred when hydrogen molecules were experimentally placed within a structure-II clathrate hydrate, which sparked much interest to determine their feasibility for energy storage [Mao et al., Science 297, 2247-2249 (2002)]. We use Path Integral Molecular Dynamics (PIMD) and Langevin equation Path Integral Ground State (LePIGS) for finite temperature and zero-temperature studies, respectively, to determine parahydrogen occupancy properties in the small dodecahedral (5¹²) and large hexakaidecahedral (5¹²6⁴) sized cages that comprise the structure-II unit cell. We look at energetic and structural properties of small clusters of hydrogen, treated as point-like particles, confined within each of the different sized clathrates, and treated as rigid, to determine energetic and structural properties in the zero-temperature limit. Our predicted hydrogen occupancy within these two cage sizes is consistent with previous literature values. We then calculate the energies as a function of temperature and merge the low temperature results calculated using finite temperature PIMD with the zero-temperature results using LePIGS, demonstrating that the two methods are compatible.

A path integral ground state approach for asymmetric top rotors with nuclear spin symmetry: Application to water chains

A path integral ground state (PIGS) approach for the simulation of asymmetric top rotors is presented. The method is based on Monte Carlo sampling of angular degrees of freedom. A symmetry-adapted rotational density matrix is used to account for nuclear spin statistics. To illustrate the method, ground-state properties of collections of para-water molecules confined to a one-dimensional lattice are computed. Those include energetic and structural observables. An advantage of the PIGS method is that expectation values can be obtained directly since the square of the wavefunction is sampled during a simulation. To benchmark the method, ground state energies and orientational distributions are computed using exact diagonalization for a single para-water molecule in an external field using a finite basis of symmetric top eigenfunctions. Benchmark results are also provided for N = 2 para-water molecules pinned to lattice sites at various distances to sample the crossover from hydrogen bonding to the dipole-dipole interaction regime. Excellent agreement between the PIGS results and the finite basis set calculations is observed. A thorough analysis of the convergence in terms of the imaginary time propagation length and systematic Trotter error is performed. The PIGS approach is then applied to a chain of N = 11 water molecules, and an equation of state is constructed in terms of the intermolecular separation. Ordering effects are also studied, and a transition between hydrogen bonding to dipole-dipole alignment is observed. The method is scalable and can also be applied in higher dimensions.

The Loss of Size Sensitivity in para-Hydrogen Clusters Due to the Strong Quantum Delocalization

para-Hydrogen (pH₂)_N clusters have been the focus of numerous computational studies. Originally motivated by the possibility of observing superfluidity, these studies also revealed rich and complex structural properties of (pH₂)_N. However, their structural analysis was typically limited to attempts to identify "magic number clusters" by computing their ground state energies E_N and the chemical potential μ_N = E_N-E_N-1 as a function of N. This was followed by structural analysis based on an ill-defined radial density profile. Surprisingly, however, there were remarkable discrepancies between the results reported in the literature for cluster sizes beyond approximately N = 25, and this ambiguity remained unsettled until now. In the present paper, we apply the diffusion Monte Carlo method to resolve inconsistencies in cluster sizes within the range (N = 24-28). Here, we try to avoid speculations based on the highly demanding energy calculations whose numerical accuracy harbors ambiguity. Instead, we focus on the direct and unambiguous structural analysis of the ground state wavefunctions, which supports the conclusion that the clusters are structurally the same in the size range considered. That is, there are no magic number clusters at least in the range N = 24-28, contrary to what some of the previous publications have suggested. This lack of size sensitivity of para-hydrogen clusters is a direct consequence of the strong quantum delocalization in these systems.

Diffusion Monte Carlo simulations of gas phase and adsorbed D2-(H2)n clusters

We have computed ground state energies and analyzed radial distributions for several gas phase and adsorbed D₂(H₂)_n and HD(H₂)_n clusters. An external model potential designed to mimic ionic adsorption sites inside porous materials is used [M. Mella and E. Curotto, J. Phys. Chem. A 121, 5005 (2017)]. The isotopic substitution lowers the ground state energies by the expected amount based on the mass differences when these are compared with the energies of the pure clusters in the gas phase. A similar impact is found for adsorbed aggregates. The dissociation energy of D₂ from the adsorbed clusters is always much higher than that of H₂ from both pure and doped aggregates. Radial distributions of D₂ and H₂ are compared for both the gas phase and adsorbed species. For the gas phase clusters, two types of hydrogen-hydrogen interactions are considered: one based on the assumption that rotations and translations are adiabatically decoupled and the other based on nonisotropic four-dimensional potential. In the gas phase clusters of sufficiently large size, we find the heavier isotopomer more likely to be near the center of mass. However, there is a considerable overlap among the radial distributions of the two species. For the adsorbed clusters, we invariably find the heavy isotope located closer to the attractive interaction source than H₂, and at the periphery of the aggregate, H₂ molecules being substantially excluded from the interaction with the source. This finding rationalizes the dissociation energy results. For D₂-(H₂)_n clusters with n≥12, such preference leads to the desorption of D₂ from the aggregate, a phenomenon driven by the minimization of the total energy that can be obtained by reducing the confinement of (H₂)₁₂. The same happens for (H₂)₁₃, indicating that such an effect may be quite general and impact on the absorption of quantum species inside porous materials.

Formulation of state projected centroid molecular dynamics: Microcanonical ensemble and connection to the Wigner distribution

A derivation of quantum statistical mechanics based on the concept of a Feynman path centroid is presented for the case of generalized density operators using the projected density operator formalism of Blinov and Roy [J. Chem. Phys. 115, 7822-7831 (2001)]. The resulting centroid densities, centroid symbols, and centroid correlation functions are formulated and analyzed in the context of the canonical equilibrium picture of Jang and Voth [J. Chem. Phys. 111, 2357-2370 (1999)]. The case where the density operator projects onto a particular energy eigenstate of the system is discussed, and it is shown that one can extract microcanonical dynamical information from double Kubo transformed correlation functions. It is also shown that the proposed projection operator approach can be used to formally connect the centroid and Wigner phase-space distributions in the zero reciprocal temperature β limit. A Centroid Molecular Dynamics (CMD) approximation to the state-projected exact quantum dynamics is proposed and proven to be exact in the harmonic limit. The state projected CMD method is also tested numerically for a quartic oscillator and a double-well potential and found to be more accurate than canonical CMD. In the case of a ground state projection, this method can resolve tunnelling splittings of the double well problem in the higher barrier regime where canonical CMD fails. Finally, the state-projected CMD framework is cast in a path integral form.

Quest for Inexpensive Hydrogen Isotopic Fractionation: Do We Need 2D Quantum Confining in Porous Materials or Are Rough Surfaces Enough? The Case of Ammonia Nanoclusters

We study the adsorption energetics and quantum properties of the molecular hydrogen isotopes H₂, D₂, and T₂ onto the surface of rigid ammonia nanoclusters with quantum simulations and accurate model potential energy surfaces (PES). A highly efficient diffusion Monte Carlo (DMC) algorithm for rigid rotors allowed us to accurately define zero-point adsorption energies for the three isotopes, as well as the degree of translational and rotational delocalization that each affords on the surface. From the data emerges that the quantum adsorption energy (E_ads) of T₂ can be up to twice the one of H₂ at 0 K, suggesting the possibility of exploiting some form of solid ammonia to selectivity separate hydrogen isotopes at low temperatures (≃20 K). This is discussed by focusing on the structural motif that may be more effective for the task. The analysis of the contributions to E_ads, however, surprisingly indicates that the average kinetic energy (E^kin) and rotation energy (E_rot^kin) of T₂ can also be, respectively, 2 times and 20 times higher than those of H₂; this finding markedly deviates from what is predicted for hydrogen molecules inside carbon nanotubes (CNT) or metallic-organic frameworks (MOF), where E^kin and E_rot^kin is higher for H₂ due to the unavoidable effects of confinement and hindrance to its rotational motion. The rationale for these differences is provided by the geometrical distributions for the rigid rotors, which reveal an increasingly stronger coupling between rotational and translational degrees of freedom upon increasing the isotopic mass. This effect has never been observed before on adsorbing surfaces (e.g., graphite) and is induced by a strongly anisotropic and anharmonic bowl-like potential experienced by the rotors.