Abstract
Objectives

In this collaborative project we aimed to link the harmonic quantum transition state theory (hQTST) code of Jonsson and Arnaldsson (1) to the DL-FIND code of Kaestner, (2) for use within ChemShell (3). We first aimed to get the rudimentary interface established, before attempting to improve the existing algorithms and parallelising the code. The goal was to yield a versatile programme package that can calculate accurate rates where tunnelling effects are important, for a variety of interesting problems in astrochemistry, biochemistry and surface science.

With the hQTST approximation rates for chemical reactions where quantum mechanical tunnelling through barriers is the dominant mechanism can be calculated quite accurately (1). The DL-FIND code is distributed as a stand-alone optimiser (2), as well as with the ChemShell package (3), which is a flexible, modular code which is interfaced to a variety of quantum chemistry and molecular mechanics (MM) codes and provides various hybrid, cost-efficient schemes applicable to biochemical, material science, and surface science problems.

Currently the hQTST only works with analytical potential energy surfaces, but in principle could work with any routine that yields forces and energies as a function of coordinates. This is an ideal linkage point for DL-FIND/ChemShell, where the energy/force expression can be handled by any quantum mechanics (QM) code, MM-code or with a combined QM/MM approach, where the reactive centre is treated with accurate QM-methods, while the effect of the surroundings is modeled with a computationally less-demanding MM approach. Therefore, the synthesis of hQTST with DL-FIND will yield a versatile tool to study reaction rates in the tunnelling regime. Tunnelling will likely play a role for activated reactions in the interstellar medium where temperatures are very low, but also even at room temperature for various enzymatic and surface reactions where hydrogen atoms are transferred over short distances.

(1) A. Arnaldsson, PhD thesis, Univ. Washington (2007)

(2) http://ccpforge.cse.rl.ac.uk/projects/dl-find

(3) http://www.cse.scitech.ac.uk/ccg/software/chemshell

Achievements

We have established the linking of the hQTST code with DL-FIND. We have fully tested this link for using the MNDO code for the energy/force evaluation, and applied this to tunnelling in the malonaldehyde tautomerisation reaction with the cheap, semi-emperical OM1 and AM1 Hamiltonians within MNDO and the B3LYP Hamiltonian with the Gaussian Code.

HQTST aims to find a maximum on the minimum action path from reactant to product. This problem is cast in a Feynman path integral approach, where the classical TS is turned into a discretised closed Feynman path. This translates into finding a first-order saddle point on a potential energy surface of dimension (number of images * degrees of freedom of each image). Originally, this search was performed with a Lanczos algorithm, where the lowest eigenmode is sought for every search step, followed by a conjugate gradient line search along the effect force, where the force along the lowest mode is inverted. It was found that with a reasonable guess for the eigenmode (e.g. the unstable mode of the classical TS, or the unstable mode of a previous converged instanton with fewer images and/or at nearby temperature), the mode does not need to be updated so often. If we update the eigenmode only every 20 steps, 20-70% of computational effort can be saved. A further 50% reduction in computational effort was achieved by using the L-BFGS optimiser instead of conjugate gradients.

It transpired that convergence of the reaction rates at temperatures well below the cross-over temperature (where tunnelling becomes by far the dominant mechanism and the quantum TS delocalises more and more) was slow with respect to the number of images of the discretised Feynman path.

We are currently further investigating how to overcome this, as it is impractical to resort to >100 images for similar systems where one needs a more accurate (and expensive) Hamiltonian. Also, the calculation of the pre-factor of the reaction rates becomes impractical for this large number of images. Since the instanton path itself seems to converge already with fewer images, we are working the implementation of a suitable interpolation scheme for the instanton path as well as calculating the pre-factor from a more modest number of images. First test cases suggest that should be possible, but more extensive testing is necessary.

We have not yet tested the code with a QM/MM example or with different QM-codes, because we did not have any other QM-codes available during this project, but we will soon test this on other compute systems available to us. In the future, we also intend to impement and test the parallelisation of the code by distributing the images among the available nodes. The code can be further parallelised by the intrinsic parallelisation of the QM- and MM-codes that are used.