Energy benchmarks for methane-water systems from quantum Monte Carlo and second-order Møller-Plesset calculations.

The quantum Monte Carlo (QMC) technique is used to generate accurate energy benchmarks for methane-water clusters containing a single methane monomer and up to 20 water monomers. The benchmarks for each type of cluster are computed for a set of geometries drawn from molecular dynamics simulations. The accuracy of QMC is expected to be comparable with that of coupled-cluster calculations, and this is confirmed by comparisons for the CH4-H2O dimer. The benchmarks are used to assess the accuracy of the second-order Møller-Plesset (MP2) approximation close to the complete basis-set limit. A recently developed embedded many-body technique is shown to give an efficient procedure for computing basis-set converged MP2 energies for the large clusters. It is found that MP2 values for the methane binding energies and the cohesive energies of the water clusters without methane are in close agreement with the QMC benchmarks, but the agreement is aided by partial cancelation between 2-body and beyond-2-body errors of MP2. The embedding approach allows MP2 to be applied without loss of accuracy to the methane hydrate crystal, and it is shown that the resulting methane binding energy and the cohesive energy of the water lattice agree almost exactly with recently reported QMC values.

Sensory nerve fibers containing calcitonin gene-related peptide in gastrocnemius, latissimus dorsi and erector spinae muscles and thoracolumbar fascia in mice.

Chronic pain is a significant burden and much is attributed to back muscles. Back muscles and their associated fasciae make important and distinct contributions to back pain. Peptidergic nociceptors innervating these structures contribute to central transmission and pain modulation by peripheral and central actions. Plastic changes that augment and prolong pain are exhibited by neurons containing calcitonin gene-related peptide (CGRP) following muscle injury. Subpopulations of neurons containing this peptide have been identified in dorsal root ganglia but the distribution of their fibers in skeletal muscles and associated fasciae has not been fully documented. This study used multiple-labeling immunofluorescence and retrograde axonal tracing to identify dorsal root ganglion cells associated with muscle, and to characterize the distribution and density of their nerve fibers in mouse gastrocnemius and back muscles and in the thoracolumbar fascia. Most nerve fibers in these tissues contained CGRP and two major subpopulations of neurons were found: those containing CGRP and substance P (SP) and those containing CGRP but not SP. Innervation density was three times higher in the thoracolumbar fascia than in muscles of the back. These studies show mouse back and leg muscles are predominantly innervated by neurons containing CGRP, an important modulator of pain signal transmission. There are two distinct populations of neurons containing this peptide and their fibers were three times more densely distributed in the thoracolumbar fascia than back muscles.

Many-body exchange-overlap interactions in rare gases and water.

Generalized-gradient approximations (GGAs) of density-functional theory can suffer from substantial many-body errors in molecular systems interacting through weak non-covalent forces. Here, the errors of a range of GGAs for the 3-body energies of trimers of rare gases and water are investigated. The patterns of 3-body errors are similar for all the systems, and are related to the form of the exchange-enhancement factor FX(x) at large reduced gradient x, which also governs 2-body exchange-overlap errors. However, it is shown that the 3-body and 2-body errors depend in opposite ways on FX(x), so that they tend to cancel in molecular aggregates. Embedding arguments are used to achieve a partial separation of contributions to 3-body error from polarization, non-local correlation, and exchange, and it emerges that exchange is a major contributor. The practical importance of beyond-2-body errors is illustrated by the energetics of the water hexamer. An analysis of exchange-energy distributions is used to elucidate why 2-body and 3-body errors of GGAs depend in opposite ways on FX(x). The relevance of the present analysis to a range of other molecular systems is noted.

Analyzing the errors of DFT approximations for compressed water systems.

