1        Introduction


1.1              Studying enzymes


Enzymes are very interesting to chemists studying catalysis because as the catalysts used by life they have been honed by millions of years of evolution to efficiently catalyse biological reactions, under biological conditions. Thus these catalysts preform their work efficiently under moderate conditions, with high efficiency and selectivity.

It is desirable to be able to understand the mechanisms of these reactions to design new catalysts, both modified proteins and completely artificial molecules, and also to design inhibitors for those reactions that have pharmacological significance.

The reaction mechanisms of enzymes have therefore been studied using many methods such as kinetic studies, substrate analogues, and isotope effects. As the mechanisms of enzymes have been elucidated it has been found that the enzymatic catalysis involves the interplay of electronic and steric effects to bind the reactants in conformations that facilitate the desired reaction, often distorting the enzymes in the process.

Computational modelling provides a way to access molecules that are not normally observable, such as transition states, and can provide a quantitative estimate of the effect that changes to the molecule might have on the reaction. As enzymes are large molecules only fast (low cost) methods such as molecular mechanics can be used to model them completely. Molecular mechanics is an empirical method, and thus is only as good as the experimental data used to refine the model, and only the broadest predictions may be made about ‘unusual’ structures such as those found in transition states, due to the absence of experimental data. Despite these limitations molecular mechanics, and dynamic simulations based on it (molecular dynamics) have provided useful information on more ‘realistic’ aspects of proteins that are difficult to find with direct methods such as x-ray crystallography, including the protein conformation, and hydration in solution.


1.2              Chemistry of metalloenzymes


During studies of enzyme activity it was discovered that certain enzymes contain metal ions that are essential to their function. The enzymes could be inactivated by removing the ions with a chelating agent, and reactivated by adding ions to the solution. These enzymes were termed metalloenzymes. The modern definition is that metalloenzymes are basically enzymes that contain one or more metal atoms. Usually this definition is restricted to those enzymes in which the metal ion is within the active site, and thus participates in the reaction.

Metalloenzymes are of particular interest because their structure, which may be viewed as a metal centre with a complex ligand, resembles that of ‘conventional’ homogeneous catalysts. In addition almost all enzymes involved in reactions of small (2-5 atom) molecules are metalloenzymes, although not all metalloenzymes are involved in the reactions of small molecules. Because of the industrial importance of reactions involving small molecules these metalloenzymes have received much attention. These reactions include many reactions that require extreme conditions using conventional catalysts, such as the reduction of molecular nitrogen to ammonia, the oxidation of methane to methanol, and the oxidation of ammonium ions to nitrite 0.

The actual chemistry of metalloenzymes is incredibly diverse 0, with many thousands of metalloenzymes being known, ranging from the simple case of a single metal ion with a well defined ligand system, through to systems containing two or three metal ions and enzymes where the metal is contained in a cluster of metal and sulfide ions. This diversity is based on fairly simple chemistry. The metals involved in metalloenzymes are usually light metals, predominantly calcium, magnesium and first row transition metals, the most notable exception to this being molybdenum (present in nitrogenase). These metals are surrounded by amino acid ligands 0; normally these are carboxylate (glutamate or aspartate), sulfide (cysteine, occasionally methionine), or nitrogen (normally histidine) ligands. In addition to amino acids it is common for iron, manganese, and cobalt ions to be chelated by a porphyrin 0.

Figure 1.1 Amino acid ligands, arrows indicate metal bonding atoms


In addition to simple amino acid coordination, some enzymes contain multiple metal ions coordinated to sulfide and sulfur amino acids forming a small cluster. As iron is the predominant metal in the clusters these are usually called iron-sulfur proteins 0.

A            B

Figure 1.2 Other Coordination systems. A) Simple metal porphyrin (M = metal). B) Fe4S4 cluster from aconitase. Three iron atoms are also bound to cysteine, the fourth is bound to water and a substrate analog (methyl isocitrate). Hydrogen was not located in the X-ray structure 1.


This group of metalloenzymes includes significant enzymes such as nitrogenase (cluster also contains molybdenum) and acetyl-CoA decarbonylase/synthase (contains nickel 2).


1.3              Modelling techniques


In order to model molecules with computers assumptions must be made, the various methods of modelling being distinguished by the number and type of these assumptions. It is not necessary to give a detailed description of the theoretical basis of the methods referred to in this document, as that information is available elsewhere 48. Therefore what follows below is simply a brief summary of the nature of the methods, and comments on their particular strengths and weaknesses as regards the modelling of metalloenzymes. Presently, the available methods can be divided into three broad groups, described below.


Molecular mechanics – these empirical methods use a classical treatment, with interatomic interactions described in terms of elastic bonds. Energy is described as a sum of functions describing the contributions of bond stretch, bend, dihedral rotation, van der Waals interactions, and electrostatic factors; the equations and parameters used to calculate the energy are usually called the force field. Variation between force fields occurs both in the interactions considered, and in the functions used to estimate them. Molecular methods are very fast, and can handle large molecules, even whole proteins. The parameters used in the functions are derived by fitting to experimental results. This link to experiment while allowing molecular mechanics methods to give good results for “normal” molecules, also is the major weakness of these methods, as they can only accurately describe molecules containing bonds which have been parameterised in this way. Deriving parameters requires extensive experimental data, which is often hard to obtain. In addition unusual bonding, such as that occurring in transition states can be difficult to describe. A particularly desirable feature of these methods for protein modelling is that a parameter set can be optimised specifically for certain types of molecules, providing greater accuracy for these types of molecules than a more general force field. The AMBER 3 force field is an example of a force field for protein modelling.


 Semi-empirical quantum mechanics – these methods use a quantum mechanical description of the bonding in the molecule, using molecular orbitals. In these methods calculation is simplified by making assumptions, and using empirical parameters to approximate certain aspects of the calculation. As a result semi-empirical calculations can be performed much faster than ab initio calculations. An additional advantage of these methods is that electron correlation effects (see ab initio) are included implicitly as a result of fitting to experimental data. The semi-empirical methods most commonly used for enzyme modelling are the AM1 and PM3 methods. The AM1 Hamiltonian 4-8 was one of the earliest semi-empirical methods with a reasonable number of elements described. The PM3 Hamiltonian 9, 10 was developed in an attempt to produce a systematic technique for semi-empirical parameterisation. The parameters were optimised to an experimental data set with minimal constraints, and as a result PM3 calculations are often more accurate than AM1, although there are certain problems resulting from the original limited data set. A significant weakness of these methods is that they can only describe the atoms that they were parameterised for. In consequence the original methods have often been extended by the addition of parameters for atoms not originally included in the fitting.


