Implementation of the self-consistent phonons method with ab initio potentials (AI-SCP).

The self-consistent phonon (SCP) method allows one to include anharmonic effects when treating a many-body quantum system at thermal equilibrium. The system is then described by an effective temperature-dependent harmonic Hamiltonian, which can be used to estimate its various dynamic and static properties. In this paper, we combine SCP with ab initio (AI) potential energy evaluation in which case the numerical bottleneck of AI-SCP is the evaluation of Gaussian averages of the AI potential energy and its derivatives. These averages are computed efficiently by the quasi-Monte Carlo method utilizing low-discrepancy sequences leading to a fast convergence with respect to the number, S, of the AI energy evaluations. Moreover, a further substantial (an-order-of-magnitude) improvement in efficiency is achieved once a numerically cheap approximation of the AI potential is available. This is based on using a perturbation theory-like (the two-grid) approach in which it is the average of the difference between the AI and the approximate potential that is computed. The corresponding codes and scripts are provided.

Molecular Spectra Calculations Using an Optimized Quasi-Regular Gaussian Basis and the Collocation Method.

We revisit the collocation method of Manzhos and Carrington [ J. Chem. Phys., 2016, 145, 224110] in which a distributed localized (e.g., Gaussian) basis is used to set up a generalized eigenvalue problem to compute the eigenenergies and eigenfunctions of a molecular vibrational Hamiltonian. Although the resulting linear algebra problem involves full matrices, the method provides a number of important advantages, namely, (i) it is very simple both conceptually and numerically, (ii) it can be formulated using any set of internal molecular coordinates, (iii) it is flexible with respect to the choice of the basis, (iv) no integrals need to be computed, and (v) it has the potential to significantly reduce the basis size through optimizing the placement and the shapes of the basis functions. In the present paper, we explore the latter aspect of the method using the recently introduced, and here further improved, quasi-regular grids (QRGs). By computing the eigenenergies of the four-atom molecule of formaldehyde, we demonstrate that a QRG-based distributed Gaussian basis is superior to the previously used choices.

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

para-Hydrogen (pH2)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 (pH2)N. However, their structural analysis was typically limited to attempts to identify "magic number clusters" by computing their ground state energies EN and the chemical potential µN = EN-EN-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.

Magic numbers, quantum delocalization, and orientational disordering in anionic hydrogen and deuterium clusters.

The Diffusion Monte Carlo (DMC) method was applied to anionic hydrogen clusters H-(H2)n (n = 1-16, 32) and their deuterated analogs using a polarizable all-atom potential energy surface (PES) developed by Calvo and Yurtsever. For the hydrogen clusters, the binding energy ΔEn appears to be a smooth function of the cluster size n, thus contradicting the previous claim that n = 12 is a "magic number" cluster. The structures of the low energy minima of the PES for these clusters belong to the icosahedral motif with the H2 molecules aligned toward the central H- ion. However, their ground state wavefunctions are highly delocalized and resemble neither the structures of the global nor local minima. Moreover, the strong nuclear quantum effects result in a nearly complete orientational disordering of the H2 molecules. For the deuterium clusters, the ground state wavefunctions are localized and the D2 molecules are aligned toward the central D- ion. However, their structures are still characterized as disordered and, as such, do not display size sensitivity. In addition, DMC simulations were performed on the mixed H-(H2)n(D2)p clusters with (n, p) = (6, 6) and (16, 16). Again, in contradiction to the previous claim, we found that the "more quantum" H2 molecules prefer to reside farther from the central H- ion than the D2 molecules.

Sampling general distributions with quasi-regular grids: Application to the vibrational spectra calculations.

We introduce a new method for sampling a general multidimensional distribution function Px using a quasiregular grid (QRG) of points xi (i = 1, …, N). This grid is constructed by minimizing a pairwise functional, ∑u(xi, xj) → min, with the short-range pair pseudopotential u(xi, xj), defined locally according to the underlying distribution P(x). While QRGs can be useful in many diverse areas of science, in this paper, we apply them to construct Gaussian basis sets in the context of solving the vibrational Schrödinger equation. Using some 2D and 3D model systems, we demonstrate that the resulting optimized Gaussian basis sets have properties superior to other choices explored previously in the literature.

Quantum-induced solid-solid transitions and melting in the Lennard-Jones LJ38 cluster.

