Polyrate 2023: A computer program for the calculation of chemical reaction rates for polyatomics. New version announcement

Polyrate is a suite of computer programs for the calculation of chemical reaction rates of polyatomic species (including atoms and diatoms as special cases) by variational transition state theory (VTST); conventional transition state theory is also supported. Polyrate can calculate the rate constants for both bimolecular reactions and unimolecular reactions , and it can be applied to reactions in the gas phase, liquid solution phase, or solid state and to reactions at gas–solid interfaces. Polyrate can perform VTST calculations on gas-phase reactions with both tight and loose transition states. For tight transition states it uses the reaction-path (RP) variational transition state theory developed by Garrett and Truhlar, and for loose transition states it uses variable-reaction-coordinate (VRC) variational transition state theory developed by Georgievskii and Klippenstein. The RP methods used for tight transition states are conventional transition state theory, canonical variational transition state theory (CVT), and microcanonical variational transition state theory (μVT) with multidimensional semiclassical approximations for tunneling and nonclassical reflection. For VRC calculations, rate constants may be calculated for canonical or microcanonical ensembles or energy- and total-angular-momentum resolved microcanonical ensembles. Pressure-dependent rate constants for elementary reactions can be computed using system-specific quantum RRK theory (SS-QRRK) with the information obtained from high-pressure-limit VTST calculation as input by using the SS-QRRK utility code. Alternatively, Polyrate 2023 may be interfaced with TUMME 2023 for a master-equation treatment of pressure dependence or to obtain phenomenological rate constants for complex mechanisms. Potential energy surfaces may be analytic functions evaluated by subroutines, or they may be implicit surfaces defined by electronic structure input files or interface subroutines containing energies, gradients, and force constants (Hessians). For the latter, Polyrate can be used in conjunction with various interfaces to electronic structure programs for direct dynamics, and it has routines designed to make such interfacing straightforward. Polyrate supports six options for direct dynamics, namely (i) straight single-level direct dynamics, (ii) zero-order interpolated variational transition state theory (IVTST-0), (iii) first-order interpolated variational transition state theory (IVTST-1), (iv) interpolated variational transition state theory by mapping (IVTST-M), (v) variational transition state theory with interpolated single-point energies (VTST-ISPE), and (vi) variational transition state theory with interpolated optimized corrections (VTST-IOC). Polyrate can handle multistructural and torsional-potential anharmonicity in conjunction with the MSTor program. Polyrate 2023 contains 112 test runs, and 46 of these are for direct dynamics calculations; 85 of the test runs are single-level runs, and 27 are dual-level calculations. Program Title: Polyrate 2023 CPC Library link to program files: https://doi.org/10.17632/vpf3y9cbhf.1 Developer's repository link: https://doi.org/10.5281/zenodo.8213313 Licensing Provisions: Apache-2.0 for the program and CC-BY-4.0 for the manual Programming language: Fortran 90 Does the new version supersede the previous version?: Yes Journal references of previous versions: Comput. Phys. Commun. 47 (1987) 91 [1]; Comput. Phys. Commun. 71 (1992) 235 [2]; Comput. Phys. Commun. 88 (1995) 341 [3]. Nature of problem: Polyrate 2023 can calculate unimolecular and bimolecular chemical reaction rates of polyatomic species (and also of atoms and diatoms as special cases) by using variational transition state theory with multidimensional tunneling [4-8] or by conventional transition state theory [9]. The program uses either reaction-path variational transition state theory for tight transition states or variable-reaction-coordinate variational transition state theory for loose transition states. Rate constants may be calculated for canonical or microcanonical ensembles with multidimensional semiclassical approximations for tunneling and nonclassical reflection and for reactions in the gas phase or solid state or at gas–solid interfaces. In addition, rates may be calculated for vibrationally adiabatic and vibrationally diabatic state-selected reactions in which one or more of the vibrational modes is restricted to the ground state or to the first excited state, while the translational, rotational, and remaining vibrational modes are treated thermally [10,11]. The program may also be used to find stationary geometries of reactants, products, and transition states and to calculate reaction paths, Arrhenius parameters, and equilibrium constants . Solution method: Energies, gradients, and Hessians are obtained from analytic potential energy functions or from electronic structure input data, where the input data may be supplied in files or by direct interface to electronic structure codes. A number of methods are available for calculating the reaction path. Nontorsional anharmonicity may be included under the assumption of independent generalized normal modes. Torsional potential anharmonicity and multiple-structure anharmonicity can be included along the reaction path using the formulas of the multi-structural torsional (MS-T) method with one input structure [12-16]. Small-curvature [2,17], large-curvature [18,19], and optimized multidimensional [20,21] semiclassical tunneling methods are available as options to calculate tunneling contributions to transmission coefficients. Pressure-dependent rate constants for elementary reactions can be computed using system-specific quantum Rice–Ramsperger–Kassel theory (SS-QRRK) [22]. Phenomenological rate constants for complex reactions can be obtained in conjunction with the TUMME 2023 program [23-25]. Interfaces to other programs: Polyrate 2023 can run as a standalone program, and many interface options are also available [24], including: CRATE, GAMESSPLUSRATE, GAUSSRATE, JAGUARATE, MC-TINKERRATE, MORATE, MULTILEVELRATE, NWCHEMRATE . Reason for new version: The new version greatly increases the functionality, stability, and accuracy; the code has been cleaned up, and some bugs have been fixed. Summary of revisions: 1. The variable-reaction-coordinate method with either single-faceted or multi-faceted dividing surfaces [26–30] has been implemented for barrierless association reactions and dissociation reactions without an intrinsic barrier, i.e., reactions without a barrier when considered in the downhill direction. This implementation is parallelized using MPI . In calculating the total number of available states for this method, the integration over the projection on the principal axis of the total angular momentum is calculated analytically. The user now has a choice of using the original variational transition state theory based on a reaction path (VTST-RP) or using variational transition state theory based on a variable reaction coordinate (VTST-VRC). 2. Options to carry out VTST-RP calculations with generalized normal modes in curvilinear coordinates, in particular in nonredundant [31–33] or redundant [34] internal coordinates, were added. All standard force constant types (stretches, nondegenerate bends, torsions, and degenerate bends) are available. Curvilinear coordinates are recommended for all calculations based on reaction paths. 3. The number-of-states calculation that is required for microcanonical variational theory (μVT) [35] is now programmed much more efficiently by the extended Beyer-Swinehart method [36,37]. Now there is no significant difference between CVT and μVT as far as computational cost. 4. The improved canonical variational theory (ICVT) option is removed. Whenever one needs a more accurate treatment of the threshold, one can simply use μVT without significantly increasing the computational cost, so there is no longer any motivation for using ICVT. 5. The first-order interpolated variational transition state theory (IVTST-1) method [38], interpolated variational transition state theory by mapping (IVTST-M) [39,40], including a spline-under-tension option, and several variational transition state theory with interpolated corrections (VTST-IC) [40] algorithms were added. IVTST-M allows arbitrary numbers of points along the reaction path and arbitrary numbers of points at which Hessians are computed, and it leads to faster convergence of direct dynamics calculations. VTST-IC methods include variational transition state theory with interpolated single-point energies (VTST-ISPE). 6. Interfaces for direct dynamics were made more convenient. The PTUMME keyword to print an interface file readable by the TUMME 2023 program [23-25] has been added. 7. An option for performing calculations with classical rather than quantal partition functions was added. 8. The code has been modified to calculate wells in the reactant and/or product channels. One can optimize the reactants, the products, and the saddle point in separate runs. Quasi-Newton optimization with BFGS updating [41–44] was added. The eigenvector following (EF) algorithm [45–47] was added for geometry optimization of transition states and other stationary states. 9. The canonical unified statistical model (CUS) [6,48–50] was added. 10. The re-orientation of the dividing surface (RODS) [51,52] and variational reaction path (VRP) [53] algorithms were implemented. These methods allow VTST calculations with larger step sizes. 11. An option was added to perform nonequilibrium solvation calculations [54]. 12. Version 4 of the large-curvature ground state tunneling method (LCG4) [55] was added as the default large-curvature tunneling (LCT) method, and it is now also used in the optimized multidimensional tunneling (OMT) [56] and microcanonical optimized multidimensional (μOMT) [20] tunneling methods. The interpolated algorithm for efficient large-curvature tunneling calculations (ILCT) [6,35] was also added. 13. New WKB vibrational anharmonicity [58,59] capability was added. 14. An option was added to scale the vibrational frequencies [60,61]. 15. A keyword was added to allow users to raise low frequencies to a certain number (e.g., 50 or 100 cm −1 ) to approximately account for anharmonicity [62,63]. 16. Several separable options for including anharmonicity in hindered internal rotations were added [64]. 17. The single-structure MS-T method [12,14] was added with coupled torsional potentials [16]. This capability is added both for a saddle point and for including torsional anharmonicity along a reaction path. Multi-structural [13] and multi-path [65–67] variational transition state theory calculations can be performed by combining the data from multiple Polyrate 2023 and MSTor outputs. 18. The SS-QRRK code for pressure dependent rate constants [22,68,69] has been included as a utility code. This method can include variational effects, multidimensional tunneling, and various vibrational and torsional anharmonicity options into the computed pressure dependent rate constants and the SS-QRRK microcanonical rate constants. The FE integral for the normalized Boltzmann-weighted number of unimolecular states above the threshold energy can be evaluated via numerical integration instead of using the approximate analytical formula [70]. 19. The standard-state free energy of activation profile is calculated along a reaction path using a standard state of 1 bar rather than an unconventional pressure. 20. The generalized free energy of activation can now be calculated using logarithms of the vibrational partition functions to avoid underflow. Similarly, a shifting constant is used to solve the overflow and underflow problem in the case that the imaginary action integral (the semiclassical tunneling integral) is too large. 21. The code was rearranged to allow parallel processing using MPI. 22. All common blocks were replaced by Fortran 90 modules. Memory for arrays is dynamically allocated. Therefore, compiling the code with a predefined system size (e.g., number of atoms) is no longer required. 23. Various bug fixes have been made, and some errors have been corrected. The restart option was improved. The code was made more portable, many options were made more convenient, defaults or keywords were improved, the output and manual were improved, and many new test runs were added. [1] A.D. Isaacson, D.G. Truhlar, S.N. Rai, R. Steckler, G.C. Hancock, B.C. Garrett, M.J. Redmon, Comput. Phys. Commun. 47 (1987) 91–102, https://doi.org/10.1016/0010-4655(87)90069-5 . [2] D.-h. Lu, T.N. Truong, V.S. Melissas, G.C. Lynch, Y.-P. Liu, B.C. Garrett, R. Steckler, A.D. Isaacson, S.N. Rai, G.C. Hancock, J.G. Lauderdale, T. Joseph, D.G. Truhlar, Comput. Phys. Commun. 71 (1992) 235-262, https://doi.org/10.1016/0010-4655(92)90012-N . [3] R. Steckler, W.-P. Hu, Y.-P. Liu, G.C. Lynch, B.C. Garrett, A.D. Isaacson, V.S. Melissas, D.-H. Lu, T.N. Truong, S.N. Rai, G.C. Hancock, J.G. Lauderdale, T. Joseph, D.G. Truhlar, Comput. Phys. Commun. 88 (1995) 341-343, https://doi.org/10.1016/0010-4655(95)00037-G . [4] D.G. Truhlar, B.C. Garrett, Acc. Chem. Res. 13 (1980) 440-448, https://doi.org/10.1021/ar50156a002 . [5] D.G. Truhlar, B.C. Garrett, Annu. Rev. Phys. Chem. 35 (1984) 159-189, https://doi.org/10.1146/annurev.pc.35.100184.001111 . [6] D.G. Truhlar, A.D. Isaacson, B.C. Garrett, in: M. Baer (Ed.), Theory of Chemical Reaction Dynamics, vol. 4, CRC Press, Boca Raton, FL, 1985, pp. 65-137. [7] A. Fernandez-Ramos, B.A. Ellingson, B.C. Garrett, D.G. Truhlar, Rev. Comput. Chem. 23 (2007) 125-232, https://doi.org/10.1002/9780470116449.ch3 . [8] J.L. Bao, D.G. Truhlar, Chem. Soc. Rev. 46 (2017) 7548-7596, https://doi.org/10.1039/C7CS00602K . [9] S. Glasstone, K.J. Laidler, H. Eyring, The Theory of Rate Processes, McGraw-Hill Book Co., New York, 1941. [10] D.G. Truhlar, A.D. Isaacson, J. Chem. Phys. 77 (1982) 3516-3522, https://doi.org/10.1063/1.444297 . [11] R. Steckler, D.G. Truhlar, B.C. Garrett, J. Chem. Phys. 84 (1986) 6712-6718, https://doi.org/10.1063/1.450673 . [12] J. Zheng, T. Yu, E. Papajak, I.M. Alecu, S.L. Mielke, D.G. Truhlar, Phys. Chem. Chem. Phys. 13 (2011) 10885-10907, https://doi.org/10.1039/C0CP02644A . [13] T. Yu, J. Zheng, D.G. Truhlar, Chem. Sci. 2 (2011) 2199-2213, https://doi.org/10.1039/C1SC00225B . [14] J. Zheng, S.L. Mielke, K.L. Clarkson, D.G. Truhlar, Comput. Phys. Commun. 183 (2012) 1803-1812, https://doi.org/10.1016/j.cpc.2012.03.007 . [15] W. Chen, J. Zheng, J.L. Bao, D.G. Truhlar, X. Xu, Comput. Phys. Commun. 288 (2023) 108740, https://doi.org/10.1016/j.cpc.2023.108740 . [16] J. Zheng, D.G. Truhlar, J. Chem. Theory Comput. 9 (2013) 1356-1367, https://doi.org/10.1021/ct3010722 . [17] Y.-P. Liu, G.C. Lynch, T.N. Truong, D.-h. Lu, D.G. Truhlar, B.C. Garrett, J. Am. Chem. Soc. 115 (1993) 2408-2415, https://doi.org/10.1021/ja00059a041 . [18] A. Fernandez-Ramos, D.G. Truhlar, J. Chem. Phys. 114 (2001) 1491-1496, https://doi.org/10.1063/1.1329893 . [19] A. Fernández-Ramos, D.G. Truhlar, J. Chem. Theory Comput. 1 (2005) 1063-1078, https://doi.org/10.1021/ct050153i . [20] Y.-P. Liu, D.-h. Lu, A. Gonzàlez-Lafont, D.G. Truhlar, B.C. Garrett, J. Am. Chem. Soc. 115 (1993) 7806-7817, https://doi.org/10.1021/ja00070a029 . [21] A. Fernández-Ramos, J.A. Miller, S.J. Klippenstein, D.G. Truhlar, Chem. Rev. 106 (2006) 4518-4584, https://doi.org/10.1021/cr050205w . [22] J.L. Bao, J. Zheng, D.G. Truhlar, J. Am. Chem. Soc. 138 (2016) 2690-2704, https://doi.org/10.1021/jacs.5b11938 . [23] R.M. Zhang, X. Xu, Donald G. Truhlar, TUMME 2023: An open-source master equation program for chemical kinetics , Zenodo, https://doi.org/10.5281/zenodo.7943284 . [24] R.M. Zhang, X. Xu, Donald G. Truhlar, Comput. Phys. Commun. 270 (2022) 108140, http://doi.org/10.1016/j.cpc.2021.108140 [25] R.M. Zhang, X. Xu, Donald G. Truhlar, Comput. Phys. Commun. 293 (2023) 108894, https://doi.org/10.1016/j.cpc.2023.108894 . [26] Y. Georgievskii, S.J. Klippenstein, J. Chem. Phys. 118 (2003) 5442-5455, https://doi.org/10.1063/1.1539035 . [27] Y. Georgievskii, S.J. Klippenstein, J. Phys. Chem. A 107 (2003) 9776-9781, https://doi.org/10.1021/jp034564b . [28] J. Zheng, S. Zhang, D.G. Truhlar, J. Phys. Chem. A 112 (2008) 11509-11513, https://doi.org/10.1021/jp806617m . [29] J.L. Bao, X. Zhang, D.G. Truhlar, Proc. Nat. Acad. Sci. USA 113 (2016) 13606-13611, https://doi.org/10.1073/pnas.1616208113 . [30] L. Zhang, D.G. Truhlar, S. Sun, Proc. Nat. Acad. Sci. USA 117 (2020) 5610-5616, https://doi.org/10.1073/pnas.1920018117 . [31] C.F. Jackels, Z. Gu, D.G. Truhlar, J. Chem. Phys. 102 (1995) 3188-3201, https://doi.org/10.1063/1.468630 . [32] K.A. Nguyen, C.F. Jackels, D.G. Truhlar, J. Chem. Phys. 104 (1996) 6491-6496, https://doi.org/10.1063/1.471369 . [33] Y.-Y. Chuang, D.G. Truhlar, J. Chem. Phys. 107 (1997) 83-89, https://doi.org/10.1063/1.474377 . [34] Y.-Y. Chuang, D.G. Truhlar, J. Phys. Chem. A 102 (1998) 242-247, https://doi.org/10.1021/jp9724028 . [35] A.D. Isaacson, M.T. Sund, S.N. Rai, D.G. Truhlar, J. Chem. Phys. 82 (1985) 1338-1340, https://doi.org/10.1063/1.448963 . [36] T. Beyer, D.F. Swinehart, Commun. ACM 16 (1973) 379, https://doi.org/10.1145/362248.362275 . [37] J.L. Bao, J. Zheng, D.G. Truhlar, J. Am. Chem. Soc. 138 (2016) 2690-2704, https://doi.org/10.1021/jacs.5b11938 . [38] A. Gonzàlez-Lafont, T.N. Truong, D.G. Truhlar, J. Chem. Phys. 95 (1991) 8875-8894, https://doi.org/10.1063/1.461221 . [39] J.C. Corchado, E.L. Coitiño, Y.-Y. Chuang, P.L. Fast, D.G. Truhlar, J. Phys. Chem. A 102 (1998) 2424-2438, https://doi.org/10.1021/jp9801267 . [40] Y.-Y. Chuang, J.C. Corchado, D.G. Truhlar, J. Phys. Chem. A 103 (1999) 1140-1149, https://doi.org/10.1021/jp9842493 . [41] C.J. Broyden, Inst. Math. App. 6 (1970) 76–90, https://doi.org/10.1093/imamat/6.1.76 . [42] R. Fletcher, Comput. J. (UK) 13 (1970) 317–322, https://doi.org/10.1093/comjnl/13.3.317 . [43] D. Goldfarb, Math. Comput. 24 (1970) 23–26, https://doi.org/10.1090/S0025-5718-1970-0258249-6 . [44] D.F. Shannon, Math. Comput. 24 (1970) 647–656, https://doi.org/10.1090/S0025-5718-1970-0274029-X . [45] A. Banerjee, N. Adams, J. Simons, R. Shepard, J. Phys. Chem. 89 (1985) 52-57, https://doi.org/10.1021/j100247a015 . [46] J. Baker, J. Comput. Chem. 7 (1986) 385-395, https://doi.org/10.1002/jcc.540070402 . [47] F. Jensen, J. Chem. Phys. 102 (1995) 6706-6718, https://doi.org/10.1063/1.469144 . [48] B.C. Garrett, D.G. Truhlar, J. Chem. Phys. 76 (1982) 1853-1858, https://doi.org/10.1063/1.443157 . [49] X.G. Zhao, D.-h. Lu, Y.-P. Liu, G.C. Lynch, D.G. Truhlar, J. Chem. Phys. 97 (1992) 6369-6383, https://doi.org/10.1063/1.463698 . [50] W.-P. Hu and D.G. Truhlar, J. Am. Chem. Soc. 117 (1995) 10726-10734, https://doi.org/10.1021/ja00148a015 . [51] J. Villà, D.G. Truhlar, Theor. Chem. Acc. 97 (1997) 317-323, https://doi.org/10.1007/s002140050267 . [52] A. González-Lafont, J. Villà, J.M. Lluch, J. Bertrán, R. Steckler, D.G. Truhlar, J. Phys. Chem. A 102 (1998) 3420-3428, https://doi.org/10.1021/jp9807672 . [53] P.L. Fast and D.G. Truhlar, J. Chem. Phys. 109 (1998) 3721-3729, https://doi.org/10.1063/1.476973 . [54] Y.-Y. Chuang, D.G. Truhlar, J. Am. Chem. Soc. 121 (1999) 10157-10167, https://doi.org/10.1021/ja991809i . [55] A. Fernandez-Ramos, D.G. Truhlar, J. Chem. Phys. 114 (2001) 1491-1496, https://doi.org/10.1063/1.1329893 . [56] D.G. Truhlar, D.-h. Lu, S.C. Tucker, X.G. Zhao, A. Gonzàlez-Lafont, T.N. Truong, D. Maurice, Y-P. Liu, G.C. Lynch, Amer. Chem. Soc. Symp. Ser. 502 (1992) 16-36, https://doi.org/10.1021/bk-1992-0502.ch002 . [57] A. Fernandez-Ramos, D.G. Truhlar, J.C. Corchado, J. Espinosa-Garcia, J. Phys. Chem. A 106 (2002) 4957-4960, https://doi.org/10.1021/jp014204t . [58] B.C. Garrett, D.G. Truhlar, J. Chem. Phys. 81 (2001) 309-317 (x), https://doi.org/10.1063/1.447386 . [59] B.C. Garrett, D.G. Truhlar, J. Phys. Chem. A 105 (2001) 2298-2307, https://doi.org/10.1021/jp003783j . [60] J. Pu, D.G. Truhlar, J. Chem. Phys. 116 (2002) 1468-1478, https://doi.org/10.1063/1.1427917 . [61] H.S. Yu, L.J. Fiedler, I.M. Alecu, D.G. Truhlar, Comput. Phys. Commun. 210 (2017) 132-138, https://doi.org/10.1016/j.cpc.2016.09.004 . [62] L.M. Pratt, D.G. Truhlar, C.J. Cramer, S.R. Kass, J.D. Thompson, J.D. Xidos, J. Org. Chem. 72 (2007) 2962-2966, https://doi.org/10.1021/jo062557o . [63] A. Marenich, W. Ding, C.J. Cramer, D.G. Truhlar, J. Phys. Chem. Lett. 3 (2012) 1437-1442, https://doi.org/10.1021/jz300416r . [64] Y.-Y. Chuang, D.G. Truhlar, J. Chem. Phys. 124 (2006) 179903, https://doi.org/10.1063/1.2190219 . [65] T. Yu, J. Zheng, D.G. Truhlar, J. Phys. Chem. A 116 (2012) 297-308, https://doi.org/10.1021/jp209146b . [66] J. Zheng, D.G. Truhlar, Faraday Discuss. 157 (2012) 59-88, https://doi.org/10.1039/C2FD20012K . [67] J.L. Bao, R. Meana-Pañeda, D.G. Truhlar, Chem. Sci. 6 (2015) 5866-5881, . J.L. Bao, D.G. Truhlar, Phys. Chem. Chem. Phys. (2016) . J.L. Bao, X. Zhang, D.G. Truhlar, Phys. Chem. Chem. Phys. (2016) . J.L. Bao, D.G. Truhlar, Phys. Chem. Chem. Phys. .