Ab initio quantum mechanics – these methods aim to reproduce the behaviour of electrons in the molecule using only basic physics. The aim is to predict the distribution of the electrons in the molecule (wave-function), and the energy of the system. There are two main approaches to this problem in current use – Hartree-Fock theory, and density functional theory. Although both methods use an iterative self-consistent field (SCF) approach to solve the problem, they differ in the mathematical treatment of the system.

Hartree-Fock techniques attempt to solve the Schrödinger equation directly, using the SCF approach to find a wavefunction that is stable in its own electric field. In contrast, density functional theory uses a more abstract approach. It has been mathematically proven 11 that there exists a functional that given only the electron density of a wavefunction can predict its energy. Although this ideal functional is not known, and indeed would probably be impossibly complex, it is possible to use simpler functionals to approximate it. Density functional calculations attempt to find the wavefunction that minimises the energy, equivalent to the solution of the Schrödinger equation.

For practical purposes the main difference between the two methods lies in their treatment of what is termed electron correlation. This effect arises because electrons do not in fact move independently of each other, but tend to stay apart – their motion is correlated. For density functional calculations this effect can be incorporated into the functional, whereas for Hartree-Fock calculations it must be computed separately. Although basic Hartree-Fock calculations are faster than density functional, when an equivalent level of correlation is included the density functional technique becomes faster. Because of the approximate nature of the functional there is no systematic way to improve the accuracy the calculation; for Hartree-Fock the accuracy of the calculation can be improved, albeit at considerable computational cost. For calculations incorporating correlation, unless an extremely high level of accuracy is required, it is usual to use density functional theory.

Common between the two techniques is the use of combinations of simple atom associated functions, termed basis functions, to form the wavefunction. The set of basis functions used to form the wavefunction is called the basis set. These functions are usually chosen to correspond to atomic orbitals, as these form the smallest comprehensive set. As the actual radial functions involved in atomic orbitals are difficult to integrate, they are usually approximated using combinations of Gaussian functions, termed primitives.

As the wavefunction is composed of these basis functions the precision with which the model can reproduce the “real” wavefunction is determined by the “flexibility” of the basis set. Increasing the number and diversity of the basis functions can increase the accuracy of a calculation, but as increasing the number of basis functions also substantially increases the computational cost, the type of basis set used represents a compromise. There are a few strategies for forming basis sets; basis sets consisting only of basis functions corresponding to atomic orbitals are termed minimal basis sets. For instance STO-3G is a minimal basis set in which Slater type orbitals (STO) are approximated by combinations of 3 Gaussian functions (3G). A common strategy for improving the performance of a basis set is to use two basis functions to represent each orbital; these are termed double zeta basis sets. Because the outer, valence orbitals are the most important to bonding it is common to use multiple basis functions only for these orbitals– these are termed split valence basis sets. 3-21G and 6-31G are examples of this type of set (6-31G = 6 primitives for core orbitals, and double valence orbitals composed of 3 and of 1 Gaussian primitives). Finally, in order to provide additional flexibility, basis functions not corresponding to atomic orbitals are often added to a basis set. These functions are classed as polarisation functions, or diffuse functions. Polarisation functions basically resemble orbitals with a higher orbital angular momentum quantum number, for instance d orbitals in a system that formally involves only p orbitals. Polarisation functions allow additional distortion (polarisation) of the electron density around the atom in wavefunction. Diffuse functions are functions that fall off slowly with distance from the nucleus; they allow the wavefunction to expand, as might happen in an anion. The basis set 6-31+G(d) contains both polarisation ((d)), and diffuse functions (+).

In addition to the basic considerations above, the modelling of heavy elements (eg transition metals) introduces additional difficulties. These elements have a large number of electrons, and as the complexity of the calculation rises rapidly with the number of electrons in the system this can seriously slow the calculation. In addition in heavy atoms where the core electrons may be travelling at very high speeds, close to the speed of light, another approximation in the treatments above becomes apparent – the Schrödinger equation is an approximation – it does not include relativistic effects. In order to address these problems it is common to approximate these troublesome, relatively inert core electrons with a fixed potential function – termed a pseudo-potential, or effective core potential (ECP). An example of this type of basis set is the LANL2DZ basis set – this is a double zeta basis set that employs ECPs; for the heavier atoms these ECPs include relativistic effects.

The summary of ab initio techniques above makes no mention of the handling of electron spin in the calculations, as this is a complex subject. Suffice it to say that the most obvious spin effect form a metalloenzyme modelling perspective (apart from the exclusion principle, obviously) occurs in systems that contain unpaired electron spins. In these systems it is necessary to treat electrons of different spins separately, as they are different in energy, using what are called unrestricted methods. The net result of this is that these calculations are much slower, and due to the possibility of two spin states lying close together can be very difficult to model.


1.4              Computational modelling of metalloenzymes


There are a number of difficulties that particularly apply to modelling metalloenzymes; the most obvious is common to almost all enzymes – the large size of the molecule. Although it is now possible to model molecules of this size with molecular mechanics, and even to simulate the thermal movements using molecular dynamics, such modelling still represents a major undertaking in itself. In order to model reactions of large molecules it is usually necessary to use a smaller model system to reduce the size of the system to something that can reasonably be modelled.

As discussed above such molecular mechanics techniques have significant weaknesses in their description of reaction transition states particularly in metalloenzymes as the metals, due to their flexible chemistry, can be difficult to parameterise. Quantum mechanical methods are therefore important in metalloenzyme modelling. Semi-empirical techniques provide a compromise and have been successfully used to investigate many metalloenzymes and are still used as an efficient way to investigate the chemistry of a system, for instance, alcohol dehydrogenase 12, or thermolysin 13.

However the limitations of the semi-empirical treatment mean that its use must usually be validated by comparisons to experiment, or ab initio calculations. An obvious failure of semi-empirical calculations is described in the chloroperoxidase section of this report. The transition states produced by semi-empirical methods may also be problematical as the parameterisation of these methods is usually focused on more conventional chemistries.

Ab initio modelling is therefore often necessary to obtain reliable results. The presence of metal atoms, involving large numbers of electrons in a small volume, makes correlation important to accurately model metalloenzymes. The use of high accuracy, correlated, ab initio techniques for investigating biological transition metal systems, has greatly increased in recent years as computer power has become cheaper, and the structures of these systems have been found. A recent review of this area 14 summarises results obtained for a variety of metalloenzymes, including oxygen activation by cytochrome C oxidase (heme/Cu), the enzyme that reduces oxygen in aerobic respiration; nitrogen activation by nitrogenase (Fe-S-Mo); the reversible oxidation of hydrogen by hydrogenase (CN-Fe-S-Ni); peroxide activation by peroxidases (heme) and other reactions.

