A stereographic projection path integral study of the coupling between the orientation and the bending degrees of freedom of water

A Monte Carlo path integral method to study the coupling between the rotation and bending degrees of freedom for water is developed. It is demonstrated that soft internal degrees of freedom that are not stretching in nature can be mapped with stereographic projection coordinates. For water, the bending coordinate is orthogonal to the stereographic projection coordinates used to map its orientation. Methods are developed to compute the classical and quantum Jacobian terms so that the proper infinitely stiff spring constant limit is recovered in the classical limit, and so that the nonconstant nature of the Riemann Cartan curvature scalar is properly accounted in the quantum simulations. The theory is used to investigate the effects of the geometric coupling between the bending and the rotating degrees of freedom for the water monomer in an external field in the 250 to 500 K range. We detect no evidence of geometric coupling between the bending degree of freedom and the orientations.     

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.     

Introducing phase transitions to quantum chemistry: from Trouton's rule to first principles vaporization entropies

In the present study, we employ quantum cluster equilibrium calculations on a small water cluster set in order to derive thermochemical equilibrium properties of the liquid phase as well as the liquid-vapor phase transition. The focus is set on the calculation of liquid phase entropies, from which entropies of vaporization at the normal boiling point of water are derived. Different electronic structure methods are compared and the influences of basis set size and of cooperative effects are discussed. In line with a previous study on the subject [B. Kirchner, J. Chem. Phys. 123, 204116 (2005)], we find that the neglect of cooperativity leads to large errors in the equilibrium cluster populations as well as in the obtained entropy values. In contrast, a correct treatment of the intermolecular many-body interaction yields liquid phase entropies and phase transition entropies being in very good agreement with the experimental reference, thus demonstrating that the quantum cluster equilibrium partition function intrinsically accounts for the shortcomings of the ideal gas partition function often employed in first principles entropy calculations. Comparing the calculated vaporization entropies to the value predicted by Trouton's rule, it is observed that for entropy calculations the consideration of intracluster cooperative effects is more important than the explicit treatment of the intercluster association even in a highly associated liquid such as water. The decomposition of entropy into contributions due to different degrees of freedom implies the need for the accurate treatment of particle indistinguishability and free volume of translation, whereas minor influences should be expected from the vibrational and rotational degrees of freedom and none from the electronic degrees of freedom.     

On the calculation of absolute free energies from molecular-dynamics or Monte Carlo data

We propose a method for calculating absolute free energies from Monte Carlo or molecular-dynamics data. The method is based on the identity that expresses the partition function Q as a Boltzmann average: 1Q=w(p,x)exp[betaH(p,x)], where w(p,x) is an arbitrary weight function such that its integral over the phase space is equal to 1. In practice, to minimize statistical errors the weight function is chosen such that the regions of the phase space where sampling statistics are poor are excluded from the average. The "ideal" weight function would be the equilibrium phase-space density exp[-betaH(p,x)]Q itself. We consider two methods for constructing the weight function based on different estimates of the equilibrium phase-space density from simulation data. In the first method, it is chosen to be a Gaussian function, whose parameters are obtained from the covariance matrix of the atomic coordinates. In the second, a clustering algorithm is used to attempt partitioning the data into clusters corresponding to different basins of attraction visited by the system. The weight function is then constructed as a superposition of Gaussians calculated for each cluster separately. We show that these strategies can be used to improve upon previous methods of estimating absolute entropies from covariance matrices.     

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.     

Efficient global biopolymer sampling with end-transfer configurational bias Monte Carlo

We develop an "end-transfer configurational bias Monte Carlo" method for efficient thermodynamic sampling of complex biopolymers and assess its performance on a mesoscale model of chromatin (oligonucleosome) at different salt conditions compared to other Monte Carlo moves. Our method extends traditional configurational bias by deleting a repeating motif (monomer) from one end of the biopolymer and regrowing it at the opposite end using the standard Rosenbluth scheme. The method's sampling efficiency compared to local moves, pivot rotations, and standard configurational bias is assessed by parameters relating to translational, rotational, and internal degrees of freedom of the oligonucleosome. Our results show that the end-transfer method is superior in sampling every degree of freedom of the oligonucleosomes over other methods at high salt concentrations (weak electrostatics) but worse than the pivot rotations in terms of sampling internal and rotational sampling at low-to-moderate salt concentrations (strong electrostatics). Under all conditions investigated, however, the end-transfer method is several orders of magnitude more efficient than the standard configurational bias approach. This is because the characteristic sampling time of the innermost oligonucleosome motif scales quadratically with the length of the oligonucleosomes for the end-transfer method while it scales exponentially for the traditional configurational-bias method. Thus, the method we propose can significantly improve performance for global biomolecular applications, especially in condensed systems with weak nonbonded interactions and may be combined with local enhancements to improve local sampling.     

