Introduction

The MD code CHARMM [1] is one of the most popular packages for biological simulations. Over recent years a number of interfaces to quantum mechanical programs have been developed, initially based on semi-empirical wavefunction and more recently including interfaces to the GAMESS-US package[2], and CADPAC. The GAMESS-UK interface follows a similar approach to these earlier QM interfaces. It provides a way of performing dynamics and minimisation calculations using a QM/MM Hamiltonian, in which the QM calculation is performed on the active site of the system. The GAMESS-UK interface code is now part of the academic CHARMM distribution, and those interested should contact the Karplus laboratory for licencing.

New Developments

QM/MM Interface Approaches

GAMESS-UK incorporates a DFT module in which an auxiliary basis fit of the charge density is used to provide an approximation to the Coulomb energy. We have used these elements of GAMESS-UK to implement a model in which the charge density of the classical system is included in the QM Hamiltonian not as a set of point charges (as is conventional for QM/MM calculations) but as a continuous charge distribution represented as a sum of Gaussian terms. This allows greater overlap between the QM and MM charge distributions without the introduction of major artefacts and thereby permits the exploration of a number of QM/MM schemes. An extensive comparison of proton affinities and rotational barriers using various QM/MM techniques including those incorporating the Gaussian blur and the associated double link atom schemes has been completed [3].

Modelling Reaction Pathways - The Replica Path Method

The replica path method involves the simultaneous optimisation of a series of geometries of the reacting system, corresponding to a series of points along the reaction pathway. The target function for the combined minimisation comprises the sum of the configurational energies, together with a series of penalty functions which ensure that the structures represent the reaction path. In the CHARMM implementation the important penalty terms include a term ensuring that the points are equally spaced in coordinate space, computed by taking the RMSD for the cartesian coordinates of successive points. Using the best fit (RMS) orientation the distance between each pair of adjacent structures is compared with the mean of all such distances using a harmonic penalty function.

A similar term, related to the angles in coordinate space (computed from the distances between successive groups of three points)

provides a penalty term when significant deviations occur from linearity, thus avoiding paths that double back on themselves. COSMAX is a parameter which determines the largest angle for which no constraint is applied (set to 0.975 in this study); the term E_angle is set to zero for cos(θ) ≤ COSMAX.

The non-equilibrium nature of the points on the pathway can be used to integrate the energetics of the reaction using a PMF formalism:

Here, RMS(i,j) is the weighted best-fit coordinates of structure i fitted onto structure j as a reference, and ∇E(i) is the negative force on the structure at point i, excluding the constraint forces. In future it is expected that the PMF approach may be used in conjunction with an MD simulation of the replicated system to access free energies of activation.

The QM/MM Replica Path Implementation

The system is divided into three regions. Part of the system may be unreplicated, so that all points on the pathway can share a common environment which is optimised to a compromise configuration. If used with care, this helps to ensure that the study of the reaction is not complicated by changes in the conformation of spectator groups, or solvent rearrangements. Within the replicated region, a subset of atoms may be specified as QM.

The parallelisation of the algorithm has been performed, using an approach that maps the QM component of the each reaction path points onto a different subset of processors, as shown below.

The implementation has been tested in the Chorismate Mutase system and the results are in press[4].

Parallelisation

GAMESS-UK has been parallelised using both MPI and the Global Array tools (see http://www.emsl.pnl.gov/docs/global). The extent of parallelisation in greater in the GA tools version, and better use is made of the global memory of the parallel machine. Both versions can be coupled with CHARMM. Over the past year the Global Array version has been ported to a number of new platforms, including the Compaq AlphaServer SC at Pittsburgh and SGI Origin 3800 at Sara.

References

[1] CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations, J. Comp. Chem. 4, 187-217 (1983), by B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and M. Karplus. http://yuri.harvard.edu/

[2] K. P. Eurenius, D. C. Chatfield, B. R. Brooks, and M. Hodoscek, Int. J. Quant. Chem. 60, 1189 (1996)

[3] Optimization of Quantum Mechanical/Molecular Mechanical Partitioning Schemes: Gaussian Delocalization of MM Charges and the Double Link Atom Method, D. Das, K.P. Eurenius, E.M. Billings, P. Sherwood, D.C. Chatfield, M.Hodošček, and B.R. Brooks J. Chem. Phys 117 (2002): 10534-10547.

[4] Exploring the quantum mechanical/molecular mechanical replica path method: a pathway optimization of the chorismate to prephenate Claisen rearrangement catalyzed by chorismate mutase. H.L. Woodcock, M. Hodoscek, P. Sherwood, Y.S. Lee, H.F. Schaefer and B.R. Brooks, Theor. Chem. Acc (2003) in press.