Ab initio calculations are very slow, and restrict the model system to no more than a few tens of atoms. In consequence all of the above calculations were performed on very small model systems. This sort of model system cannot reproduce the steric effects of the protein ligand; in order to achieve a combination of the abilities of ab initio techniques with less expensive methods hybrid methods have been developed. In these methods “difficult” parts of the molecule are modelled at the ab initio level, while less critical parts are modelled using another technique, usually molecular mechanics. There are many different techniques used to perform the combination. These methods have been in use for some years, and continue to be used. For example these techniques have been used on the important zinc metalloenzymes such as carbonic anhydrase 15, thermolysin 13 (a protease), alcohol dehydrogenase 12, and dihydrofolate reducatse 16 (not an exhaustive list).

Of the hybrid techniques, the ONIOM (Out Own N-layered molecular Orbital molecular Mechanics) technique 17-22 is particularly attractive as it uses a simple subtractive scheme to provide a general framework in which to combine any methods.

Figure 1.3       ONIOM partitioning and capping, only capping hydrogens shown

For example the ONIOM technique for a two-layer system may be summarised thus:

  1. Atoms in the system are partitioned – the atoms in the high level must be specified.
  2. Two systems are generated, the original, large system, and a small system containing the high-level atoms. Any bonds broken by the reduction are capped by adding atoms, usually hydrogen, to fill the valency. This process is illustrated in Figure 1.3.
  3. Three energies are calculated

a)      The energy of the large system using the low level technique

b)      The energy of the capped small system, using the low level technique

c)      The energy of the capped small system using the high level technique

  1. The ONIOM energy of the system is then the result of the simple subtraction
    ONIOM energy = a – b + c

The real power of the ONIOM comes from being able to calculate first and second derivatives of the energy at the ONIOM level, and thus being able to carry out geometry optimisations on the large system for only moderate cost.

ONIOM is a recent development and, although there is much ONIOM enzyme work in progress, little has been published as yet. The only paper relating the use of ONIOM on a metal containing enzyme that was found described the modelling of catechol O-methyltransferase 23, a magnesium-containing enzyme. The metal was not included in the model system used in that study.


1.5              The project


The intention of this project was to use the ONIOM method to investigate the chemistry of a metalloenzyme catalysed reaction, particularly the role of the protein environment.

Two enzymes were investigated: chloroperoxidase, and methanol dehydrogenase. These are described in chapters 2 and 3.

2        Chloroperoxidase


2.1              Introduction


Peroxidases are a large family of redox-enzymes that catalyse the oxidation of a substrate by hydrogen peroxide. Their chemistry varies slightly, but all of the enzymes contain a metal porphyrin, with an amino acid ligand also coordinated. The reaction is thought to involve intermediates coordinated opposite to the amino acid, completing an octahedral geometry. Most frequently the metal is iron and the amino acid ligand, termed the proximal ligand, is histidine. Variations from this basic chemistry include manganese peroxidases, in which the metal is manganese, and certain special function peroxidases in which the proximal ligand is cysteine, eg chloroperoxidase.  Peroxidases, like oxidases, function primarily as detoxification enzymes, oxidising non-polar toxins, such as pharmaceuticals, to produce more easily excreted polar compounds.

Figure 1.1 Peroxidase redox cycles 24. Ring represents porphyrin. (A and S represent generalised substrates)

The reactions catalysed by peroxidases consist of two steps, the first step being the reaction of the "resting state" of the enzyme with hydrogen peroxide, forming water and an oxidised form of the enzyme, commonly referred to as compound 1; this step has been studied very intensively. The second step, the reaction of compound 1 with the substrate is more varied and less studied. The redox chemistry (Figure 2.1) of iron peroxidases is now accepted to consist of the oxidation of the resting state (iron (III) porphyrin), by hydrogen peroxide to compound 1 (oxo-iron (IV) porphyrin cation-radical 25), followed by the oxidation of the substrate, reducing the enzymes back to the resting state.


2.2              Chemistry and formation of compound 1


As mentioned in the introduction most of the study of peroxidase chemistry has focused on compound 1 and its formation from the enzyme and hydrogen peroxide. There would seem to be two reasons for this, firstly, compound 1 is not only formed in peroxidases, but is also formed from oxygen as an intermediate in the reactions of many oxidase enzymes, such as p450 monooxygenases. The oxidase reaction involves the cleavage of the oxygen double bond 14; studies of peroxidases have therefore focused on comparing the cleavage of the peroxide single bond, to the cleavage of the stronger oxygen double bond. Secondly the first step can be studied with fairly small model systems that include little of the surrounding protein structure, making ab initio modelling feasible.

The presently accepted mechanism for compound 1 formation consists of three steps (Figure 2.2).

Figure 2.2       Formation of compound 1 in peroxidases. B represents a residue acting as a base 25.

The first step is the coordination of hydrogen peroxide to the heme iron in the distal position; density functional calculations 25 indicate that this complex is a doublet, with the spin concentrated on the iron atom. In the next stage of the reaction the proton on the iron coordinated oxygen is transferred to the other oxygen, this step is thought to be catalysed by a specific amino acid residue in the surrounding protein and probably consists of two proton transfers, the residue acting as a proton carrier. The third step is the cleavage of the now activated peroxide bond to form water and compound 1; this is probably a concerted process, occurring simultaneously to the second proton transfer. The importance of the protein environment, i.e. the correct placement of the proton carrier for proton transfer has been investigated with ab intio 25 and molecular dynamics 26 calculations. As a final proof, molecular engineering of myoglobin, which does not naturally show peroxidase activity, has produced artificial peroxidases. By using genetic engineering techniques, the amino acid sequence of myoglobin was altered to place a histidine residue in the appropriate position, and when tested the altered protein showed good peroxidase activity 27.

The electronic structure of compound 1 is more complex than the first intermediate; recent calculations favour a quartet state, in which there is approximately one unpaired electron spin each on oxygen, iron, and the porphyrin ring. A doublet state (unpaired spin mainly on the porphyrin) lies very close in energy, so close in fact that altering the electrostatic environment of the compound, as occurs in a protein, can cause the stable state to change to the doublet 28. In addition in certain proteins the charge/spin on the porphyrin may be delocalised to surrounding amino acid residues 25.


2.3              Epoxidation reaction