Path integral Monte Carlo approach for weakly bound van der Waals complexes with rotations: algorithm and benchmark calculations

A path integral Monte Carlo technique suitable for the treatment of doped helium clusters with inclusion of the rotational degrees of freedom of the dopant is introduced. The extrapolation of the results to the limit of infinite Trotter number is discussed in detail. Benchmark calculations for small weakly bound (4)He(N)--OCS clusters are presented. The Monte Carlo results are compared with those of basis set calculations for the He--OCS dimer. A technique to analyze the orientational imaginary time correlation function is suggested. It allows one to obtain information regarding the effective rotational constant for a doped helium cluster based on a model for the rotational Hamiltonian. The renormalization of the effective rotational constant for (4)He(N)--OCS clusters derived from the orientational imaginary time correlation function is in good agreement with experimental results.     

Path integral methods for rotating molecules in superfluids

We present a path integral Monte Carlo (PIMC) methodology for quantum simulation of molecular rotations in superfluid environments such as helium and para-hydrogen that combines the sampling of rotational degrees of freedom for a molecular impurity with multilevel Metropolis sampling of Bose permutation exchanges for the solvating species. We show how the present methodology can be applied to the evaluation of imaginary time rotational correlation functions of the molecular impurity, from which the effective rotational constants can be extracted. The combined rotation/permutation sampling approach allows for the first time explicit assessment of the effect of Bose permutations on molecular rotation dynamics, and the converse, i.e., the effect of molecular rotations on permutation exchanges and local superfluidity. We present detailed studies showing that the effect of Bose permutations in the solvating environment is more significant for the dynamics of heavy than light molecules in helium, and that Bose permutation exchanges are slightly enhanced locally by molecular rotation. Finally, the examples studied here reveal a size dependence of rotational excitations for molecules possessing a strongly anisotropic interaction with helium in 4HeN clusters between N approximately 20 and N approximately 10(3).     

Quantum mechanical contributions to the structure of liquid water

Path integral Monte Carlo methods are used to study the effect of quantization of the orientational degrees of freedom of water (H₂O), using the ST2 model. A comparison of the classical and quantum atom—atom intermolecular correlation functions show that significant quantitative effects are manifest in the results.     

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.     

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.     

Monte Carlo method for computing density of states and quench probability of potential energy and enthalpy landscapes

The thermodynamics and kinetics of a many-body system can be described in terms of a potential energy landscape in multidimensional configuration space. The partition function of such a landscape can be written in terms of a density of states, which can be computed using a variety of Monte Carlo techniques. In this paper, a new self-consistent Monte Carlo method for computing density of states is described that uses importance sampling and a multiplicative update factor to achieve rapid convergence. The technique is then applied to compute the equilibrium quench probability of the various inherent structures (minima) in the landscape. The quench probability depends on both the potential energy of the inherent structure and the volume of its corresponding basin in configuration space. Finally, the methodology is extended to the isothermal-isobaric ensemble in order to compute inherent structure quench probabilities in an enthalpy landscape.     

Quantum monte carlo methods for constrained systems

The torsional ground state for ethane, the torsional, rotational, and mixed torsional and rotational ground state of propane are computed with a version of diffusion Monte Carlo adapted to handle the geometric complexity of curved spaces such as the Ramachandra space. The quantum NVT ensemble average for the mixed torsional and rotational degrees of freedom of propane is computed, using a version of Monte Carlo path integral, also adapted to handle curved spaces. These three problems are selected to demonstrate the generality and the applicability of the approaches described. The spaces of coordinates can be best constructed from the parameters of continuous Lie groups, and alternative methods based on vector spaces, where extended Lagrangian terms would be too cumbersome to implement. We note that the geometric coupling between the torsions and the rotations of propane produces a substantial effect on the ground state energy of propane, and that the quantum effects on the energy of propane are quite large even well above room temperature. © 2014 Wiley Periodicals, Inc.     

Computational methods in coupled electron-ion Monte Carlo simulations.

In the last few years, we have been developing a Monte Carlo simulation method to cope with systems of many electrons and ions in the Born-Oppenheimer approximation: the coupled electron-ion Monte Carlo method (CEIMC). Electronic properties in CEIMC are computed by quantum Monte Carlo rather than by density functional theory (DFT) based techniques. CEIMC can, in principle, overcome some of the limitations of the present DFT-based ab initio dynamical methods. The new method has recently been applied to high-pressure metallic hydrogen. Herein, we present a new sampling algorithm that we have developed in the framework of the reptation quantum Monte Carlo method chosen to sample the electronic degrees of freedom, thereby improving its efficiency. Moreover, we show herein that, at least for the case of metallic hydrogen, variational estimates of the electronic energies lead to an accurate sampling of the proton degrees of freedom.     