We report an extensive study of the errors of density functional theory (DFT) approximations for compressed water systems. The approximations studied are based on the widely used PBE and BLYP exchange-correlation functionals, and we characterize their errors before and after correction for 1- and 2-body errors, the corrections being performed using the methods of Gaussian approximation potentials. The errors of the uncorrected and corrected approximations are investigated for two related types of water system: first, the compressed liquid at temperature 420 K and density 1.245 g/cm(3) where the experimental pressure is 15 kilobars; second, thermal samples of compressed water clusters from the trimer to the 27-mer. For the liquid, we report four first-principles molecular dynamics simulations, two generated with the uncorrected PBE and BLYP approximations and a further two with their 1- and 2-body corrected counterparts. The errors of the simulations are characterized by comparing with experimental data for the pressure, with neutron-diffraction data for the three radial distribution functions, and with quantum Monte Carlo (QMC) benchmarks for the energies of sets of configurations of the liquid in periodic boundary conditions. The DFT errors of the configuration samples of compressed water clusters are computed using QMC benchmarks. We find that the 2-body and beyond-2-body errors in the liquid are closely related to similar errors exhibited by the clusters. For both the liquid and the clusters, beyond-2-body errors of DFT make a substantial contribution to the overall errors, so that correction for 1- and 2-body errors does not suffice to give a satisfactory description. For BLYP, a recent representation of 3-body energies due to Medders, Babin, and Paesani [J. Chem. Theory Comput. 9, 1103 (2013)] gives a reasonably good way of correcting for beyond-2-body errors, after which the remaining errors are typically 0.5 mE(h) ≃ 15 meV/monomer for the liquid and the clusters.

Energy benchmarks for water clusters and ice structures from an embedded many-body expansion.

We show how an embedded many-body expansion (EMBE) can be used to calculate accurate ab initio energies of water clusters and ice structures using wavefunction-based methods. We use the EMBE described recently by Bygrave et al. [J. Chem. Phys. 137, 164102 (2012)], in which the terms in the expansion are obtained from calculations on monomers, dimers, etc., acted on by an approximate representation of the embedding field due to all other molecules in the system, this field being a sum of Coulomb and exchange-repulsion fields. Our strategy is to separate the total energy of the system into Hartree-Fock and correlation parts, using the EMBE only for the correlation energy, with the Hartree-Fock energy calculated using standard molecular quantum chemistry for clusters and plane-wave methods for crystals. Our tests on a range of different water clusters up to the 16-mer show that for the second-order Møller-Plesset (MP2) method the EMBE truncated at 2-body level reproduces to better than 0.1 mE(h)/monomer the correlation energy from standard methods. The use of EMBE for computing coupled-cluster energies of clusters is also discussed. For the ice structures Ih, II, and VIII, we find that MP2 energies near the complete basis-set limit reproduce very well the experimental values of the absolute and relative binding energies, but that the use of coupled-cluster methods for many-body correlation (non-additive dispersion) is essential for a full description. Possible future applications of the EMBE approach are suggested.

Communication: energy benchmarking with quantum Monte Carlo for water nano-droplets and bulk liquid water.

We show the feasibility of using quantum Monte Carlo (QMC) to compute benchmark energies for configuration samples of thermal-equilibrium water clusters and the bulk liquid containing up to 64 molecules. Evidence that the accuracy of these benchmarks approaches that of basis-set converged coupled-cluster calculations is noted. We illustrate the usefulness of the benchmarks by using them to analyze the errors of the popular BLYP approximation of density functional theory (DFT). The results indicate the possibility of using QMC as a routine tool for analyzing DFT errors for non-covalent bonding in many types of condensed-phase molecular system.

First-principles energetics of water clusters and ice: a many-body analysis.

Standard forms of density-functional theory (DFT) have good predictive power for many materials, but are not yet fully satisfactory for cluster, solid, and liquid forms of water. Recent work has stressed the importance of DFT errors in describing dispersion, but we note that errors in other parts of the energy may also contribute. We obtain information about the nature of DFT errors by using a many-body separation of the total energy into its 1-body, 2-body, and beyond-2-body components to analyze the deficiencies of the popular PBE and BLYP approximations for the energetics of water clusters and ice structures. The errors of these approximations are computed by using accurate benchmark energies from the coupled-cluster technique of molecular quantum chemistry and from quantum Monte Carlo calculations. The systems studied are isomers of the water hexamer cluster, the crystal structures Ih, II, XV, and VIII of ice, and two clusters extracted from ice VIII. For the binding energies of these systems, we use the machine-learning technique of Gaussian Approximation Potentials to correct successively for 1-body and 2-body errors of the DFT approximations. We find that even after correction for these errors, substantial beyond-2-body errors remain. The characteristics of the 2-body and beyond-2-body errors of PBE are completely different from those of BLYP, but the errors of both approximations disfavor the close approach of non-hydrogen-bonded monomers. We note the possible relevance of our findings to the understanding of liquid water.

