Hybrid quantum/classical path integral approach for simulation of hydrogen transfer reactions in enzymes

A hybrid quantum/classical path integral Monte Carlo (QC-PIMC) method for calculating the quantum free energy barrier for hydrogen transfer reactions in condensed phases is presented. In this approach, the classical potential of mean force along a collective reaction coordinate is calculated using umbrella sampling techniques in conjunction with molecular dynamics trajectories propagated according to a mapping potential. The quantum contribution is determined for each configuration along the classical trajectory with path integral Monte Carlo calculations in which the beads move according to an effective mapping potential. This type of path integral calculation does not utilize the centroid constraint and can lead to more efficient sampling of the relevant region of conformational space than free-particle path integral sampling. The QC-PIMC method is computationally practical for large systems because the path integral sampling for the quantum nuclei is performed separately from the classical molecular dynamics sampling of the entire system. The utility of the QC-PIMC method is illustrated by an application to hydride transfer in the enzyme dihydrofolate reductase. A comparison of this method to the quantized classical path and grid-based methods for this system is presented.     

Rigid quantum Monte Carlo simulations of condensed molecular matter: water clusters in the n=2-->8 range

The numerical advantage of quantum Monte Carlo simulations of rigid bodies relative to the flexible simulations is investigated for some simple systems. The results show that if high frequency modes in molecular condensed matter are predominantly in the ground state, the convergence of path integral simulations becomes nonuniform. Rigid body quantum parallel tempering simulations are necessary to accurately capture thermodynamic phenomena in the temperature range where the dynamics are influenced by intermolecular degrees of freedom; the stereographic projection path integral adapted for quantum simulations of asymmetric tops is a significantly more efficient strategy compared with Cartesian coordinate simulations for molecular condensed matter under these conditions. The reweighted random series approach for stereographic path integral Monte Carlo is refined and implemented for the quantum simulation of water clusters treated as an assembly of rigid asymmetric tops.     

Rotational fluctuation of molecules in quantum clusters. I. Path integral hybrid Monte Carlo algorithm

In this paper, we present a path integral hybrid Monte Carlo (PIHMC) method for rotating molecules in quantum fluids. This is an extension of our PIHMC for correlated Bose fluids [S. Miura and J. Tanaka, J. Chem. Phys. 120, 2160 (2004)] to handle the molecular rotation quantum mechanically. A novel technique referred to be an effective potential of quantum rotation is introduced to incorporate the rotational degree of freedom in the path integral molecular dynamics or hybrid Monte Carlo algorithm. For a permutation move to satisfy Bose statistics, we devise a multilevel Metropolis method combined with a configurational-bias technique for efficiently sampling the permutation and the associated atomic coordinates. Then, we have applied the PIHMC to a helium-4 cluster doped with a carbonyl sulfide molecule. The effects of the quantum rotation on the solvation structure and energetics were examined. Translational and rotational fluctuations of the dopant in the superfluid cluster were also analyzed.     

Thermodynamic properties of small rare gas clusters by path integral Monte Carlo simulations: a preliminary study

A strategy for reducing the risk of non-ergodic simulations in Monte Carlo calculations of the thermodynamic properties of clusters is discussed with the support of some examples. The results obtained attest the significance of the approach for the low-temperature regime, as non-ergodic sampling of potential energy surfaces is a particularly insidious occurrence. Fourier path integral Monte Carlo techniques for taking into account quantum effects are adopted, in conjunction with suitable tricks for improving the procedure reliability. Applications are restricted to Lennard-Jones clusters of rare-gas systems.     

Path integral hybrid Monte Carlo algorithm for correlated Bose fluids

Path integral hybrid Monte Carlo (PIHMC) algorithm for strongly correlated Bose fluids has been developed. This is an extended version of our previous method [S. Miura and S. Okazaki, Chem. Phys. Lett. 308, 115 (1999)] applied to a model system consisting of noninteracting bosons. Our PIHMC method for the correlated Bose fluids is constituted of two trial moves to sample path-variables describing system coordinates along imaginary time and a permutation of particle labels giving a boundary condition with respect to imaginary time. The path-variables for a given permutation are generated by a hybrid Monte Carlo method based on path integral molecular dynamics techniques. Equations of motion for the path-variables are formulated on the basis of a collective coordinate representation of the path, staging variables, to enhance the sampling efficiency. The permutation sampling to satisfy Bose-Einstein statistics is performed using the multilevel Metropolis method developed by Ceperley and Pollock [Phys. Rev. Lett. 56, 351 (1986)]. Our PIHMC method has successfully been applied to liquid helium-4 at a state point where the system is in a superfluid phase. Parameters determining the sampling efficiency are optimized in such a way that correlation among successive PIHMC steps is minimized.     

A path integral influence functional for excess electron in fluids: Density-functional formulation