The solid-solid and melting transitions that occur in Lennard-Jones LJ n clusters have been both fascinating and challenging for the computational physics community over the last several decades. A number of attempts to extend these studies to the quantum case have also been made. Particularly interesting is the exploration of the parallel between the thermally induced and quantum-induced transitions. Yet, both numerically accurate and systematic studies of the latter are still lacking. In this paper, we apply the diffusion Monte Carlo method to the especially difficult case of LJ38. Starting with the truncated octahedral global minimum configuration, as the de Boer quantum delocalization parameter Λ increases, the system undergoes two consecutive solid-solid transitions, switching to anti-Mackay configurations. At sufficiently large values of Λ, the cluster is completely "melted," which is manifested by delocalization of the ground state wavefunction over a very large number of minima that represent several structural motifs.

Nuclear Quantum Effects and Thermodynamic Properties for Small (H2O)1-21X- Clusters (X- = F-, Cl-, Br-, I-).

We carried out accurate diffusion Monte Carlo (DMC) studies for small (H2O) NX- clusters ( N = 1-5; X- = F-, Cl-, Br-, I-) and their D2O isotopologues. We found remarkably good agreement (i.e., ∼0.1 kcal/mol or better) with no exceptions between the DMC solvation energies and the corresponding harmonic approximation (HA) estimates, due, apparently, to massive error cancellations. This is surprising, in particular, because HA does not account for a substantial (i.e., ∼ 3%) increase of the mean O-O distances, caused by the anharmonicity in conjunction with the nuclear quantum effects, although the other distances in the system are affected to a much lesser extent. This agreement for the solvation energies motivated us to extend the current study to larger ( N = 6-21) clusters to explore their thermodynamic properties using the harmonic superposition method (HSM). The HSM results for the solvation free energies in turn reveal that at finite temperatures the nuclear quantum effects (including the isotope effects) in these systems are miniscule.

Quantum Melting and Isotope Effects from Diffusion Monte Carlo Studies of p-H2 Clusters.

We present a rigorous characterization of the ground state structures of p-H2 clusters and their isotopologues using diffusion Monte Carlo combined with the inherent structures analysis. For the N = 19 cluster we explore the effect of "quantum melting" by quantifying the contributions of local minima to the ground state as a function of continuously varying particle mass. Doubling the cluster size leads to an enormous increase of its complexity: the ground state of (p-H2)38 is highly delocalized over a large number of minima representing all the funnels of the potential energy surface. The ground state of (o-D2)38 is also delocalized, but over a smaller subset of minima, which exclusively belong to the same disordered motif.

Monitoring Water Clusters "Melt" Through Vibrational Spectroscopy.

Characterizing structural and phase transformations of water at the molecular level is key to understanding a variety of multiphase processes ranging from ice nucleation in the atmosphere to hydration of biomolecules and wetting of solid surfaces. In this study, state-of-the-art quantum simulations with a many-body water potential energy surface, which exhibits chemical and spectroscopic accuracy, are carried out to monitor the microscopic melting of the water hexamer through the analysis of vibrational spectra and appropriate structural order parameters as a function of temperature. The water hexamer is specifically chosen as a case study due to the central role of this cluster in the molecular-level understanding of hydrogen bonding in water. Besides being in agreement with the experimental data available for selected isomers at very low temperature, the present results provide quantitative insights into the interplay between energetic, entropic, and nuclear quantum effects on the evolution of water clusters from "solid-like" to "liquid-like" structures. This study thus demonstrates that computer simulations can now bridge the gap between measurements currently possible for individual isomers at very low temperature and observations of isomer mixtures at ambient conditions.

Binding energies from diffusion Monte Carlo for the MB-pol H2O and D2O dimer: A comparison to experimental values.

The diffusion Monte Carlo (DMC) method is applied to compute the ground state energies of the water monomer and dimer and their D2O isotopomers using MB-pol; the most recent and most accurate ab inito-based potential energy surface (PES). MB-pol has already demonstrated excellent agreement with high level electronic structure data, as well as agreement with some experimental, spectroscopic, and thermodynamic data. Here, the DMC binding energies of (H2O)2 and (D2O)2 agree with the corresponding values obtained from velocity map imaging within, respectively, 0.01 and 0.02 kcal/mol. This work adds two more valuable data points that highlight the accuracy of the MB-pol PES.

Assessing the Performance of the Diffusion Monte Carlo Method As Applied to the Water Monomer, Dimer, and Hexamer.