Study of the polymer interaction with a surface by entropic sampling

By means of a variant of the Monte Carlo method (entropic sampling within the Wang-Landau algorithm) the models of the interaction of a neutral polymer with a flat surface are studied. The method yields distribution functions over the energy and the distance from the polymer to the surface. Based on these distributions, excess entropies of the systems and their thermal properties are calculated: internal energy, heat capacity, average radius of gyration, average chain end-to-end distance, and average distance from the polymer to the surface. Continuous and lattice models are considered.     

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.     

Estimation of the absolute internal-rotation entropy of molecules with two torsional degrees of freedom from stochastic simulations

A method of statistical estimation is applied to the problem of evaluating the absolute entropy of internal rotation in a molecule with two torsional degrees of freedom. The configurational part of the entropy is obtained as that of the joint probability density of an arbitrary form represented by a two-dimensional Fourier series, the coefficients of which are statistically estimated using a sample of the torsional angles of the molecule obtained by a stochastic simulation. The internal rotors in the molecule are assumed to be attached to a common frame, and their reduced moments of inertia are initially calculated as functions of the two torsional angles, but averaged over all the remaining internal degrees of freedom using the stochastic-simulation sample of the atomic configurations of the molecule. The torsional-angle dependence of the reduced moments of inertia can be also averaged out, and the absolute internal-rotation entropy of the molecule is obtained in a good approximation as the sum of the configurational entropy and a kinetic contribution fully determined by the averaged reduced moments of inertia. The method is illustrated using Monte Carlo simulations of isomers of stilbene and halogenated derivatives of propane. The two torsional angles in cis-stilbene are found to be much more strongly correlated than those in trans-stilbene, while the degree of the angular correlation in propane increases strongly on substitution of hydrogen atoms with chlorine.     

Replica exchange and expanded ensemble simulations as Gibbs sampling: simple improvements for enhanced mixing

The widespread popularity of replica exchange and expanded ensemble algorithms for simulating complex molecular systems in chemistry and biophysics has generated much interest in discovering new ways to enhance the phase space mixing of these protocols in order to improve sampling of uncorrelated configurations. Here, we demonstrate how both of these classes of algorithms can be considered as special cases of Gibbs sampling within a Markov chain Monte Carlo framework. Gibbs sampling is a well-studied scheme in the field of statistical inference in which different random variables are alternately updated from conditional distributions. While the update of the conformational degrees of freedom by Metropolis Monte Carlo or molecular dynamics unavoidably generates correlated samples, we show how judicious updating of the thermodynamic state indices--corresponding to thermodynamic parameters such as temperature or alchemical coupling variables--can substantially increase mixing while still sampling from the desired distributions. We show how state update methods in common use can lead to suboptimal mixing, and present some simple, inexpensive alternatives that can increase mixing of the overall Markov chain, reducing simulation times necessary to obtain estimates of the desired precision. These improved schemes are demonstrated for several common applications, including an alchemical expanded ensemble simulation, parallel tempering, and multidimensional replica exchange umbrella sampling.     

A critical comparison of approximation methods and models for equilibrium properties of low-barrier hydrogen bonds

Recent experimental evidence has led to the conclusion that short, strong hydrogen bonds can stabilize transition states of enzyme catalyzed biochemical reactions. Evidence for such hydrogen bonds is the low value of the isotopic fractionation factor, phi, which is defined as the equilibrium constant for the generic reaction, R-H + DOH <--> R-D + HOH, where H is the hydrogen atom participating in the low-barrier hydrogen bond in a molecule R-H. In this work we assess two approximation methods for computing the isotopic fractionation factors for single and multidimensional systems containing a low-barrier hydrogen bond. These methods are WKB and an approach that corrects the classical partition function via a quantum correction factor. We find that the latter approach is universally accurate and applicable in both single and multidimensional systems containing a low-barrier hydrogen bond. We also assess two different models for the coupling of a molecule's low-barrier hydrogen bond to other degrees of freedom, both internal and external to the molecule, and show that each leads to a lowering of the fractionation factor.     

Application of the Monte Carlo method to modeling of kinetic processes with the participation of polyatomic molecules and clusters: 1. Kinetics of energy transfer during collisions of polyatomic molecules

The Monte Carlo method was used to model the collisional energy transfer for polyatomic molecules within the framework of the statistical theory of reactions. A model describing energy transfer through the formation of a statistical collisional complex was suggested. It was assumed that the total energy of the complex was randomized in the course of collisions and statistically distributed among the internal and translational degrees of freedom. The method was verified by comparing the equilibrium distribution functions for the vibrational, rotational, and total energies of the molecule. The mean energy portion and the root-mean-square energy portion transferred per collision, as functions of the total molecular energy, were determined. The relaxation parameters of the population distribution over energy after a sharp increase in the bath-gas temperature were calculated.