There are many oxidation reactions catalysed by peroxidases. From a synthetic perspective the most interesting reactions are those that involve the addition of an oxygen atom to the substrate, including the oxidation of thioethers to sulfoxides, alkanes to alcohols, and alkenes to epoxides. Acting in this way the enzymes are often termed peroxygenases. These reactions are stereoselective; in some epoxidations the enantiomeric excess can exceed 95% 27, 29.

Although other enzymes such as oxidases can catalyse these reactions, peroxidases have more synthetic utility as the concentration of hydrogen peroxide can be easily controlled, providing a way to regulate the reaction.

Of these reactions the epoxidation has been the most studied, as chiral epoxides are very useful synthetic intermediates, and peroxidase oxidation provides a means to prepare epoxides that can be difficult to access by other methods.

Table 2.1 contains a list of some of the epoxidation reactions catalysed by peroxidases.


Table 2.1 Peroxidase catalysed epoxidations


Yield of epoxide

(By product)

%EE (isomer)




Not given

80 (1R)

Engineered myoglobin (L29H/H64L)



40% (aldehyde)




2-heptene (cis)





2-heptene (trans)


Not stated



3-heptene (cis)

40% (allylic hydroxylation)




3-heptene (trans)

No epoxide formed










A number of mechanisms have been proposed for the epoxidation reactions, Figure 2.3 lists possible intermediates 30.

Figure 2.3 Possible intermediates in epoxidation of alkenes by compound 1


 In the reaction of compound 1 analogs of free porphyrins, ESR 31, and kinetic 32 data indicates an intermediate of type E. No published work on the intermediates involved in the enzyme reaction was found, and it is entirely possible that the steric effects of the surrounding enzyme could substantially change the characteristics of the reaction.


2.4              Techniques


The input structures for modelling were retrieved from protein data bank structures for cytochrome C peroxidase 33, chloroperoxidase34 and engineered myoglobin 35,36. In order to reasonably model and manipulate the structures it was necessary to reduce the size of the model system. This was achieved by removing residues that contained no heavy (non-hydrogen) atoms closer than a specified distance to the iron atom; this process was automated using software developed by the author (Appendix 2). For the molecular mechanics docking study this distance was 9Å, the system then being further reduced by manually removing the residues on the proximal face of the heme group (opposite the reactant) leaving only a thiomethanol ligand in chloroperoxidase. In this way the number of atoms in the system was reduced to around 260 (including hydrogen). Hydrogen was added to the structure using the ‘add hydrogen to biomolecule’ function of the Gaussian GUI editor utility Gaussview. SPARTAN  (SGI version 5.03) was used to perform molecular mechanics optimisations, in which the non-hydrogen atoms from the protein were frozen to prevent the system from spreading excessively in the absence of the surrounding protein. This approximation of the protein steric environment was then used for calculations on transition state analogs. For the concerted symmetrical oxygen transfer (mechanism A, Fig 2.3), styrene epoxide was used as a reasonable analog of the transition state. For this structure calculations were done for both R and S configurations of the epoxide. The metallocycle intermediate (fig 2.3B), and an analog of the open chain intermediates (Fig 2.3D, E) were also prepared by changing the bonding of the optimised (S)-epoxide structure. All of these structures were optimised using the minimise option of SPARTAN (molecular mechanics, Merck force field). The energies estimated for the inner parts of the molecule, close to the metal, are expected to be inaccurate. In order to get a comparable index of strain between mechanisms it was attempted to compare only the strain due to the phenyl group. This was calculated using a second optimisation in which the phenyl group is removed, and the oxygen and two carbons of the epoxide group were frozen (the hydrogens remained free). The estimated phenyl group strain was then simply the difference in the strain energies between these two models.


Semi-empirical calculations were performed using the SPARTAN program, specifying the PM3(tm) Hamiltonian. These calculations were performed on a chloroperoxidase model system consisting only of a simple iron porphyrin, with a proximal thiomethanol ligand and distal oxygen.

Hybrid QM/MM calculations were performed using the ONIOM method as implemented by Gaussian 98 37. For the initial ONIOM calculations a small system was used in which only the heme, distal oxygen, and a fragment of the proximal amino acid residue were included. For chloroperoxidase this was the same as the system used for semi-empirical calculations; in addition some calculations were performed on a model system for proximal histidine peroxidases. In this system the proximal ligand was nitrogen-coordinated imidazole.

The ONIOM partitioning scheme used (similar to that used as an example in Fig 1.3) is illustrated (Figure 2.4); this is similar to a scheme that was used in the study of heme/oxygen systems 38.


Figure 2.4       ONIOM partitioning scheme, ab initio atoms highlighted. a) Chloroperoxidase model, atoms with numbers < 20 are ab initio. b) Myoglobin model, atoms with numbers < 24 are ab initio.


Universal force field 39 (UFF) molecular mechanics was used for the low level calculation, as this field has parameters for iron. For ab initio calculations a general basis set was used, consisting of an all electron basis set for iron 40, the 6-31G (d) basis set for nitrogen and oxygen, and STO-3G for carbon and hydrogen. Some calculations were also performed using a similar basis set except the iron used the LANL2DZ basis set, and associated effective core potential (ECP) 41; some calculations were also performed using STO-3G for all high level atoms.


2.5              Calculations

Coordinates calculated for all of the relevant systems are given in Appendix 1


Semi-empirical calculations PM3(tm) calculations on models for compound 1, compound 2, and the ferric resting state of chloroperoxidase all gave structures in which the porphyrin ring is substantially distorted from planar (Figure 2.5). The reason for this distortion is not clear. As the distortion seemed to indicate major problems with the semi-empirical representation of the system, no further calculations at the semi-empirical level were performed.


             Figure 2.5      Two views of the PM3(tm) optimised structure for a compound 1 model of chloroperoxidase (quartet)

Despite extensive attempts it was not possible to obtain SCF convergence for the model systems. Due to this difficulty no ONIOM results were obtained. A summary of the attempts is given in Table 2.2.



No. of SCF iterations

SCF Final energy (Hartree)








Vshift = default (200)






Vshift = default (200)






Vshift = 400






Vshift = 400






Vshift = 600






Vshift = 1000






Quadratic convergence

94 (S) 96(N)





Table 2.2         Attempts at convergence UHF/STO-3G:UFF ONIOM system, all for quartet state. System S is minimal chloroperoxidase model, N is engineered myoglobin minimal system


Simple molecular mechanics calculations were performed on a system prepared as described in section 2.4. The resultant geometries are pictured below (Figure 2.6)



A                                                             B                                             C

D                                                             E                                              F