The diffusion Monte Carlo (DMC) method is applied to the water monomer, dimer, and hexamer using q-TIP4P/F, one of the most simple empirical water models with flexible monomers. The bias in the time step (Δτ) and population size (Nw) is investigated. For the binding energies, the bias in Δτ cancels nearly completely, whereas a noticeable bias in Nw remains. However, for the isotope shift (e.g, in the dimer binding energies between (H2O)2 and (D2O)2), the systematic errors in Nw do cancel. Consequently, very accurate results for the latter (within ∼0.01 kcal/mol) are obtained with moderate numerical effort (Nw ∼ 10(3)). For the water hexamer and its (D2O)6 isotopomer, the DMC results as a function of Nw are examined for the cage and prism isomers. For a given isomer, the issue of the walker population leaking out of the corresponding basin of attraction is addressed by using appropriate geometric constraints. The population size bias for the hexamer is more severe, and to maintain accuracy similar to that of the dimer, Nw must be increased by ∼2 orders of magnitude. Fortunately, when the energy difference between the cage and prism is taken, the biases cancel, thereby reducing the systematic errors to within ∼0.01 kcal/mol when using a population of Nw = 4.8 × 10(5) walkers. Consequently, a very accurate result for the isotope shift is also obtained. Notably, both the quantum and isotope effects for the prism-cage energy difference are small.

On the ground state calculation of a many-body system using a self-consistent basis and quasi-Monte Carlo: an application to water hexamer.

Given a quantum many-body system, the Self-Consistent Phonons (SCP) method provides an optimal harmonic approximation by minimizing the free energy. In particular, the SCP estimate for the vibrational ground state (zero temperature) appears to be surprisingly accurate. We explore the possibility of going beyond the SCP approximation by considering the system Hamiltonian evaluated in the harmonic eigenbasis of the SCP Hamiltonian. It appears that the SCP ground state is already uncoupled to all singly- and doubly-excited basis functions. So, in order to improve the SCP result at least triply-excited states must be included, which then reduces the error in the ground state estimate substantially. For a multidimensional system two numerical challenges arise, namely, evaluation of the potential energy matrix elements in the harmonic basis, and handling and diagonalizing the resulting Hamiltonian matrix, whose size grows rapidly with the dimensionality of the system. Using the example of water hexamer we demonstrate that such calculation is feasible, i.e., constructing and diagonalizing the Hamiltonian matrix in a triply-excited SCP basis, without any additional assumptions or approximations. Our results indicate particularly that the ground state energy differences between different isomers (e.g., cage and prism) of water hexamer are already quite accurate within the SCP approximation.

Filter diagonalization method for processing PFG NMR data.

Obtaining diffusion coefficients from PFG NMR diffusion (a.k.a DOSY) data is, in the general case, an ill-posed problem. Numerous methods for processing such data have therefore been developed, each with different constraints and assumptions. The Regularized Resolvent Transform (RRT) is a proven and robust method for spectral inversion. In earlier papers RRT, albeit very slow, was argued to be superior for DOSY processing to a related algorithm, the Filter Diagonalization Method (FDM). Here FDM is revisited and a new regularization method is implemented, which drastically improves the performance and provides spectra of comparable or better quality to those provided by RRT. Both the RRT and the FDM for DOSY processing have been implemented as options in the free and open source DOSY Toolbox.

Mapping the phase diagram for neon to a quantum Lennard-Jones fluid using Gibbs ensemble simulations.

In order to address the issue of whether neon liquid in coexistence with its gas phase can be mapped to a quantum Lennard-Jones (LJ) fluid, we perform a series of simulations using Gibbs ensemble Monte Carlo for a range of de Boer quantum parameters Λ=ℏ/(σ√(mε)). The quantum effects are incorporated by implementing the variational gaussian wavepacket method, which provides an efficient numerical framework for estimating the quantum density at thermal equilibrium. The computed data for the LJ liquid is used to produce its phase diagram as a function of the quantum parameter, 0.065 ≤ Λ ≤ 0.11. These data are then used to fit the experimental phase diagram for neon liquid. The resulting parameters, ε = 35.68 ± 0.03 K and σ = 2.7616 ± 0.0005 Å (Λ = 0.0940), of the LJ pair potential are optimized to best represent liquid neon in coexistence with its gas phase for a range of physically relevant temperatures. This multi-temperature approach towards fitting and assessing a pair-potential is much more consistent than merely fitting a single data point, such as a melting temperature or a second virial coefficient.

Self-consistent phonons revisited. II. A general and efficient method for computing free energies and vibrational spectra of molecules and clusters.

