Showing papers on "Path integral molecular dynamics published in 2013"

Efficient methods and practical guidelines for simulating isotope effects

Abstract: The shift in chemical equilibria due to isotope substitution is frequently exploited to obtain insight into a wide variety of chemical and physical processes. It is a purely quantum mechanical effect, which can be computed exactly using simulations based on the path integral formalism. Here we discuss how these techniques can be made dramatically more efficient, and how they ultimately outperform quasi-harmonic approximations to treat quantum liquids not only in terms of accuracy, but also in terms of computational cost. To achieve this goal we introduce path integral quantum mechanics estimators based on free energy perturbation, which enable the evaluation of isotope effects using only a single path integral molecular dynamics trajectory of the naturally abundant isotope. We use as an example the calculation of the free energy change associated with H/D and 16O/18O substitutions in liquid water, and of the fractionation of those isotopes between the liquid and the vapor phase. In doing so, we demonstrate and discuss quantitatively the relative benefits of each approach, thereby providing a set of guidelines that should facilitate the choice of the most appropriate method in different, commonly encountered scenarios. The efficiency of the estimators we introduce and the analysis that we perform should in particular facilitate accurate ab initio calculation of isotope effects in condensed phase systems.

Quantum Drude oscillator model of atoms and molecules: many-body polarization and dispersion interactions for atomistic simulation.

Abstract: Treating both many-body polarization and dispersion interactions is now recognized as a key element in achieving the level of atomistic modeling required to reveal novel physics in complex systems. The quantum Drude oscillator (QDO), a Gaussian-based, coarse grained electronic structure model, captures both many-body polarization and dispersion and has linear scale computational complexity with system size, hence it is a leading candidate next-generation simulation method. Here, we investigate the extent to which the QDO treatment reproduces the desired long-range atomic and molecular properties. We present closed form expressions for leading order polarizabilities and dispersion coefficients and derive invariant (parameter-free) scaling relationships among multipole polarizability and many-body dispersion coefficients that arise due to the Gaussian nature of the model. We show that these ``combining rules'' hold to within a few percent for noble gas atoms, alkali metals, and simple (first-row hydride) molecules such as water; this is consistent with the surprising success that models with underlying Gaussian statistics often exhibit in physics. We present a diagrammatic Jastrow-type perturbation theory tailored to the QDO model that serves to illustrate the rich types of responses that the QDO approach engenders. QDO models for neon, argon, krypton, and xenon, designed to reproduce gas phase properties, are constructed and their condensed phase properties explored via linear scale diffusion Monte Carlo (DMC) and path integral molecular dynamics (PIMD) simulations. Good agreement with experimental data for structure, cohesive energy, and bulk modulus is found, demonstrating a degree of transferability that cannot be achieved using current empirical models or fully ab initio descriptions.

A Surface-Specific Isotope Effect in Mixtures of Light and Heavy Water

Abstract: Isotope fractionation between different phases is a subtle but very important phenomenon that is related to the quantum nature of light nuclei, and that has important consequences for geochemistry, hydrology, and biology. Here we present a joint experimental/theoretical investigation of the differential segregation of hydrogen and deuterium at the liquid/vapor interface of mixtures of light and heavy water. We use both vibrational sum-frequency spectroscopy and path integral molecular dynamics simulations to quantitatively assess this phenomenon. The experimental and theoretical results indicate that the last layer of water molecules at the liquid/vapor interface is enriched in hydrogen. We discuss in detail the extent, the physical origin, and the implications of this surface-specific isotope effect.

Identifying and characterising the different structural length scales in liquids and glasses: an experimental approach