Assessing the accuracy of quantum Monte Carlo and density functional theory for energetics of small water clusters.

We present a detailed study of the energetics of water clusters (H(2)O)(n) with n ≤ 6, comparing diffusion Monte Carlo (DMC) and approximate density functional theory (DFT) with well converged coupled-cluster benchmarks. We use the many-body decomposition of the total energy to classify the errors of DMC and DFT into 1-body, 2-body and beyond-2-body components. Using both equilibrium cluster configurations and thermal ensembles of configurations, we find DMC to be uniformly much more accurate than DFT, partly because some of the approximate functionals give poor 1-body distortion energies. Even when these are corrected, DFT remains considerably less accurate than DMC. When both 1- and 2-body errors of DFT are corrected, some functionals compete in accuracy with DMC; however, other functionals remain worse, showing that they suffer from significant beyond-2-body errors. Combining the evidence presented here with the recently demonstrated high accuracy of DMC for ice structures, we suggest how DMC can now be used to provide benchmarks for larger clusters and for bulk liquid water.

The kinetics of homogeneous melting beyond the limit of superheating.

Molecular dynamics simulation is used to study the time-scales involved in the homogeneous melting of a superheated crystal. The interaction model used is an embedded-atom model for Fe developed in previous work, and the melting process is simulated in the microcanonical (N, V, E) ensemble. We study periodically repeated systems containing from 96 to 7776 atoms, and the initial system is always the perfect crystal without free surfaces or other defects. For each chosen total energy E and number of atoms N, we perform several hundred statistically independent simulations, with each simulation lasting for between 500 ps and 10 ns, in order to gather statistics for the waiting time τ(w) before melting occurs. We find that the probability distribution of τ(w) is roughly exponential, and that the mean value <τ(w)> depends strongly on the excess of the initial steady temperature of the crystal above the superheating limit identified by other researchers. The mean <τ(w)> also depends strongly on system size in a way that we have quantified. For very small systems of ~100 atoms, we observe a persistent alternation between the solid and liquid states, and we explain why this happens. Our results allow us to draw conclusions about the reliability of the recently proposed Z method for determining the melting properties of simulated materials and to suggest ways of correcting for the errors of the method.

Ab initio statistical mechanics of surface adsorption and desorption. II. Nuclear quantum effects.

We show how the path-integral formulation of quantum statistical mechanics can be used to construct practical ab initio techniques for computing the chemical potential of molecules adsorbed on surfaces, with full inclusion of quantum nuclear effects. The techniques we describe are based on the computation of the potential of mean force on a chosen molecule and generalize the techniques developed recently for classical nuclei. We present practical calculations based on density functional theory with a generalized-gradient exchange-correlation functional for the case of H(2)O on the MgO (001) surface at low coverage. We note that the very high vibrational frequencies of the H(2)O molecule would normally require very large numbers of time slices (beads) in path-integral calculations, but we show that this requirement can be dramatically reduced by employing the idea of thermodynamic integration with respect to the number of beads. The validity and correctness of our path-integral calculations on the H(2)O/MgO(001) system are demonstrated by supporting calculations on a set of simple model systems for which quantum contributions to the free energy are known exactly from analytic arguments.

Spatially local parallel tempering for thermal-equilibrium sampling.