The self-consistent phonons (SCP) method provides a consistent way to include anharmonic effects when treating a many-body quantum system at thermal equilibrium. The system is then described by an effective temperature-dependent harmonic Hamiltonian, which can be used to estimate the system's properties, such as its free energy or its vibrational spectrum. The numerical bottleneck of the method is the evaluation of Gaussian averages of the potential energy and its derivatives. Several algorithmic ideas/tricks are introduced to reduce the cost of such integration by orders of magnitude, e.g., relative to that of the previous implementation of the SCP approach by Calvo et al. [J. Chem. Phys. 133, 074303 (2010)]. One such algorithmic improvement is the replacement of standard Monte Carlo integration by quasi-Monte Carlo integration utilizing low-discrepancy sequences. The performance of the method is demonstrated on the calculation of vibrational frequencies of pyrene. It is then applied to compute the free energies of five isomers of water hexamer using the WHBB potential of Bowman and co-workers [J. Chem. Phys. 134, 094509 (2011)]. The present results predict the hexamer prism being thermodynamically most stable, with the free energy of the hexamer cage being about 0.2 kcal mol(-1) higher at all temperatures below T = 200 K.

Self-consistent phonons revisited. I. The role of thermal versus quantum fluctuations on structural transitions in large Lennard-Jones clusters.

The theory of self-consistent phonons (SCP) was originally developed to address the anharmonic effects in condensed matter systems. The method seeks a harmonic, temperature-dependent Hamiltonian that provides the "best fit" for the physical Hamiltonian, the "best fit" being defined as the one that optimizes the Helmholtz free energy at a fixed temperature. The present developments provide a scalable O(N) unified framework that accounts for anharmonic effects in a many-body system, when it is probed by either thermal (ℏ → 0) or quantum fluctuations (T → 0). In these important limits, the solution of the nonlinear SCP equations can be reached in a manner that requires only the multiplication of 3N × 3N matrices, with no need of diagonalization. For short range potentials, such as Lennard-Jones, the Hessian, and other related matrices are highly sparse, so that the scaling of the matrix multiplications can be reduced from O(N(3)) to ~O(N). We investigate the role of quantum effects by continuously varying the de-Boer quantum delocalization parameter Λ and report the N-Λ (T = 0), and also the classical N-T (Λ = 0) phase diagrams for sizes up to N ~ 10(4). Our results demonstrate that the harmonic approximation becomes inadequate already for such weakly quantum systems as neon clusters, or for classical systems much below the melting temperatures.

A fast variational Gaussian wavepacket method: size-induced structural transitions in large neon clusters.

The variational Gaussian wavepacket (VGW) approximation provides an alternative to path integral Monte Carlo for the computation of thermodynamic properties of many-body systems at thermal equilibrium. It provides a direct access to the thermal density matrix and is particularly efficient for Monte Carlo approaches, as for an N-body system it operates in a non-inflated 3N-dimensional configuration space. Here, we greatly accelerate the VGW method by retaining only the relevant short-range correlations in the (otherwise full) 3N × 3N Gaussian width matrix without sacrificing the accuracy of the fully coupled VGW method. This results in the reduction of the original O(N(3)) scaling to O(N(2)). The fast-VGW method is then applied to quantum Lennard-Jones clusters with sizes up to N = 6500 atoms. Following Doye and Calvo [JCP 116, 8307 (2002)] we study the competition between the icosahedral and decahedral structural motifs in Ne(N) clusters as a function of N.

Thermal Gaussian molecular dynamics for quantum dynamics simulations of many-body systems: application to liquid para-hydrogen.

A new method, here called thermal Gaussian molecular dynamics (TGMD), for simulating the dynamics of quantum many-body systems has recently been introduced [I. Georgescu and V. A. Mandelshtam, Phys. Rev. B 82, 094305 (2010)]. As in the centroid molecular dynamics (CMD), in TGMD the N-body quantum system is mapped to an N-body classical system. The associated both effective Hamiltonian and effective force are computed within the variational Gaussian wave-packet approximation. The TGMD is exact for the high-temperature limit, accurate for short times, and preserves the quantum canonical distribution. For a harmonic potential and any form of operator Â, it provides exact time correlation functions C(AB)(t) at least for the case of B, a linear combination of the position, x, and momentum, p, operators. While conceptually similar to CMD and other quantum molecular dynamics approaches, the great advantage of TGMD is its computational efficiency. We introduce the many-body implementation and demonstrate it on the benchmark problem of calculating the velocity time auto-correlation function for liquid para-hydrogen, using a system of up to N = 2592 particles.

Single-molecule imaging of platinum ligand exchange reaction reveals reactivity distribution.

Single-molecule fluorescence microscopy provided information about the real-time distribution of chemical reactivity on silicon oxide supports at the solution-surface interface, at a level of detail which would be unavailable from a traditional ensemble technique or from a technique that imaged the static physical properties of the surface. Chemical reactions on the surface were found to be uncorrelated; that is, the chemical reaction of one metal complex did not influence the location of a future chemical reaction of another metal complex.