RPMDrate: Bimolecular chemical reaction rates from ring polymer molecular dynamics

Abstract

We present RPMDrate, a computer program for the calculation of gas phase bimolecular reaction rate coefficients using the ring polymer molecular dynamics (RPMD) method. The RPMD rate coefficient is calculated using the Bennett–Chandler method as a product of a static (centroid density quantum transition state theory (QTST) rate) and a dynamic (ring polymer transmission coefficient) factor. The computational procedure is general and can be used to treat bimolecular polyatomic reactions of any complexity in their full dimensionality. The program has been tested for the H+H2, H+CH4, OH+CH4 and H+C2H6 reactions. Program summaryProgram title: RPMDrateCatalogue identifier: AENW_v1_0Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AENW_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: MIT licenseNo. of lines in distributed program, including test data, etc.: 94512No. of bytes in distributed program, including test data, etc.: 1395674Distribution format: tar.gzProgramming language: Fortran 90/95, Python (version 2.6.x or later, including any version of Python 3, is recommended).Computer: Not computer specific.Operating system: Any for which Python, Fortran 90/95 compiler and the required external routines are available.Has the code been vectorized or parallelized?: The program can efficiently utilize 4096+ processors, depending on problem and available computer. At low temperatures, 110 processors are reasonable for a typical umbrella integration run with an analytic potential energy function and gradients on the latest x86-64 machines.RAM: 256 MbClassification: 16.12.External routines:•NumPy (http://numpy.scipy.org, version 1.5.0 or later is recommended).•FFTW3 (http://www.fftw.org, version 3.3 or later is recommended).Nature of problem:The RPMDrate program calculates thermal bimolecular rate coefficients of thermally activated atom–diatom and more complex bimolecular chemical reactions in the gas phase.Solution method:The RPMD rate is calculated using the Bennett–Chandler factorization as a product of a static (centroid density quantum transition state theory (QTST) rate) and a dynamic (transmission coefficient) factor. A key feature of this procedure is that it does not require that one calculate the absolute quantum mechanical partition function of the reactants or the transition state. The centroid density QTST rate is calculated from the potential of mean force along the reaction coordinate using umbrella integration. The reaction coordinate is taken to be an interpolating function that connects two dividing surfaces: one located in the asymptotic reactant valley and the other located in the transition state region. The Hessian of the collective reaction coordinate is obtained analytically. The transmission coefficient is calculated from the RPMD simulations with the hard constraint along the reaction coordinate.Restrictions:The applicability of RPMDrate is restricted to global potential energy surfaces with gradients. In the current release, they should be provided by Python callable objects.Unusual features: Simple and user-friendly input system provided by Python syntax.Additional comments:Test calculations for the H+H2 reactions were performed using the Boothroyd–Keogh–Martin–Peterson-2 (BKMP2) potential energy surface (PES) [1]. PESs for the H+CH4,OH+CH4 and H+C2H6 reactions are taken from the online POTLIB library [2]. PESs are included within the distribution package as Fortran subroutines. Implementations of the colored-noise, generalized Langevin equation (GLE) thermostats [3–5] have been included in the current release.The distribution contains example data and a detailed manual describing the use of RPMDrate.Running time:Highly dependent on the input parameters. The running time of RPMDrate depends mainly on the complexity of the potential energy surface and number of ring polymer beads. For the H+H2, H+CH4, and OH+CH4 test calculations given (with 128 ring polymer beads and analytic gradients), the running time is approximately 1800, 3600 and 4000 processor hours, respectively, on the Silicon Mechanics nServ A413 servers.