Parallel tempering (PT) is a set of techniques for accelerating thermal-equilibrium sampling in systems where the exploration of configuration space is hindered by energy barriers. With standard PT algorithms, the computational effort scales unfavorably with system size, so that it is difficult to apply them to large systems. We propose local PT algorithms, for which the computational effort is proportional to the number of degrees of freedom. We demonstrate the effectiveness of the new algorithms on two one-dimensional model systems, showing that results for selected observables are correctly reproduced, and that practical linear scaling is achieved. We show also that the algorithms are readily applied to systems in higher dimensions. We note the prospects for studying large extended systems, including surfaces and interfaces.

Melting properties of a simple tight-binding model of transition metals. I. The region of half-filled d-band.

We present calculations of the free energy, and hence the melting properties, of a simple tight-binding model for transition metals in the region of d-band filling near the middle of a d-series, the parameters of the model being designed to mimic molybdenum. The melting properties are calculated for pressures ranging from ambient to several megabars. The model is intended to be the simplest possible tight-binding representation of the two basic parts of the energy: first, the pairwise repulsion due to Fermi exclusion; and second, the d-band bonding energy described in terms of an electronic density of states that depends on structure. In addition to the number of d-electrons, the model contains four parameters, which are adjusted to fit the pressure dependent d-band-width and the zero-temperature pressure-volume relation of Mo. We show that the resulting model reproduces well the phonon dispersion relations of Mo in the body-centered-cubic structure, as well as the radial distribution function of the high-temperature solid and liquid given by earlier first-principles simulations. Our free energy calculations start from the free energy of the liquid and solid phases of the purely repulsive pair potential model, without d-band bonding. The free energy of the full tight-binding model is obtained from this by thermodynamic integration. The resulting melting properties of the model are quite close to those given by earlier first-principles work on Mo. An interpretation of these melting properties is provided by showing how they are related to those of the purely repulsive model.

High-precision calculation of Hartree-Fock energy of crystals.

When using quantum chemistry techniques to calculate the energetics of bulk crystals, there is a need to calculate the Hartree-Fock (HF) energy of the crystal at the basis-set limit. We describe a strategy for achieving this, which exploits the fact that the HF energy of crystals can now be calculated using pseudopotentials and plane-wave basis sets, an approach that permits basis-set convergence to arbitrary precision. The errors due to the use of pseudopotentials are then computed from the difference of all-electron and pseudopotential total energies of atomic clusters, extrapolated to the bulk-crystal limit. The strategy is tested for the case of the LiH crystal, and it is shown that the HF cohesive energy can be converged with respect to all technical parameters to a precision approaching 0.1 mE(h) per atom. This cohesive energy and the resulting HF value of the equilibrium lattice parameter are compared with literature values obtained using Gaussian basis sets.

Atomic structure of misfit dislocations at InAs/GaAs(110).

Heteroepitaxy of InAs on GaAs(110) leads to the formation of subsurface misfit dislocations to relieve strain. These dislocations have been observed both with transmission electron microscopy (TEM) and scanning tunnelling microscopy (STM), and show regular spacing. Electronic structure calculations of the structure of the core of the dislocations, as well as their location within the epitaxial layer, are presented. The most stable location is found to be at the interface, with the core centred over In. Calculated strain profiles and the thickness at which dislocations should form are in good agreement with available experimental data.

Ab initio statistical mechanics of surface adsorption and desorption. I. H2O on MgO (001) at low coverage.

We present a general computational scheme based on molecular dynamics (MD) simulation for calculating the chemical potential of adsorbed molecules in thermal equilibrium on the surface of a material. The scheme is based on the calculation of the mean force in MD simulations in which the height of a chosen molecule above the surface is constrained and subsequent integration of the mean force to obtain the potential of mean force and hence the chemical potential. The scheme is valid at any coverage and temperature, so that in principle it allows the calculation of the chemical potential as a function of coverage and temperature. It avoids all statistical mechanical approximations, except for the use of classical statistical mechanics for the nuclei, and assumes nothing in advance about the adsorption sites. From the chemical potential, the absolute desorption rate of the molecules can be computed, provided that the equilibration rate on the surface is faster than the desorption rate. We apply the theory by ab initio MD simulation to the case of H2O on MgO (001) in the low-coverage limit, using the Perdew-Burke-Ernzerhof (PBE) form of exchange correlation. The calculations yield an ab initio value of the Polanyi-Wigner frequency prefactor, which is more than two orders of magnitude greater than the value of 10(13) s(-1) often assumed in the past. Provisional comparison with experiment suggests that the PBE adsorption energy may be too low, but the extension of the calculations to higher coverages is needed before firm conclusions can be drawn. The possibility of including quantum nuclear effects by using path-integral simulations is noted.