Figure 2.6       Space filling views of molecular mechanics docking calculations, styrene epoxidation A) reduced 9Å chloroperoxidase system. The void in which the substrate was placed is occupied by three water molecules. B) Coordinated styrene epoxide (S isomer). C) Heavy atoms frozen as in B, styrene phenyl group removed, and hydrogens optimised (steric reference). D) Open chain type intermediate. E) Metallocycle type intermediate. F) Coordinated styrene epoxide (R isomer).


The void in which the phenyl group was placed is most obvious in Figure 2.6C, where the phenyl group has been removed; in the crystal structure this void is occupied by three water molecules (Figure 2.6A). Slight changes in the orientation of the styrene phenyl group are apparent in these figures, mostly as changes in the visibility of atoms eclipsed by the phenyl group. Strain energy differences are given in table 2.3.





Open chain


Phenyl strain energy (kJ/mol)






Table 2.3         Strain due to phenyl group for various intermediates/transition state analogs (see section 2.3 and figure 2.3).


2.6              Conclusions


The bent geometries returned by the semi-empirical optimisations of compound 1 and the related species, are clearly in error as only planar geometries have been observed and predicted for these compounds. This error suggests that the PM3(tm) Hamiltonian used by SPARTAN is inadequate for these compounds.


A great deal of work is required before ab initio methods can be applied successfully to these open shell systems. It is clear from the work performed and published accounts 25,28 that this system can be troublesome in this regard. It may be some years before open shell methods can be routinely applied to these systems.


Without a better estimate of the energy of the active region of the reaction the strain energies are only an approximate guide as to the likelihood of a reaction mechanism. It is notable that the coordinated 1-(S)-epoxide had the lowest strain, followed by the (R)- epoxide. This would seem to indicate that the symmetrical oxygen transfer mechanism (Figure 2.3A) is the best fit in the active site. No reference to the configuration of the enantiomer favoured in the epoxidation of styrene by chloroperoxidase was found, so the significance of the slightly lower energy of the S-isomer is not known. In addition it must be noted that this study used the geometry from an X-ray crystal structure. In solution the conformation of the enzyme may be sufficiently different to change the order of these energies, or to allow the phenyl group to sit in an entirely different position.


Overall, it was found that calculations based solely on semi-empirical methods are unsuitable for modelling this system. Current ab initio methods are unable to be applied. In the absence of these methods information was obtained form molecular mechanics models.


3        Quinone Alcohol Dehydrogenase


3.1              Introduction


Alcohol dehydrogenases are enzymes that catalyse the oxidation of an alcohol to an aldehyde with the simultaneous reduction of a hydrogen carrier compound. Reactions of this type are very important in many biochemical processes, and play a key role in energy metabolism. Alcohol dehydrogenases have synthetic utility as a way to reduce an aldehyde or oxidise an alcohol, both for the gentle conditions of the reaction and more interestingly for the stereoselective nature of many of the reactions. There are many different classes of alcohol dehydrogenase, but most of them fall into two basic groups, depending on the hydrogen carrier either NAD (occasionally NADP) or PQQ (pyrrolquinoline quinone) (Figure 3.1).


Figure 3.1       Redox states of NAD and PQQ


Most eukaryote dehydrogenases belong to the former class, liver alcohol dehydrogenase being a much-studied example 12, 42. The second class consists mainly of prokaryote enzymes of which methanol dehydrogenase 47 is a well-studied example, as it is important for C1 metabolic pathways such as those in methanogens. Another well-studied example is soluble glucose dehydrogenase 46, used in glucose test strips 47. Both classes of alcohol dehydrogenase are metalloenzymes, the first class having a zinc ion in the active site, the second calcium. The mechanism of the first class has been extensively studied, and the currently accepted mechanism consists of hydride transfer from a zinc-coordinated alkoxide to the nicotinamide 12. This reaction has been used as a model for a wide range of biological hydride transfer reactions and has consequently been studied very intensively, and the mechanism has been refined to include effects such as hydride tunnelling 42. The latter group is less studied and the mechanism is less well understood. In these enzymes the calcium is coordinated directly to the PQQ, and to two amino acids 43 (Figure 3.2).


A        B

Figure 3.2 Coordination of Calcium in methanol dehydrogenase. A) Schematic. B) Ball and stick from an X-ray structure 43 showing only those residues coordinated to calcium (dark grey)


Some model studies have been performed with calcium, or another metal, coordinated to PQQ or a related compound, and it has been found that in the presence of a base and calcium ions PQQ can oxidise methanol to formaldehyde 44. In the redox reaction the PQQ orthoquinone is reduced to a 'quinol' state.


3.2              Proposed mechanisms


Survey of the literature suggests that there are three mechanisms favoured for the hydride transfer.

Figure 3.3 Possible mechanisms for methanol oxidation by [CaPQQ]. Mechanisms are illustrated using the larger system (system 2)



Figure 3.4       A) Numbering scheme for PQQ. B) Water bridged O4 mechanism


Two of the mechanisms involve hydride transfer from calcium-coordinated methoxide to PQQ, either to oxygen (O4) forming the quinol directly 45 (Figure 3.3a), or to carbon (C5) to give an intermediate product (Figure 3.3b) that then gives the quinol by tautomerisation 46. In the third mechanism, the reacting species is actually a hemi-acetal, formed by nucleophilic addition of methanol to C5, which then reacts via an addition-elimination mechanism to form formaldehyde, and the quinol 45 (Figure 3.3c). However computational modelling has indicated that this method would have the same transition state as the first mechanism 45, and as the intermediate for this reaction is higher in energy than the coordinated methoxide form, this mechanism is unlikely to be important, and has not been studied in this project. In addition to the two remaining mechanisms above, some calculations were also performed for a water-bridged mechanism (Figure 3.4B). The numbering scheme used when referring to atoms in PQQ is given in Figure 3.4A.


3.3              Techniques


Simple semi-empirical calculations were performed using the SPARTAN program (SGI version 5.03), specifying the PM3 Hamiltonian. Also most of the structures optimised at the ab initio level were pre-optimised at this level.

Ab initio calculations were performed using Gaussian 98 37. A general basis set was used, consisting of the LANL2DZ basis set, and associated effective core potential (ECP), for calcium, and the 6‑31+G(d) basis set for other atoms. Density functional theory was used, the B3LYP functional being used for all calculations, unless otherwise specified. Additionally some calculations were performed using STO-3G basis sets for some atoms.


ONIOM (see 1.3) calculations used the ONIOM feature of Gaussian, the ab initio calculations were performed using the same settings as in the simple ab initio case, and semi-empirical calculations used the AM1 Hamiltonian as implemented by Gaussian 98 37, 4-8.