Abstract: The structure of several network-forming liquids and glasses is considered, where a focus is placed on the detailed information that is made available by using the method of neutron diffraction with isotope substitution (NDIS). In the case of binary network glass-forming materials with the MX2 stoichiometry (e.g. GeO2, GeSe2, ZnCl2), two different length scales at distances greater than the nearest-neighbour distance manifest themselves by peaks in the measured diffraction patterns. The network properties are influenced by a competition between the ordering on these “intermediate” and “extended” length scales, which can be manipulated by changing the chemical identity of the atomic constituents or by varying state parameters such as the temperature and pressure. The extended-range ordering, which describes the decay of the pair-correlation functions at large-r, can be represented by making a pole analysis of the Ornstein–Zernike equations, an approach that can also be used to describe the large-r behaviour of the pair-correlation functions for liquid and amorphous metals where packing constraints are important. The first applications are then described of the NDIS method to measure the detailed structure of aerodynamically-levitated laser-heated droplets of “fragile” glass-forming liquid oxides (CaAl2O4 and CaSiO3) at high-temperatures (∼2000 K) and the structure of a “strong” network-forming glass (GeO2) under pressures ranging from ambient to ∼8 GPa. The high-temperature experiments show structural changes on multiple length scales when the oxides are vitrified. The high-pressure experiment offers insight into the density-driven mechanisms of network collapse in GeO2 glass, and parallels are drawn with the high-pressure behaviour of silica glass. Finally, the hydrogen-bonded network of water is considered, where the first application of the method of oxygen NDIS is used to measure the structures of light versus heavy water and a difference of ≃0.5% is found between the O–D and O–H intra-molecular bond lengths. The experimental data are best matched by using path integral molecular dynamics simulations with a flexible anharmonic water model, and the results support a competing quantum effects model for water in which its structural and dynamical properties are governed by an offset between intra-molecular and inter-molecular quantum contributions.

Enabling Forbidden Processes: Quantum and Solvation Enhancement of Nitrate Anion UV Absorption

Abstract: We present simulated electronic absorption spectra of isolated and solvated nitrate anion in the UV region, focusing primarily on the absorption into the first absorption band around 300 nm. This weak absorption band in this spectral region is responsible for the generation of NOx in the polar areas or OH• radicals in the hydrosphere. The 300 nm absorption band is symmetrically strongly forbidden and coupling of at least two vibrational modes is needed to allow the transition in the isolated nitrate anion. Further symmetry breaking is provided by solvation. In this study we model the absorption spectra of nitrate–water clusters using the combined reflection principle path integral molecular dynamics (RP-PIMD) method. Condensed phase UV spectra are modeled within a cluster-continuum model. The calculated spectra are compared with experimental bulk phase measurements and reasonable agreement is found. We also provide a benchmarking of the DFT functionals to be used for a description of the electronically ex...

Classical and quantum ordering of protons in cold solid hydrogen under megabar pressures

Abstract: A combination of state-of-the-art theoretical methods has been used to obtain an atomic-level picture of classical and quantum ordering of protons in cold high-pressure solid hydrogen. We focus mostly on phases II and III of hydrogen, exploring the effects of quantum nuclear motion on certain features of these phases (through a number of ab initio path integral molecular dynamics (PIMD) simulations at particular points on the phase diagram). We also examine the importance of van der Waals forces in this system by performing calculations using the optB88-vdW density functional, which accounts for non-local correlations. Our calculations reveal that the transition between phases I and II is strongly quantum in nature, resulting from a competition between anisotropic inter-molecular interactions that restrict molecular rotation and thermal plus quantum fluctuations of the nuclear positions that facilitate it. The transition from phase II to III is more classical because quantum nuclear motion plays only a secondary role and the transition is determined primarily by the underlying potential energy surface. A structure of P21/c symmetry with 24 atoms in the primitive unit cell is found to be stable when anharmonic quantum nuclear vibrational motion is included at finite temperatures using the PIMD method. This structure gives a good account of the infra-red and Raman vibron frequencies of phase II. We find additional support for a C2/c structure as a strong candidate for phase III, since it remains transparent up to 300 GPa, even when quantum nuclear effects are included. Finally, we find that accounting for van der Waals forces improves the agreement between experiment and theory for the parts of the phase diagram considered, when compared to previous work which employed the widely-used Perdew–Burke–Ernzerhof exchange–correlation functional.

Investigations of the very short hydrogen bond in the crystal of nitromalonamide via Car-Parrinello and path integral molecular dynamics.