In this paper, we propose a path integral influence functional from a solvent to determine a self-correlation function of a quantum particle in classical simple fluid. It is shown that the influence functional is related to a grand potential functional of the pure solvent under a three-dimensional external field arising from a classical isomorphic polymer, on which the quantum particle is mapped. The influence functional can be calculated from the self-correlation function, the solute-solvent and the solvent-solvent pair correlation function. The obtained equation of the self-correlation function is applied to an excess electron problem in fluid helium. The Fourier path-integral Monte Carlo method is employed to perform the path integral of the electron. The solute-solvent pair correlation function is estimated from a reference interaction site model integral equation. These results obtained form our proposed influence functional and from that proposed by Chandler, Singh, and Richardson are compared with those provided by a path integral Monte Carlo simulation with the explicit helium solvent.     

The partial averaging Fourier path Integral approach based on the harmonic reference path

The Monte Carlo Fourier path integral approach has proved to be quite useful in calculating equilibrium thermodynamic properties. One of its advantages is that it can be formulated in such a way as to include higher order terms using the partial averaging technique, which includes the contribution from higher terms usually neglected by the discretized path integral approach. In the original approach, the Feynman path integral is evaluated via a free-particle reference state. Here, using a new expression for the Feynman paths expanded around a harmonic reference path, we derive an alternative formulation for the density matrix element and its corresponding partial averaging expression. Received: 16 September 1998 / Accepted: 30 October 1998 / Published online: 1 February 1999     

Boltzmann bias grand canonical Monte Carlo

We derive an efficient method for the insertion of structured particles in grand canonical Monte Carlo simulations of adsorption in very confining geometries. We extend this method to path integral simulations and use it to calculate the isotherm of adsorption of hydrogen isotopes in narrow carbon nanotubes (two-dimensional confinement) and slit pores (one-dimensional confinement) at the temperatures of 20 and 77 K, discussing its efficiency by comparison to the standard path integral grand canonical Monte Carlo algorithm. We use this algorithm to perform multicomponent simulations in order to calculate the hydrogen isotope selectivity for adsorption in narrow carbon nanotubes and slit pores at finite pressures. The algorithm described here can be applied to the study of adsorption of real oligomers and polymers in narrow pores and channels.     

The dynamo library for molecular simulations using hybrid quantum mechanical and molecular mechanical potentials

The Dynamo module library has been developed for the simulation of molecular systems using hybrid quantum mechanical (QM) and molecular mechanical (MM) potentials. Dynamo is not a program package but is a library of Fortran 90 modules that can be employed by those interested in writing their own programs for performing molecular simulations. The library supports a range of different types of molecular calculation including geometry optimizations, reaction‐path determinations and molecular dynamics and Monte Carlo simulations. This article outlines the general structure and capabilities of the library and describes in detail Dynamo's semiempirical QM/MM hybrid potential. Results are presented to indicate three particular aspects of this implementation—the handling of long‐range nonbonding interactions, the nature of the boundary between the quantum mechanical and molecular mechanical atoms and how to perform path‐integral hybrid‐potential molecular dynamics simulations. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 1088–1100, 2000     

Accurate vibrational-rotational partition functions and standard-state free energy values for H2O2 from Monte Carlo path-integral calculations

Accurate quantum mechanical partition functions and absolute free energies of H(2)O(2) are determined using a realistic potential energy surface [J. Koput, S. Carter, and N. C. Handy, J. Phys. Chem. A 102, 6325 (1998)] for temperatures ranging from 300 to 2,400 K by using Monte Carlo path integral calculations with new, efficient polyatomic importance sampling methods. The path centroids are sampled in Jacobi coordinates via a set of independent ziggurat schemes. The calculations employed enhanced-same-path extrapolation of trapezoidal Trotter Fourier path integrals, and the paths were constructed using fast Fourier sine transforms. Importance sampling was also used in Fourier coefficient space, and adaptively optimized stratified sampling was used in configuration space. The free energy values obtained from the path-integral calculations are compared to separable-mode approximations, to the Pitzer-Gwinn approximation, and to values in thermodynamic tables. Our calculations support the recently proposed revisions to the JANAF tables.     

Higher order and infinite Trotter-number extrapolations in path integral Monte Carlo

Improvements beyond the primitive approximation in the path integral Monte Carlo method are explored both in a model problem and in real systems. Two different strategies are studied: The Richardson extrapolation on top of the path integral Monte Carlo data and the Takahashi-Imada action. The Richardson extrapolation, mainly combined with the primitive action, always reduces the number-of-beads dependence, helps in determining the approach to the dominant power law behavior, and all without additional computational cost. The Takahashi-Imada action has been tested in two hard-core interacting quantum liquids at low temperature. The results obtained show that the fourth-order behavior near the asymptote is conserved, and that the use of this improved action reduces the computing time with respect to the primitive approximation.     

Connection between the observable and centroid structural properties of a quantum fluid: application to liquid para-hydrogen

It is shown that the discrepancy between path integral Monte Carlo [M. Zoppi et al., Phys. Rev. B 65, 092204 (2002)] and path integral centroid molecular dynamics [F. J. Bermejo et al., Phys. Rev. Lett. 84, 5359 (2000)] calculations of the static structure factor of liquid para-hydrogen can be explained based on a deconvolution equation connecting centroid and physical radial distribution functions. An explicit expression for the kernel of the deconvolution equation has been obtained using functional derivative techniques. In the superposition approximation, this kernel is given by the functional derivative of the effective potential with respect to the pairwise classical potential. Results of path integral Monte Carlo calculations for the radial distribution function and the static structure factor of liquid para-hydrogen are presented.     