Figure 3.5       Reactant and product for small system (system 1)

The initial work concentrated on the use of smaller model systems for the polycyclic PQQ system. A reduced system (Figure 3.5, system 1) was modelled to reduce the time needed for the ab initio calculations. This system is the same as the high level partition of the ONIOM calculations, what is referred to as the large system can be seen in the mechanism examples Figure 3.3.

To test the reliability of the Gaussian semi-empirical calculations some purely semi‑empirical calculations were performed using Gaussian; both the AM1 and PM3 37, 9-10 Hamiltonians were used.

A couple of calculations were performed in which some attempt was made to reproduce the effect of the protein environment on the redox chemistry of PQQ. To do this a model system was prepared from an X-ray structure of Methylophilus methylotrophus methanol dehydrogenase 43. A structure was generated containing all residues that came within 3Å of calcium or an atom in PQQ (generated automatically from the protein data bank file, using software designed as part of this project, see 2.4 and appendix 2). This was then manually reduced to leave only the relevant functionalities of changed and polar residues. The non-hydrogen atoms of these residues were then frozen (except PQQ and calcium), and the structure optimised using the PM3 Hamiltonian in SPARTAN.


Coordinates calculated for all of the relevant systems are given in Appendix 1

3.4              Semi-empirical Calculations


A large number of semi-empirical calculations were performed on the smaller model system and on the large PQQ system. The energies obtained for these calculations are contained in Table 3.1



Oxidised PQQ coordinated methoxide

Reduced PQQ (O4), coordinated formaldehyde

Reduced PQQ (C5), coordinated formaldehyde

Transition structure (O4)

Small system



Not calculated


Large system







Transition structure (C5)

Oxidised CaPQQ system

Reduced CaPQQ system (O4)

Reduced CaPQQ system (C5)

Small system

Not calculated



Not calculated

Large system





Table 3.1         PM3 heat of formation of optimised structures (kJ/mol)


It was not possible to locate a transition state for the water-bridged mechanism in these calculations. The most significant aspect of these calculations is that they include the only transition structure found in this study. Using these energies the reaction activation energies may be estimated: 76.2kJ/mol for the large system C5 type structure, and 162.4kJ/mol for the O4 type structure (128.8kJ/mol small system) This may be compared to a published ab initio energy for this (O4) transition state for a similar (large) system of 118kJ/mol 45. The small system energy is probably lower because the O4 transition structure is very distorted, and the smaller system is less rigid.

The geometries obtained are shown in Figure 3.6 & 3.7 (below); the most notable aspects of these geometries are the highly asymmetric coordination environment of calcium in the coordinated methoxide structures, and the position of calcium substantially above the plane of the PQQ. The geometries for the large system are generally more planar.


3.5              Small system ab initio Calculations


The small system ab initio optimisation gave structures (Figure 2.6) with significant differences form those obtained for the PM3 optimisation. These differences were most marked for the coordinated methoxide system. In the ab initio geometries the structures are generally more planar and the coordination around calcium is more symmetrical. Also the geometry of the methoxide oxygen has changed considerably; instead of the bent shape seen in the semi-empirical system, this oxygen is now nearly linear.

Despite many attempts no transition structures were located for this system.



PM3 no ligand                                      PM3 methoxide                                     PM3 O4 transition structure

PM3 formaldehyde                              B3LYP no ligand                                  B3LYP methoxide

 B3LYP formaldehyde

Figure 2.6       Geometries for small system optimisations


3.6              Two level ONIOM semi-empirical / ab initio calculations


The optimised ONIOM (B3LYP: AM1) structures are shown below (Figure 2.7). The ONIOM structure appears quite good, incorporating the geometry changes noted in the progression from semi-empirical to ab initio for the small system, while additionally giving a more planar PQQ geometry in the coordinated methoxide structure.


PM3 no ligand                                      PM3 methoxide                     PM3 O4 transition structure

PM3 C5 transition structure               PM3 formaldehyde (O4 product)       PM3 C5 product

ONIOM no ligand                ONIOM methoxide                               ONIOM formaldehyde (O4 product)

Figure 3.7       Optimised geometries for large system


3.7              Large system ab initio calculations


To check the effectiveness of the ONIOM partitioning scheme in reproducing the chemistry of the system some single point ab initio calculations were performed using the large system at the ONIOM optimised geometry. These results are included in Table 3.2.


Optimised structures


Small PM3

Large PM3

Small DFT


Reaction energy change (O4 mechanism), kJ/mol






Single point B3LYP large system, ONIOM geometry

Basis set (see caption)



Scheme 1

Scheme 2

Reaction energy change (O4 mechanism), kJ/mol





Table 3.2         Reaction energy change for O4 reduction. For the basis sets – normal: the LANL2DZ:6-31+G(d) basis set for all atoms. Scheme 1: normal for atoms that would be in the ONIOM high level, STO-3G for the others. Scheme 2: as Scheme 1, except the STO-3G region now includes the three carbons tertiary carbons at the rear of the orthoquinone ring (1a, 3a, 9a using the labelling scheme of Figure 3.4A).


If the results of the reaction study are compared using only the energy change for the O4 mechanism, as in Table 3.2, it can be seen that the reaction energies are far from consistent between methods. However, if the ONIOM result is ignored, the PM3 energies are fairly close to the ab initio ones; small system 130kJ/mol PM3, vs. 118kJ/mol DFT, and 46kJ/mol PM3 vs. 48kJ/mol for the large system (although ONIOM geometry may be questionable). This agreement is surprising as the geometries produced by the various methods differ greatly.

The density functional energy difference is very different from that predicted by the ONIOM method, suggesting that the ONIOM technique used was not a good representation of the system. As this seemed to implicate the semi-empirical calculations some test calculations using the semi-empirical capabilities of Gaussian were performed. The optimised geometries are shown in Figure 3.8.


Figure 3.8       AM1 (Gaussian 98) optimised structures for coordinated methoxide (a), and coordinated formaldehyde product (b).


These geometries highlight the problem – in the methoxide structure the methoxide has disassociated from the calcium, and in the formaldehyde structure formaldehyde is coordinated through the carbon. Examination of the output files for these calculations reveals that the charge for calcium is always exactly 2; this is also the case for PM3 calculations. This indicates that the calcium has no electrons – this may indicate that it has no orbitals. The error therefore arises from the abnormal charge distribution caused by the absence of these orbitals. This error presumably escapes cancellation in the ONIOM subtraction because the effect is different in the small and large systems.

The distances of the calcium ligands for the various systems are given in Table 3.3.