Abstract: In this paper are presented the results of theoretical studies of the structure in proton motion in a very short O···O and two weak N–H···O intramolecular hydrogen bonds in the nitromalonamide crystal. The dynamics of proton motion in hydrogen bonds were investigated in the NVT ensemble at 298 K using the Car–Parrinello and the path integral molecular dynamics. A very large delocalization of proton in the slightly asymmetrical single well of free energy potential of O–H···O intramolecular hydrogen bond was noted especially in the path integral simulation where quantum effects are taken into account. This hydrogen bond is very strong with the estimated energy of hydrogen bond ca. −27 kcal/mol. The nature of intra- and intermolecular interactions was studied by means of quantum theory of atoms in molecules. The infrared spectra were calculated and compared with available experimental data. CPMD vibrational results appear to be in good agreement with the experimental ones.

Simulation of quantum zero-point effects in water using a frequency-dependent thermostat

Abstract: Molecules like water have vibrational modes with a zero-point energy well above room temperature. As a consequence, classical molecular dynamics simulations of their liquids largely underestimate the energy of modes with a higher zero-point temperature, which translates into an underestimation of covalent interatomic distances due to anharmonic effects. Zero-point effects can be recovered using path integral molecular dynamics simulations, but these are computationally expensive, making their combination with ab initio molecular dynamics simulations a challenge. As an alternative to path integral methods, from a computationally simple perspective, one would envision the design of a thermostat capable of equilibrating and maintaining the different vibrational modes at their corresponding zero-point temperatures. Recently, Ceriotti et al. [Phys. Rev. Lett. 102, 020601 (2009)] introduced a framework to use a custom-tailored Langevin equation with correlated noise that can be used to include quantum fluctuations in classical molecular dynamics simulations. Here we show that it is possible to use the generalized Langevin equation with suppressed noise in combination with Nose-Hoover thermostats to efficiently impose a zero-point temperature on independent modes in liquid water. Using our simple and inexpensive method, we achieve excellent agreement for all atomic pair correlation functions compared to the path integral molecular dynamics simulation.

Boundary based on exchange symmetry theory for multilevel simulations. II. Multiple time scale approach

Abstract: The QM/MM BEST method presented in the first article of this series [M. Shiga and M. Masia, J. Chem. Phys. 139, 044120 (2013)] has been applied herein to simulate the whole series of hydrated alkali ions. In this article we show how to overcome the sampling bottleneck for QM/MM simulations by using our method with multiple time scale algorithm (MTS-BEST). We extend the use of MTS-BEST to ab initio QM/MM path integral molecular dynamics simulations, thus demonstrating that one could obtain a complete quantum description of the primary subsystem based on first principles. We highlight that the MTS-BEST approach could be generally applied to hybrid multiscale simulation of diffusive systems, thus extending its relevance to a broad class of simulation techniques beyond QM/MM. We show that it is important to account for electron correlation to better reproduce the hydration structural properties such as the ion-water radial distribution functions, and the anisotropic angular distributions around the ion.

INVITED ARTICLE Electronically coarse-grained molecular dynamics using quantum Drude oscillators

Abstract: Standard molecular dynamics (MD) simulations generally make use of a basic description of intermolecular forces which consists of fixed, pairwise, atom-centred Coulomb, van der Waals and short-range repulsive terms. Important interactions such as many-body polarisation and many-body dispersion which are sensitive to changes in the environment are usually neglected, and their effects treated effectively within mean-field approximations to reproduce a single thermodynamic state point or physical environment. This leads to difficulties in modelling the complex interfaces of interest today where the behaviour may be quite different from the regime of parameterisation. Here, we describe the construction and properties of a Gaussian coarse-grained electronic structure, which naturally generates many-body polarisation and dispersion interactions. The electronic structure arises from a fully quantum mechanical treatment of a set of distributed quantum Drude oscillators (QDOs), harmonic atoms which interact with each other and other moieties via electrostatic (Coulomb) interactions; this coarse-grained approach is capable of describing many-body polarisation and dispersion but not short-range interactions which must be parametrised. We describe how on-the-fly forces due to this exchange-free Gaussian model may be generated with linear scale in the number of atoms in the system using an adiabatic path integral molecular dynamics for quantum Drude oscillators technique (APIMD-QDO). We demonstrate the applicability of the QDO approach to realistic systems via a study of the liquid–vapour interface of water.

