Introduction

The ChemShell computational chemistry environment [1,2] is used to perform combined quantum mechanical (QM) and molecular mechanical (MM) calculations. In a QM/MM calculation, accurate (but computationally expensive) QM methods are used to treat the reactive site of a chemical system while a cheaper classical force field is used to describe the surrounding environment, where a description of the electronic structure is not required. In this way the overall computational cost is reduced compared to a full QM calculation without sacrificing overall accuracy.

The QM/MM approach is particularly useful in the study of heterogeneous catalysis. ChemShell supports embedded cluster calculations [2,3] in which the bulk material or surface is represented by a finite MM cluster model, optionally with additional point charges to mimic the effect of bulk electrostatics. The active site is modelled by a high-level QM calculation with electrostatic embedding, where the cluster environment is represented by point charges in the QM Hamiltonian.

On HECToR the GAMESS-UK code [4,5] is typically used for the QM calculations and GULP for the MM calculations [6], both of which are linked in directly to ChemShell as libraries for high parallel performance. The DL-FIND geometry optimisation library within ChemShell [7,8] can be used to characterise reactions taking place on the surface by optimising reactants, products, reaction paths and transition states.

The aim of this distributed Computational Science and Engineering (dCSE) project is to improve the performance of ChemShell for large-scale geometry optimisation by implementing microiterative techniques in DL-FIND. The computational time of a QM/MM optimisation is usually dominated by the quantum mechanical component (100s of atoms), even though the surrounding MM region is typically much larger (1000s of atoms). In a standard geometry optimisation each step requires a QM and MM evaluation at the new geometry. This means that the geometries of the QM and MM regions have to relax at the same rate. In a microiterative optimisation scheme [9], the system is divided into an inner region containing (at a minimum) the QM atoms, and an outer region containing the rest of the system. After each optimisation step of the inner region (the `macroiterative' cycle), the outer region is fully optimised (the `microiterative' cycle). By optimising in this way, the number of QM evaluations is reduced significantly at the cost of increasing the number of MM evaluations of the outer region. As MM evaluations are usually much cheaper, this reduces the overall computational time.

The microiterative QM/MM minimisation scheme implemented in DL-FIND follows an earlier implementation in the HDLCOpt module in ChemShell [9]. HDLCOpt was designed specifically for biological systems and is not suitable for materials chemistry. It has been largely superseded by DL-FIND, but until this project DL-FIND lacked microiterative functionality. This project was therefore both a refactoring exercise (transferring functionality from the HDLCOpt module to DL-FIND, with the intention of deprecating HDLCOpt entirely) and an optimisation of the code (speeding-up DL-FIND calculations for a range of algorithms). The DL-FIND implementation also goes beyond the HDLCOpt version by extending the microiterative approach to a wider range of optimisation methods.

9 months of development effort were deployed on the project over a 12 month period from April 2012 to March 2013. The code development objectives were:

• Implementation of the microiterative optimisation framework in DL-FIND, involving the restructuring of the code to introduce macroiterative and microiterative cycles.
• Modification of the ChemShell/DL-FIND interface to support microiterative optimisation.
• Microiterative energy minimisation, in which the outer region is minimised after each step taken in the inner region.
• Microiterative transition state optimisation using the P-RFO and dimer methods, in which the inner region is optimised to a saddle point while the surrounding environment is relaxed to a minimum.
• Microiterative reaction path optimisation, where the nudged elastic band method (involving multiple geometries along the path) is applied to the inner region, while the environment of each image is relaxed.

All methods have been successfully implemented and tested on HECToR Phase 3 for the case of hydrogen dissociation over Li-doped MgO, as discussed below.