Methoxide coordinated





Small DFT










Small PM3





Large PM3





Ref 45















Formaldehyde coordinated (O4H)





Small DFT










Small PM3





Large PM3





Ref 45





Table 3.3         Distances (Å), of ligand atoms surrounding calcium, ref 45 = large system rhf/3-21G(d).

These distances show that there are substantial differences between the ab initio and PM3 treatment of the calcium coordination sphere.


3.8              Protein environment


Although the problems with semi-empirical ONIOM calculations, and limited time prevented the use of the ONIOM technique to investigate the effect of the surrounding protein on the reaction some preliminary semi-empirical calculations were performed using SPARTAN. The geometries for these are illustrated in Figure 3.9.


Figure 3.9         Charge environment calculations. A) O4 reduced, B) C5 reduced.

These calculations were only performed for the O4 and C5 reduction products; calculations were attempted for the oxidised form but SCF convergence could not be achieved. The energy difference between the two reduced forms is strongly affected by the environment, causing the most stable reduced form to change from the quinol O4 species (by 21kJ/mol), to the C5 species (by 29kJ/mol). That the less conjugated species may be more stable is surprising, but this stability may provide the explanation to a question that has surrounded PQQ dehydrogenases – in some crystal forms the C5 carbon appears to be tetrahedral (e.g. see Figure 3.10). Many explanations have been suggested for this, such as non-planar tautomers of the semi-quinone 43, 45, as it was assumed that the C5 reduced form would be unstable. If this form is indeed more stable then the tetrahedral geometry is simply due to this compound.


Figure 3.10     Steric environment of PQQ in of Methylophilus methylotrophus methanol dehydrogenase, showing PQQ between tryptophan and the cyclic disulfide. From X-ray structure 43, hydrogens added by Gauss view. Note tetrahedral C5




3.9              Conclusions


The use of the ONIOM method to combine semi-empirical and ab initio calculations in an attempt to reproduce the effects of the large pi system appears to have been at best only partly successful. Although the geometries produced appeared reasonable, the energies, when compared to ab initio results, indicate that the treatment was inaccurate, a problem that was traced to an error in the semi-empirical treatment of calcium. This illustrates that although ONIOM uses a subtractive method to cancel out the effects of the less accurate method in the high level partition, this is not always successful as these errors may influence the system outside the high level region.

Although semi-empirical transition states are not as reliable as ab initio, those obtained indicate that the C5 reaction mechanism should be preferred (76 vs. 162 kJ/mol) for the isolated model system. This of course assumed that the energies are reasonably representative, despite the errors indicated by comparison to the published ab initio results. Until recently this mechanism was not much considered for methanol dehydrogenase, but for steric reasons this mechanism is the only one available in glucose dehydrogenase, and the strong homology of the two enzymes has been used to argue the merits of this mechanism for the other enzyme 46.

As the environment in the protein seems to favour the C5 reduced form of PQQ, it is reasonable to assume that this would facilitate the reaction that forms this compound providing a strong case for the plausibility of the C5 reduction mechanism.

The semi-empirical calculations indicate that the C5 alcohol tautomer may be the most stable form of reduced PQQ in the protein environment provides a plausible explanation for the tetrahedral geometry of PQQ C5 in certain preparations of the enzymes. A more comprehensive study might try to include the effect of non-polar residues as well, as the steric confinement of PQQ by these residues could further alter the energetics of the system.

4        Summary and conclusions


As time is limited it is not always possible to meet the aims set at the inception of the project. It was intended to utilize the ONIOM methodology to provide insights into the role of the protein ligand in metalloenzyme reactions. As a consequence of the difficulties encountered with chloroperoxidase, only limited application of the ONIOM method to methanol dehydrogenase was achieved, and this was of dubious accuracy. Despite this some insight into the effect of the environment was gained using a semi-empirical model.


4.1              Chloroperoxidase


Due to the failure of semi-empirical and ab initio modelling on these systems, no reliable information about the reaction energetics or transition state geometries was obtained. Because of this it was not possible to make any determinations as to the preferred mechanism of the epoxidation reaction. Despite this a limited amount of information about the strain associated with various reaction configurations was obtained from molecular mechanics models.

The modelling of molecules with delocalised radicals is often problematical, and in this case proved impossible in the time available.


4.2              Quinone Alcohol Dehydrogenase


Although the convergence problems that plagued the modelling of the peroxidase system were not encountered with this system, the very limited time remaining prevented a complete investigation of the system. Despite the time limitations it was still possible to perform ONIOM calculations, in which it was attempted to use an ab initio/semi-empirical scheme to reproduce the effect of the large pi system without excessive computational cost. The results of this system seemed to be questionable, and further investigation indicated that there were problems relating to the implementation of the semi-empirical calculation.