Ab initio melting curve of molybdenum by the phase coexistence method.

Ab initio calculations of the melting curve of molybdenum for the pressure range 0-400 GPa are reported. The calculations employ density functional theory (DFT) with the Perdew-Burke-Ernzerhof exchange-correlation functional in the projector augmented wave (PAW) implementation. Tests are presented showing that these techniques accurately reproduce experimental data on low-temperature body-centered cubic (bcc) Mo, and that PAW agrees closely with results from the full-potential linearized augmented plane-wave implementation. The work attempts to overcome the uncertainties inherent in earlier DFT calculations of the melting curve of Mo, by using the "reference coexistence" technique to determine the melting curve. In this technique, an empirical reference model (here, the embedded-atom model) is accurately fitted to DFT molecular dynamics data on the liquid and the high-temperature solid, the melting curve of the reference model is determined by simulations of coexisting solid and liquid, and the ab initio melting curve is obtained by applying free-energy corrections. The calculated melting curve agrees well with experiment at ambient pressure and is consistent with shock data at high pressure, but does not agree with the high-pressure melting curve deduced from static compression experiments. Calculated results for the radial distribution function show that the short-range atomic order of the liquid is very similar to that of the high-T solid, with a slight decrease of coordination number on passing from solid to liquid. The electronic densities of states in the two phases show only small differences. The results do not support a recent theory according to which very low dT(m)dP values are expected for bcc transition metals because of electron redistribution between s-p and d states.

Density functional calculations of surface free energies.

We propose a general method of thermodynamic integration to find the free energy of a surface, where our integration parameter is taken to be the strain on the unit cell of the system (which in the example presented in this paper is simply the extension of the unit cell along the normal to the surface), and the integration is performed over the thermal average stress from a molecular dynamics run. In order to open up a vacuum gap in a continuous and reversible manner, an additional control interaction has been introduced. We also use temperature integration to find a linear relation for the temperature dependence of the free surface energy. These methods have been applied to the titanium dioxide (110) surface, using first principles density functional theory. A proof of principle calculation for zero temperature shows excellent agreement between the integral calculation and the difference in energy calculated by the DFT program. Calculations that have been performed at 295 and 1000 K give excellent agreement between the two integration methods.

Extension of molecular electronic structure methods to the solid state: computation of the cohesive energy of lithium hydride.

We describe a simple strategy for calculating the cohesive energy of certain kinds of crystal using readily available quantum chemistry techniques. The strategy involves the calculation of the electron correlation energies of a hierarchy of free clusters, and the cohesive energy E(coh) is extracted from the constant of proportionality between these correlation energies and the number of atoms in the limit of large clusters. We apply the strategy to the LiH crystal, using the MP2 and CCSD(T) schemes for the correlation energy, and show that for this material E(coh) can be obtained to an accuracy of approximately 30 meV per ion pair. Comparison with the experimental value, after correction for zero-point energy, confirms this accuracy.

Absolute rate of thermal desorption from first-principles simulation.

We present a technique for computing by first-principles simulation the absolute desorption rate γ of adsorbate molecules from a surface for any coverage and temperature. The technique is valid when the thermal equilibration rate on the surface is faster than γ, and is based on an exact expression for γ in terms of the difference of non-configurational chemical potentials of gas-phase and adsorbed molecules. This difference is expressed in terms of a potential of mean force, which is computed by constrained first-principles molecular dynamics. The technique is applied to D(2)O on the MgO(001) surface at low coverage, using the generalized gradient approximation (GGA) for exchange-correlation energy. Comparisons with experimental temperature programmed desorption data allow an assessment of the accuracy of the GGA for the adsorption of D(2)O on MgO(001).