Efficient estimators for quantum instanton evaluation of the kinetic isotope effects: application to the intramolecular hydrogen transfer in pentadiene

The quantum instanton approximation is used to compute kinetic isotope effects for intramolecular hydrogen transfer in cis-1,3-pentadiene. Due to the importance of skeleton motions, this system with 13 atoms is a simple prototype for hydrogen transfer in enzymatic reactions. The calculation is carried out using thermodynamic integration with respect to the mass of the isotopes and a path integral Monte Carlo evaluation of relevant thermodynamic quantities. Efficient "virial" estimators are derived for the logarithmic derivatives of the partition function and the delta-delta correlation functions. These estimators require significantly fewer Monte Carlo samples since their statistical error does not increase with the number of discrete time slices in the path integral. The calculation treats all 39 degrees of freedom quantum mechanically and uses an empirical valence bond potential based on a molecular mechanics force field.     

Iterative Monte Carlo formulation of real-time correlation functions

We present an iterative Monte Carlo path integral methodology for evaluating thermally averaged real-time correlation functions. Standard path integral Monte Carlo methods are used to sample paths along the imaginary time contour. Propagation of the density matrix is performed iteratively on a grid composed of the end points of the sampled paths. Minimally oscillatory propagators are constructed using energy filtering techniques. A single propagation yields the values of the correlation function at all intermediate time points. Model calculations suggest that the method yields accurate results over several oscillation periods and the statistical error grows slowly with increasing propagation time.     

Path integral molecular dynamics methods: Application to neon

Feynman's path integral formulation of quantum statistical mechanics, which has commonly been applied be Monte Carlo methods, is now also implemented by traditional molecular dynamics simulations of the microcanonical ensemble and in the Nosé-Hoover method simulating the isothermal-isobaric ensemble. In this article these two methods are applied to solid and liquid neon, in which quantum effects are not negligible. The validity of the procedure is shown by comparison with Monte Carlo and Brownian Dynamics computer simulations and with experiment. © 1995 by John Wiley & Sons, Inc.     

Path integral Monte Carlo method for ab initio calculation

The Feynman path integral method is applied to the many-electron problem. We first give new closure relations in terms of ordinary complex and real numbers, which could be derived from an arbitrary complete set of state vectors. Then, in the path integral form, the partition function of the system and the ensemble average of energy are explicitly expressed in terms of these closure relations. It is impossible to evaluate the path integral by direct numerical integrations because of its huge amount of integration variables. Therefore, we develop an algorithm by the Monte Carlo method with constraints corresponding to the normalization condition of states to calculate the required integral. Finally, the ensemble average of energy for the hydrogen molecule is explicitly evaluated by the quantum Monte Carlo method and results are compared with the result obtained by the ordinary full configuration interaction (CI) method. © 1996 John Wiley & Sons, Inc.     

Nuclear quantum effect on the hydrogen-bonded structure of guanine�Ccytosine pair

The structure of Watson?CCrick type guanine?Ccytosine (G?CC) base pair has been studied by classical hybrid Monte Carlo (HMC) and quantum path integral hybrid Monte Carlo (PIHMC) simulations on the semiempirical PM6 potential energy surface. For the three NH?X hydrogen-bonded moieties, the intramolecular NH bonds are found systematically longer while the H?X distance shorter in the PIHMC simulation than in the HMC simulation. We found that the hydrogen bonded length N?X correlates with the H?X distance, but not with the NH distance. A correlation is also between the neighboring hydrogen bonds in the G?CC base pair.     

Extrapolation of the time-step bias in diffusion quantum Monte Carlo by a differential sampling technique

A new method of eliminating the finite-time-step error inherent in diffusion quantum Monte Carlo is presented, utilizing an improved version of the existing differential techniques. An implementation is described and results of several small but representative calculations are discussed. The pertinent computation requirements on these systems were reduced by up to a factor of five by the new algorithm. It is speculated that this method may be easily applied to other quantum Monte Carlo and discretized path integral Monte Carlo techniques having related finite step-size errors with a possibility of obtaining similar good results.     

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.     

Path ensembles and path sampling in nonequilibrium stochastic systems

Markovian models based on the stochastic master equation are often encountered in single molecule dynamics, reaction networks, and nonequilibrium problems in chemistry, physics, and biology. An efficient and convenient method to simulate these systems is the kinetic Monte Carlo algorithm which generates continuous-time stochastic trajectories. We discuss an alternative simulation method based on sampling of stochastic paths. Utilizing known probabilities of stochastic paths, it is possible to apply Metropolis Monte Carlo in path space to generate a desired ensemble of stochastic paths. The method is a generalization of the path sampling idea to stochastic dynamics, and is especially suited for the analysis of rare paths which are not often produced in the standard kinetic Monte Carlo procedure. Two generic examples are presented to illustrate the methodology.