Temperature dependence on the structure of Zundel cation and its isotopomers

Abstract: Temperature dependence on the structural fluctuations of Zundel cation, H5O2+, and its isotopomers, D5O2+ and T5O2+, have been studied using path integral molecular dynamics simulations in which nuclear quantum effect is fully taken into account. It has been found that the fluctuations of hydrogen-oxygen and oxygen-oxygen distances, which are relevant to the hydrogen bonded structure, grow drastically as the temperature increases within the range of investigation between 100 K and 900 K. The fluctuation with respect to the position of non-bonded hydrogen also increases substantially as the temperature increases. The temperature dependence on the fluctuation is greater for D5O2+ or T5O2+ than that of H5O2+, since the zero-point effect of the former is less than the latter.

Nuclear quantum effects on protonated lysine with an asymmetric low barrier hydrogen bond: an ab initio path integral molecular dynamics study

Abstract: The nuclear quantum effect on the short and asymmetric hydrogen bond of protonated lysine (LysH+) at room temperature is explored by ab initio path integral molecular dynamics (PIMD) simulation. From static electronic structure calculations, the barrier height of proton transfer in LysH+ is 1.1 kcal mol−1, which is much lower than that of typical hydrogen bonds. The hydrogen-bonded proton is delocalized in between two nitrogen atoms in the PIMD simulation including both thermal and nuclear quantum effects, while the proton is localized on a nitrogen atom in a conventional ab initio molecular dynamics simulation including thermal effects alone. We found that the proton transfer barrier found in the static calculation and conventional ab initio simulation is completely washed out in the PIMD simulation. Meanwhile, the proton distribution at the Nζ atom was larger than that at the N atom, as found in the static calculation and conventional ab initio molecular dynamics simulation. We clarified that an asymmetric low barrier hydrogen bond exists in LysH+ at room temperature from our PIMD simulation.

Low-temperature anharmonicity of barium titanate: A path-integral molecular-dynamics study

Abstract: We investigate the influence of quantum effects on the dielectric and piezoelectric properties of barium titanate in its (low-temperature) rhombohedral phase, and show the strongly anharmonic character of this system even at low temperature. For this purpose, we perform path-integral molecular-dynamics simulations under fixed pressure and fixed temperature, using an efficient Langevin thermostat-barostat, and an effective Hamiltonian derived from first-principles calculations. The quantum fluctuations are shown to significantly enhance the static dielectric susceptibility ($\ensuremath{\approx}$ by a factor of 2) and the piezoelectric constants, reflecting the strong anharmonicity of this ferroelectric system even at very low temperature. The slow temperature-evolution of the dielectric properties observed below $\ensuremath{\approx}$100 K is attributed (i) to zero-point energy contributions and (ii) to harmonic behavior if the quantum effects are turned off.

On the origin and convergence of a Post-Quantization Constrained propagator for path integral simulations of rigid bodies.

Abstract: We present a new methodological procedure, based on Post-Quantization Constraints (PQC), to obtain approximate density matrices and energy estimators for use in path integral molecular dynamics and Monte Carlo simulations. The approach serves as a justification of the use of “RATTLE & SHAKE” type methods for path integrals. A thorough discussion of the underlying geometrical concepts is given. Two standard model systems, the particle on a ring and the three-dimensional linear rotor, are used to illustrate and benchmark the approach. In these two cases, matrix elements of the newly defined propagator are explicitly computed in both “angular coordinate” and “angular momentum” bases. A detailed analysis of the convergence properties of the density matrix, and energy estimator with respect to their “exact” counterparts, is presented along with numerical illustrations. We conclude that the use of a PQC-type propagator is justified and practical.

Ab Initio Molecular Dynamics Study of the Very Short O-H···O Hydrogen Bonds in the Condensed Phases.