Using a different semi-empirical system (SPARTAN) it was possible to obtain transition states for the reaction in a model system, and gain some insight into the effect of a model protein environment. The hydride-to-carbon (C5) mechanism had a much lower transition state energy than the alternative hydride to oxygen mechanism (O4). This is in line with recent reports favoring the C5 mechanism by analogy with glucose dehydrogenase 46. When some aspects of the protein environment were included is was found that the less aromatic C5 reduced tautomer of reduced PQQ was more stable than the aromatic quinol O4 species. This result may provide an explanation for the tetrahedral geometry 43 of C5 observed in some X-ray structures.



  1. For the an introduction to biological metal chemistry try the following refs:
    1. DF Shriver, PW Atkins, CH Langford “Inorganic Chemistry”, 2nd Edn. Oxford University Press, Ch 19 “Bioinorganic chemistry”
    2. KD Karlin, Science, 261, 701-708 (1993)
    3. S Kuby (Ed.) “A Study of Enzymes” CRC Press Vol. II (particularly Ch2) (1991)
  2. H.Lauble,C.D.Stout, Bovine (Bos Taurus) Aconitase  (E.C. Complexed With 2-Alpha-Methyl-Isocitrate, Protein Data Bank structure 1ami, No journal reference given (1994)
  3. SW Ragsdale & CG Riordan, J. Biol. Inorg. Chem. 1(6):489-493 (1996)
  4. W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz Jr., D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell and P. A. Kollman, J. Am. Chem. Soc. 117, 5179 (1995).
  5. M. Dewar and W. Thiel, J. Am. Chem. Soc. 99, 4499 (1977).
  6. L. P. Davis et. al., J. Comp. Chem. 2, 433 (1981).
  7. M. J. S. Dewar, M. L. McKee and H. S. Rzepa, J. Am. Chem. Soc. 100, 3607 (1978).
  8. M. J. S. Dewar, E. G. Zoebisch and E. F. Healy, "AM1: A New General Purpose Quantum Mechanical Molecular Model," J. Am. Chem. Soc. 107, 3902-3909 (1985).
  9. M. J. S. Dewar and C. H. Reynolds, J. Comp. Chem. 2, 140 (1986).
  10. J. J. P. Stewart, J. Comp. Chem. 10, 209 (1989).
  11. J. J. P. Stewart, J. Comp. Chem. 10, 221 (1989).
  12. P. Hohenberg & W. Kohn, Phys. Rev. A, 136, 864-871 (1964)
  13. PK Agarwal, SP Webb, SH Hammes-Schiffer, J. Am. Chem. Soc. 122, 4803‑4812 (2000)
  14. S Antonczak, G Monard, MF Ruis-Lòpez, JL Rivail, J. Am. Chem. Soc. 120, 8825‑8833 (1998)
  15. PEM Siegbahn & MRA Blomberg, Chem. Rev. 100 421-427 (2000)
  16. KM Merz & L Banci, J. Am. Chem. Soc. 119, 863-871 (1997)
  17. R Castillo, J Andrés, V Moliner, J. Am. Chem. Soc. 121 12140-12147 (1999)
  18. T. Matsubara, S. Sieber and K. Morokuma, "A Test of the New "Integrated MO + MM" (IMOMM) Method for the Conformational Energy of Ethane and n-Butane," Journal of Quantum Chemistry 60, 1101 (1996).
  19. I. Komaromi, S. Dapprich, K. S. Byun, K. Morokuma and M. J. Frisch, "A New ONIOM Implementation for the Calculation of Energies, Gradients and Higher Derivatives Using Mechanical and Electronic Embedding II.," Theo. Chem. Act. , in prep. (1998).
  20. S. Dapprich, I. Komaromi, K. S. Byun, K. Morokuma and M. J. Frisch, "A New ONIOM Implementation for the Calculation of Energies, Gradients and Higher Derivatives Using Mechanical and Electronic Embedding I.," Theo. Chem. Act. , in prep. (1998).
  21. M. Svensson, S. Humbel, R. D. J. Froese, T. Matsubara, S. Sieber and K. Morokuma, J. Phys. Chem. 100, 19357 (1996).
  22. S. Humbel, S. Sieber and K. Morokuma, J. Chem. Phys. 105, 1959 (1996).
  23. F. Maseras and K. Morokuma, J. Comp. Chem. 16, 1170 (1995).
  24. K Kahn & TC Bruice, J. Am. Chem. Soc. 122, 46-51 (2000)
  25. S Ozaki, Y Inada, Y Watanabe, J. Am. Chem. Soc. 120, 8020-8025 (1998)
  26. M Wirstam, MRA Blomberg, PEM Siegbahn, J. Am. Chem. Soc. 121 10178‑10185 (1999)
  27. M Filizola & Gh Loew, J. Am. Chem. Soc. 122, 18-25 (2000)
  28. S Ozaki, T Matsui, MP Roach, Y Watanabe, Coordination Chemistry Reviews, 198, 39-59 (2000)
  29. GH Loew & DL Harris, Chem. Rev. 100 407-419 (2000)
  30. A Zaks & DR Dodds, J. Am. Chem. Soc. 117, 10419-10424 (1995)
  31. A Tuynman, JL Spelberg, IM Kooter, HE Schoemaker, R Waver, J. Biol. Chem. 275 3025-3030 (2000)
  32. Z Gross &S Nimri, J. Am. Chem. Soc. 117, 8021-8022 (1995)
  33. JT Groves, Z Gross, MK Stern, Inorg. Chem. 33, 5065-5072 (1994)
  34. D.B.Goodin,D.E.McRee, Protein Databank structure 1cca, Biochemistry, 32, 3319 (1993)
  35. M.Sundaramoorthy,T.L.Poulos, Protein Databank structure 1cpo, Structure (London), 3, 1367, (1995)
  36. T.Matsui, S.-I.Ozaki, Y.Watanabe, E.C.Liong, G.N.Phillips, Protein Databank structure 1ofj, J. Am. Chem. Soc. 118, 9784-9785 (also ref 27) (1996)
  37. T.Matsui, S.-I.Ozaki, Y.Watanabe, E.C.Liong, G.N.Phillips, Protein Databank structure 1ofk, J. Am. Chem. Soc. 119, 6666-6667 (also ref 27) (1997)
  38.  Gaussian 98, Revision A.7, M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, V. G. Zakrzewski, J. A. Montgomery, Jr., R. E. Stratmann, J. C. Burant, S. Dapprich, J. M. Millam, A. D. Daniels, K. N. Kudin, M. C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G. A. Petersson, P. Y. Ayala, Q. Cui, K. Morokuma, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. Cioslowski, J. V. Ortiz, A. G. Baboul, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, J. L. Andres, C. Gonzalez, M. Head-Gordon, E. S. Replogle, and J. A. Pople, Gaussian, Inc., Pittsburgh PA, 1998.


  1. JD Maréchal, G Barea, F Maseras, A Lledós, L Mouaward, D Pérahia, J. Comp. Chem. 21, 282-294 (2000)
  2. A.K. Rappe, C. J. Casewit, K. S. Colwell, W. A. Goddard III and W. M. Skiff, J. Am. Chem. Soc. 114, 10024 (1992).
  3.  CW Bauschlichen, J. Phys. Chem. 97, 3709-3711 (1993)
  4. PJ Hay & WR Wadt, J. Phys. Chem. 82, 299 (1985)
  5. J Rucker & JP Klinman, J. Am. Chem. Soc. 121, 1997-2006 (1999)
  6. Z Xia, Y He, W Dai, SA White, GD Boyd, FS Mathews, Protein Databank structure 1b2n, Biochemistry, 38, 1214‑1220 (1999)
  7. S Itoh, H Kawakami, S Fukuzumi, Biochemistry, 37, 6562-6571 (1998)
  8. YJ Zheng & TC Bruice, Proc. Natl. Acad. Sci. 94, 11881-11886 (1997)
  9. A Oubrie, HJ Rozeboom, KH Kalk, AJJ Olsthoorn, JA Duine, BW Dijkatra, EMBO Journal, 18(19), 5187-5194 (1999)
  10. JA Duine, Journal of Bioscience and Bioengineering, 88(3), 231-236 (1999)
  11.  There are many works that give an introductory view of computational chemistry. For quick reference try:
    Paul van Ragué Schlyer (Ed) “Encyclopedia of Computational Chemistry” John Wiley & Sons (1998)