Abstract: In this paper are presented the results of theoretical studies of the structure and proton motion in very short O···O intramolecular hydrogen bonds in two molecular crystals. A comparison was conducted between 3-cyano-2,4-pentanedione (I) and 4-cyano-2,2,6,6-tetramethyl-3,5-heptanedione (II) in the solid state. The dynamics of proton motion in the O–H···O hydrogen bond were investigated in he NVT ensemble at 298 and 50 K, respectively, for crystals I and II using Car–Parrinello and path integral molecular dynamics. Very large delocalization of the bridging proton was noted especially in the path integral simulation where quantum effects are taken into account. The infrared spectrum was calculated, and a comparative vibrational analysis was performed. CPMD vibrational results appear to be in qualitative agreement with the experimental ones.

Geometric isotope effects on small chloride ion water clusters with path integral molecular dynamics simulations

Abstract: The geometric isotope effects on the structures of hydrated chloride ionic hydrogen bonded clusters are explored by carrying out path integral molecular dynamics simulations. First, an outer shell coordinate is selected to display the rearrangement of single and multi hydration shell cluster structures. Next, to show the competition of intramolecular and intermolecular nuclear quantum effects, the intramolecular OH ∗ stretching and intermolecular ion–water wagging motions are studied for single and multi shell structures, respectively. The results indicate that the intermolecular nuclear quantum effects stabilize the ionic hydrogen bonds in single shell structures, while they are destabilized through the competition with intramolecular nuclear quantum effects in multi shell structures. In addition, the correlations between ion–water stretching motion and other cluster vibrational coordinates are discussed. The results indicate that the intermolecular nuclear quantum effects on the cluster structures are strongly related to the cooperation of the water–water hydrogen bond interactions.

Variational path integral molecular dynamics study of a water molecule

Abstract: In the present study, a variational path integral molecular dynamics method developed by the author [Chem. Phys. Lett. 482, 165 (2009)] is applied to a water molecule on the adiabatic potential energy surface. The method numerically generates an exact wavefunction using a trial wavefunction of the target system. It has been shown that even if a poor trial wavefunction is employed, the exact quantum distribution is numerically extracted, demonstrating the robustness of the variational path integral method.

Force-field functor theory: classical force-fields which reproduce equilibrium quantum distributions

Abstract: Feynman and Hibbs were the first to variationally determine an effective potential whose associated classical canonical ensemble approximates the exact quantum partition function. We examine the existence of a map between the local potential and an effective classical potential which matches the exact quantum equilibrium density and partition function. The usefulness of such a mapping rests in its ability to readily improve Born-Oppenheimer potentials for use with classical sampling. We show that such a map is unique and must exist. To explore the feasibility of using this result to improve classical molecular mechanics, we numerically produce a map from a library of randomly generated one-dimensional potential/effective potential pairs then evaluate its performance on independent test problems. We also apply the map to simulate liquid para-hydrogen, finding that the resulting radial pair distribution functions agree well with path integral Monte Carlo simulations. The surprising accessibility and transferability of the technique suggest a quantitative route to adapting Born-Oppenheimer potentials, with a motivation similar in spirit to the powerful ideas and approximations of density functional theory.

Spectroscopic fingerprints of toroidal nuclear quantum delocalization via ab initio path integral simulations.

Abstract: We investigate the quantum-mechanical delocalization of hydrogen in rotational symmetric molecular systems. To this purpose, we perform ab initio path integral molecular dynamics simulations of a methanol molecule to characterize the quantum properties of hydrogen atoms in a representative system by means of their real-space and momentum-space densities. In particular, we compute the spherically averaged momentum distribution n(k) and the pseudoangular momentum distribution n(kθ). We interpret our results by comparing them to path integral samplings of a bare proton in an ideal torus potential. We find that the hydroxyl hydrogen exhibits a toroidal delocalization, which leads to characteristic fingerprints in the line shapes of the momentum distributions. We can describe these specific spectroscopic patterns quantitatively and compute their onset as a function of temperature and potential energy landscape. The delocalization patterns in the projected momentum distribution provide a promising computational tool to address the intriguing phenomenon of quantum delocalization in condensed matter and its spectroscopic characterization. As the momentum distribution n(k) is also accessible through Nuclear Compton Scattering experiments, our results will help to interpret and understand future measurements more thoroughly. © 2012 Wiley Periodicals, Inc.

Response to Commentary on “Force-field functor theory: classical force-fields which reproduce equilibrium quantum distributions”

Abstract: Dr. Ruggero Vaia has made several claims (Vaia, 2013) about our recent submission entitled Force-field functor theory: classical force-fields which reproduce equilibrium quantum distributions (Babbush et al., 2013) which we would like to address. Dr. Vaia points out that the idea of an effective potential which reproduces quantum distributions is a well-known observation. We agree and point out that our paper contains many citations to the relevant literature. The contribution we make is to describe a new perspective on how the effective potential might be obtained in a spirit similar to Density Functional Theory (DFT). To make progress on a DFT-like theory, a proof of the existence and uniqueness of force-field functors is necessary. The main goal of our manuscript is to develop such proofs and to demonstrate and test the idea of functor on various test problems. After establishing uniqueness and existence proofs, we introduce an intentionally simple example of a functor in order to illustrate the utility of our broader approach. In particular, we use the Jensen–Peierls inequality to construct a variational approximation in which the mapping is linear. This procedure closely parallels the effective potential of Feynman and Hibbs (1965) in which the same approximation is used to model the effective potential as a linear convolution of the physical potential with a Gaussian of variance σ2=βℏ212 m. Dr. Vaia correctly shows that for a one-dimensional system, our analytical approximation reduces to a Gaussian smearing of the potential with exactly twice the width of the Feynman–Hibbs approximation. This is because the Feynman–Hibbs potential is a centroid effective potential which does not attempt to reproduce the correct particle distribution. While we make reference to this in our paper (we wrote “our [linear mapping] can be imagined as a Gaussian smearing of V(q)”), we do not discuss the explicit form of the mapping in one-dimension as this obscures the more general point we are trying to make. Linearity is not an essential feature of our approach; rather, it is an approximation which underlies a particular functor that we construct as a pedagogical example of how the functor approach could be used. In fact, the procedure we demonstrate to invert the linear mapping is significantly more accurate than the explicit Gaussian convolution which Dr. Vaia discusses. We used the bijective nature of the mapping (implied by our uniqueness proof) to empirically fit the optimal linear functor. This was accomplished by computing the exact effective potentials for 1000 random physical potentials and then using a least squares procedure with cross-validation to determine optimal matrix elements which map basis vectors of the physical potential to basis vectors of the effective potential. This cross-validation approach allowed us to test the accuracy of our method. We show that our functor works surprising well (especially considering its simplicity) and demonstrate its impressive performance on independent test problems including the simulation of liquid hydrogen at temperatures as low as 14 K. Though the error in this approximation grows as temperature decreases, the same is true of all path integral approaches since the number of required integration time slices grows without bound as temperature approaches zero. Degrading in quality at T = 0 does not render path integral molecular dynamics, or our technique useless. Dr. Vaia goes on to compare the linear approximation to the Feynman–Kleinert and Giachetti–Tognetti effective potentials and mentions that those approximations are more robust. While this is true, we are not trying to demonstrate the superiority of a particular explicit form for the functor under any particular approximation and we are certainly not interested in comparing the performance of this simple example to other paradigm methods in one dimension. The linear functor is for our theory analogous to what Local Density Approximation (LDA) is for DFT: a simplistic but illustrative example of the type of theory which could emerge from our uniqueness theorem. We have provided ample evidence for a new approach to obtaining effective potentials, in a similar spirit to DFT.

Low-temperature isomers of the water hexamer as predicted by q-TIP4P/F and TTM3-F potentials

Abstract: The hexamer is the smallest non-planar water cluster exhibiting several isomeric forms whose equilibria are governed by the subtle interplay between energetic and entropic contributions affected significantly by nuclear quantum effects. The hexamer thus can serve as a sensitive test for model potential energy surfaces (PESs), probing both the relative depth and shape of the minima corresponding to the different isomers. Here we report on the equilibrium populations of the hexamer isomers in the temperature range from 30 K to 150 K as predicted by two popular flexible force fields (q-TIP4P/F and TTM3-F) and compare them against the results obtained recently for an ab initio-based potential (WHBB). The calculations were performed using path-integral molecular dynamics combined with the replica exchange method.