Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Chem Biol Interact. Author manuscript; available in PMC 2011 September 6.
Published in final edited form as:
PMCID: PMC2894988

Reaction Profiles of the Interaction between Sarin and Acetylcholinesterase and the S203C Mutant: Model Nucleophiles and QM/MM Potential Energy Surfaces


The phosphonylation mechanism of AChE and the S203C mutation by sarin (GB) is evaluated using two reaction schemes: a small model nucleophile (ethoxide, CH3CH2O) and quantum mechanical/molecular mechanical (QM/MM) simulations. Calculations utilizing small model nucleophiles indicate that the reaction barrier for addition to GB is the rate-limiting step for both ethoxide and ethyl thiolate (CH3CH2S); moreover, the activation barrier for addition to the phosphorus center of GB by ethyl thiolate is significantly larger (13.2 kcal/mol) than for ethoxide (8.3 kcal/mol). The decomposition transition state for both nucleophiles was determined to be ~1 kcal/mol. QM/MM simulations for AChE suggest a similar reaction mechanism for phosphonylation of the catalytic S203; however, the relative energetics are altered significantly compared to the isolated system. QM/MM results indicate that formation of the penta-coordinate intermediate is the rate–limiting step in the enzymatic system, with an activation barrier of 3.6 kcal/mol. Hydrogen-bonding interactions between the fluoride leaving group of GB with Y124 in AChE are observed throughout the reaction profile. The S203C mutation alters the relative energetics of the reaction, increasing the energy barrier for formation of the penta-coordinate intermediate to a value of 4.7 kcal/mol; moreover, the penta-coordinate intermediate (as product) is stabilized by an additional 6 kcal/mol when compared to wild-type AChE.

Keywords: Nerve Agent, Acetylcholinesterase, Phosphonylation, QM/MM


Nerve agents, a subfamily of organophosphorus compounds (OPs), have a history of utilization as chemical weapons and insecticides. There are two main families of chemical weapons, described as the G-series, and V-series nerve agents. The main structural difference between these two families is that G-series nerve agents contain a small fluoride leaving group, while the V-series nerve agents possess a larger amino-thiolate leaving group. These compounds exert their toxicity via the inhibition of acetylcholinesterase (AChE), leading to a buildup of acetylcholine and resulting in overstimulation of cholinergic receptors. The active site of acetylcholinesterase is the typical serine hydrolase catalytic triad, comprised of Ser203, His447, and Glu334 (numbering for human AChE).1 OPs inhibit AChE through phosphylation of the catalytic serine and the formation of a stable OP-AChE adduct (Figure 1), thereby resulting in a P–O bond that is resistant to cleavage by water molecules in the gorge. Although enzymatic function can be regenerated by cleaving the P–O(Ser) bond using a stronger nucleophile, most commonly an oxime,2 reactivation of the native serine is complicated by inefficient reactivation by oximes for specific OPs; frequently, the efficacy of an oxime is correlated to a specific OP nerve agent. In addition, the OP adduct can undergo dealkylation (aging) which renders the enzyme unreactivatable by common methods.3

Figure 1
Proposed Mechanism for Inhibition of AChE by OPs.

Due to the difficulty in regenerating cholinesterase activity, a detailed study of the enzymolysis of OPs was undertaken previously with the hopes of creating a mutant enzyme capable of catalyzing the hydrolysis of nerve agents.4,5,6 One of the more interesting mutations, a substitution of the catalytic serine for cysteine (S203C) was found to render the protein inactive against CMP,7 an OP with a coumarin leaving group. A chemically similar mutation (Ser/Cys to selenocysteine) in other enzymes has resulted in increases of 100–500× activity over the corresponding wild-type enzyme.8 It is, therefore, not clear why the substitution of a stronger nucleophile in the form of cysteine resulted in an inactive enzyme for the S203C mutant.

The hydrolysis of OPs has a long history of study, both experimentally and computationally, that suggests the hydrolysis of nerve agents proceeds via an addition-elimination pathway, forming a stable penta-coordinated phosphorus intermediate.9,10,11,12,13 Along this pathway, there are two main transition state structures: the first corresponds to the association between nucleophile and OP, which forms a pentacoordinate intermediate, and the second transition state corresponds to cleavage of the OP-leaving group bond, followed by dissociation of the complex. Depending on the OP, there may also be additional transition states that involve reorientations of the O-alkyl substituent. Predicted solvent effects on the potential energy surface vary between studies, with some studies suggesting solvation increases the association barrier while decreasing the decomposition barrier,9 and another predicting solvation will decrease both barriers.10 The majority of these studies have not optimized the OP in an implicit polarizable continuum model (PCM) solvent representation; instead, they have computed single-point energy evaluations with implicit solvation using the gas-phase geometries, which may be the cause of this disagreement. However, the majority of computational studies indicate that the association between nucleophile and OP is the rate-limiting step, and have accurately reproduced the observed experimental ~9 kcal/mol energy barrier of GB.14

Several studies have expanded beyond aqueous hydrolysis in order to study the reaction at the catalytic site, using small fragments of the active site,15 as well as the full enzymatic environment through quantum mechanical, molecular mechanical (QM/MM) calculations.16,17,18 These studies predict the presence of a stable penta-coordinate intermediate, similar to the isolated systems. However, studies that utilize a full enzymatic approach suggest that, in the active site, decomposition of the intermediate is the rate-limiting step, at least for the acylation/deacylation of AChE.16

This study aims to determine the potential energy surface for phosphonylation of the AChE active site by GB, as well as to clarify the roles that solvation and active-site residues play in this reaction. Results will be compared to the potential energy surface for aqueous hydrolysis of GB by both ethoxide (CH3CH2O) and ethyl thiolate (CH3CH2S). In addition to wild-type (wt) AChE, the previously studied S203C mutant7 is investigated in order to elucidate the cause of the experimentally observed inactivity.

Computational Methods

Calculations were carried out using ethoxide as a model nucleophile. Calculations were completed with the Gaussian 0319 suite of programs for the gas-phase optimizations, and in Gaussian 0920 for the condensed phase optimizations. The B3LYP hybrid density functional21 and a 6-31+G* basis set22 were utilized for all geometry optimizations, as double zeta basis sets have been shown to be sufficient for geometry optimizations along the OP hydrolysis reaction coordinate.12 For condensed-phase optimizations, Tomasi’s polarizable continuum model23 was utilized, with atomic radii taken from Cramer’s and Truhlar’s SMD solvation method.24 Transition-state geometries were optimized, and vibrational frequency calculations confirmed each structure was a single-order saddle point on the potential energy surface. Each transition state was propagated by 10% along the imaginary vibrational frequency in both directions, and the resulting geometries were minimized to obtain the local minima connected by the transition state. Single-point energy evaluations were carried out on relevant structures using the B3LYP method with a larger 6-311+G** basis set for more accurate energetics. Energies presented herein will be ΔG(298K) with thermal and entropic corrections to the free energy as obtained from the B3LYP/6-31+G* vibrational frequency calculations.

For the QM/MM simulations, the structure for human acetylcholinesterase was constructed using the crystal structure 1B41.25 Missing loops in the protein were modeled based on homology to electric eel AChE.26 Missing side-chain atoms were replaced using the xLEaP module of AMBER,27 and the structure was then minimized. The two crystallographic waters bridging Glu202 were retained in the structure. While the role of these two waters appears to be structural, in that they stabilize the orientation of the Glu202 sidechain, the two water molecules may play an indirect role in the potential energy surface. Glu202 has previously been shown to greatly destabilize the acylation transition state for acetylcholine hydrolysis by 7.4 kcal/mol through electrostatic interactions.28 To corroborate this possibility, initial testing revealed that the relative energetics are largely unchanged between treating the two waters using QM and MM treatments, suggesting that their roles in the studied potential energy surface are mainly structural, and therefore the MM treatment is suitable. Protonation states of titratable residues were determined for a pH of 7.0 using the program PDB2PQR.29 The orientation of the GB–AChE intermediate was added to the catalytic site based on the orientation of GB in the GB-AChE product from crystal structure 2JGG,30 and by comparison to our model system’s geometries.

QM/MM optimizations were then carried out using ChemShell31 as an interface to Turbomole (v5.10)32 for the QM calculations, and ChemShell’s internal version of DL_POLY33 for the MM treatment using the CHARMM force field.34 Relaxed potential energy surface scans were carried out by varying the reaction coordinate, defined as the difference between the P–O(Ser) and P–F bonds, and optimizing all other degrees of freedom at each step of the fixed reaction coordinate. The QM layer was selected to include Glu202, Ser203, His447, Glu334, and the OP, and was treated using RI-BLYP35 with an SV(P) basis set, which is comparable to the 6-31+G* basis set used in calculations for the isolated system.36 All atoms within 15 Å around the GB–Ser203 adduct were allowed to optimize in the MM treatment. The rest of the protein was held rigid. An electrostatic embedding scheme coupled the MM layer to the QM system, allowing the QM system to polarize due to electrostatic influences of the protein environment. Single-point energy calculations were carried out on each point using BLYP with a TZVP basis set. An identical method was utilized to determine the potential energy surface of the S203C mutant, after replacing the Ser203 with a cysteine residue. Due to the inherent difficulty of calculating vibrational frequencies of an enzyme, the QM/MM energetics will be presented as ΔE values in units of kcal/mol.

Results and Discussion

Model systems

A three-step reaction profile was determined for the hydrolysis of GB by ethoxide, with transition states corresponding to the association of ethoxide with GB, a conformational rearrangement of the O-alkyl group of GB, and dissociation of the fluoride leaving group. This reaction profile is consistent with the association-elimination pathway described previously in the literature (Figure 2). The gas-phase barrier for association of ethoxide with GB was determined to be 5.5 kcal/mol, while the aqueous (PCM) barrier increased to 8.3 kcal/mol. The P–O bond distance in the aqueous system is 2.5 Å. A small barrier for rearrangement of the O-alkyl substituent was observed in both phases: 2.7 kcal/mol in the gas phase and 1.3 kcal/mol in aqueous (PCM) solution. The barrier for loss of fluoride was comparable to O-alkyl rearrangement with barriers of 2.3 kcal/mol in gas phase, and 1.0 kcal/mol in solution, occurring at a bond distance of 2.1 Å in the PCM calculations. These calculations, carried out with full condensed-phase (PCM) optimizations and frequency calculations, determine conclusively that solvent does, in fact, increase the barrier for formation of the penta-coordinate intermediate with GB. While the barrier for decomposition of the intermediate is decreased slightly, there is only a slight change when moving from gas phase calculations to PCM solvation. For ethoxide with GB, formation of the penta-coordinate intermediate is the rate-limiting step, and the calculated 8.3 kcal/mol barrier is comparable to the observed barrier of ~9 kcal/mol for aqueous hydrolysis.

Figure 2
Potential energy surface (ΔG(298K), kcal/mol, B3LYP/6311+G**//B3LYP/6-31+G*) for the hydrolysis of GB by ethoxide in gas phase (open circles) and aqueous (solid triangles) solution (PCM).

The potential energy surface for ethyl thiolate as the nucleophile bears a resemblance to hydrolysis of GB by ethoxide, although the relative energetics are altered significantly (Figure 3). The barrier for association between ethyl thiolate and GB increases to 13.2 kcal/mol in solution, and the transition state occurs at a P–S bond distance of 2.7 Å. Almost no stabilization is observed moving from the association transition state to the first stable penta-coordinate intermediate, as opposed to the roughly 10 kcal/mol stabilization for ethoxide. A 1.0 kcal/mol barrier for O-alkyl rearrangement was observed. As is the case for hydrolysis by ethoxide, decomposition of the intermediate occurs with a barrier of only 0.9 kcal/mol. These results suggest that aqueous hydrolysis utilizing ethyl thiolate is less favorable than aqueous hydrolysis by ethoxide.

Figure 3
Potential energy surface (ΔG(298K), kcal/mol, B3LYP/6311+G**//B3LYP/6-31+G*) for the aqueous (PCM) hydrolysis of GB by ethyl thiolate. *The PC structure is omitted due to failures in convergence of the PCM-optimization.

QM/MM Treatment

The reaction profile for phosphonylation of the AChE catalytic serine by GB follows a similar trend as the isolated system. The most notable change, however, is the elimination of the O-alkyl group rotation from the reaction profile. At the OP-AChE penta-coordinate intermediate geometry, the protein environment, specifically Glu202 and His447, prevent GB from interacting with the catalytic serine if the O-isopropyl group is not oriented toward the fluorine of GB. Therefore, the QM/MM-derived potential energy surface contains only two transition states, corresponding to formation and decomposition of the pentacoordinate intermediate.

In the QM/MM computations with AChE, the barrier for association of the S203 nucleophile with GB was calculated to be 3.6 kcal/mol (Figure 4), and the transition state was located at a P–O distance of 2.7 Å, comparable to the geometry determined for aqueous hydrolysis (Figure 5a). It should be noted that our QM/MM simulations predict the migration of the hydroxyl hydrogen of S203 to His447 is concerted with the approach of GB. A barrier of 10.5 kcal/mol for decomposition of the penta-coordinate intermediate and loss of the fluoride leaving group was observed, with the maximum being at a P–F distance of 1.8 Å. Although the P–F bond length is significantly shorter than the 2.1 Å determined for aqueous hydrolysis with the ethoxide model system, the observed hydrogen bonding between the fluoride leaving group and the Y124 residue in AChE is most likely the cause of the geometric perturbation (Figure 5a). While the decomposition barrier is significantly larger than the values obtained for the model system, these results are consistent with reported data for the GA-AChE complex studied previously using a similar methodology.18 For S203 as the nucleophile, the rate-limiting step, however, remains formation of the penta-coordinate intermediate.

Figure 4
QM/MM (ΔE, kcal/mol, RI-BLYP/TZVP//RI-BLYP/SV(P)) potential energy surfaces for phosphonylation of AChE (solid squares) and the S203C AChE mutant (open circles) by GB.
Figure 5
QM/MM-optimized geometries of GB in the active site of AChE. (A) The transition state for addtion of GB to wt AChE. (B) The overlap of the GB intermediate with wt AChE (elemental coloring) and the S203C mutant (green) showing the shift of intermediate ...

The S203C mutation causes two perturbations to the potential energy surface relative to wt AChE. The barrier for association of GB with S203C is increased slightly to 4.5 kcal/mol. However, there is a much larger stabilization of the penta-coordinate intermediate (Figure 4), even though aqueous hydrolysis suggests the there would be very little stabilization of the intermediate. This perturbation of the potential energy surface is most likely the result of stabilization of the intermediate geometry in the active site of AChE due to the surrounding amino acid environment, including the oxyanion hole. The aqueous hydrolysis results suggest that the nucleophile is decoupled from the intermediate decomposition, as is evident by the similarity in energetic barriers for decomposition for the two nucleophiles. Furthermore, due to this stabilization, the effects of the S203C mutation should only be on the formation of the intermediate, and loss of the cysteine moiety as a leaving group would only facilitate decomposition. Thus, it appears that the S203C mutant might be viable for OP turnover, and perhaps proper folding of the tertiary structure or binding of CMP, and not activity, was to blame for previous failed attempts7 with the S203C mutant of AChE.


The reaction profiles for the phosphonylation of AChE and the S203C mutant by the nerve agent GB were evaluated using an isolated model system as well as with QM/MM calculations involving the full enzymatic system. In the isolated system, the change from an oxygen to a sulfur nucleophile results in a 4.9 kcal/mol increase in the barrier for formation of the penta-coordinate intermediate; however, the same modification in the active site of AChE only results in a 1 kcal/mol increase. This suggests that the main effects of altering the nucleophile are isolated to the formation transition state and intermediate structures. However, substituting a sulfur nucleophile for the common oxygen nucleophile does not change the overall scheme of the reaction: an addition-elimination pathway. While the decomposition of the intermediate is predicted to encounter only a small barrier during aqueous hydrolysis, QM/MM simulations predict a much larger (10.5 kcal/mol) barrier. The S203C mutation is predicted to increase the barrier for intermediate formation to 4.5 kcal/mol, as well as making intermediate formation reaction more exothermic than for wt AChE; however, formation of the penta-coordinate intermediate remains the rate-limiting step. Our QM/MM simulations also suggest that the proton on the His447 δ–nitrogen does not migrate to Glu334 during the hydrolysis mechanism, resulting in a positively charged imidazolium species. In addition to the previously known interactions with the oxyanion hole, it appears that the sidechain of Y124 can take part in hydrogen bonding with the fluoride leaving group of GB. These calculations suggest that the S203C mutant of AChE should be active against OPs, although the association between GB and S203C–AChE is predicted to be slower based on energetic barriers.


This research was supported by grants from the NIH (U54-NS058183) and the U.S. Army (W91ZLK-06-C-0012). The authors are grateful to the Ohio Supercomputer Center for a generous grant of computational resources.


Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


1. Sussman JL, Harel M, Frolow F, Oefner C, Goldman A, Toker L, Silman I. Atomic Structure of Acetylcholinesterase from Torpedo californica: A Prototypic Acetylcholine-Binding Protein. Science. 1991;253:872–879. [PubMed]
2. Worek F, Szinicz L, Eyer P, Thiermann H. Evalulation of oxime efficacy in nerve agent poisoning: Development of a kinetic-based dynamic model. Toxicol Appl Pharmacol. 2005;209:193–202. [PubMed]
3. Michel HO, Hackley BE, Jr, Berkowitz L, List G, Hackley EB, Gilliam W, Paukan M. Aging and dealkylation of soman (pinacolylmethyl- phosphonofluoridate)- inactivated eel cholinesterase. Arch Biochem Biophys. 1967;121:29–34. [PubMed]
4. Oppenoorth FJ, van Asperen CJ. Allelic genes in the housefly modified enzymes that cause organophosphate resistance. Science. 1960;29:298–299. [PubMed]
5. Newcomb RD, Campbell PM, Ollis DL, Cheah E, Russell RJ, Oakeshott JG. A single amino acid substitution converts a carboxylesterase to an organophosphorus hydrolase and confers insecticide resistance on a blowfly. PNAS. 1997;94:7464–7468. [PubMed]
6. Masson P, Nachon F, Broomfield CA, Lenz DE, Verdier L, Schopfer LM, Lockridge O. A collaborative endeavor to design cholinesterase-based catalytic scavengers against toxic organophosphorus esters. Chem Biol Interact. 2008;175:273–280. [PubMed]
7. Shafferman A, Kronman C, Flashner Y, Leitner M, Grosfeld H, Ordentlich A, Gozes Y, Cohen S, Ariel N, Barak D. Mutagenesis of human acetylcholinesterase. Identification of residues involved in catalytic activity and polypeptide folding. J Biol Chem. 1992;267:17640–17648. [PubMed]
8. Johansson L, Gafvelin G, Arnér ESJ. Selenocysteine in proteins-properties and biotechnical use. Biochimica et Biophysica Acta. 2005;1726:1–13. [PubMed]
9. Zheng F, Zhan CG, Ornstein RL. Theoretical studies of reaction pathways and energy barriers for alkaline hydrolysis of phosphotriesterase substrates paraoxon and related toxic phosphofluoridate nerve agents. Perkin Trans. 2001;2:2355–2363.
10. Wang J, Gu J, Leszczynski J. Phosphonylation Mechanisms of Sarin and Acetylcholinesterase: A Model DFT Study. J Phys Chem B. 2006;110:7567–7573. [PubMed]
11. Florián J, Warshel A. Phosphate Ester Hydrolysis in Aqueous Solution: Associative versus Dissociative Mechanisms. J Phys Chem B. 1998;102:719–734.
12. Seckute J, Menke JL, Emnett RJ, Paterson EV, Cramer CJ. Ab Initio Molecular Orbital and Density Functional Studies on the Solvolysis of Sarin and O,S-Dimethyl Methylphosphonothiolate, a VX-like Compound. J Org Chem. 2005;70:8649–8660. [PubMed]
13. Åqvist J, Kolmodin K, Florian J, Warshel A. Mechanistic alternatives in phosphate monoester hydrolysis: what conclusions can be drawn from available experimental data? Chem & Biol. 1999;6:R71–R80. [PubMed]
14. Larsson L. The Alkaline Hydrolysis of isoPropoxy-methyl-phosphoryl Fluoride (Sarin) and some Analogues. Acta Chem Scand. 1957;11:1131–1142.
15. Majumdar D, Roszak S, Lszczynski J. Probing the Acetylcholinesterase Inhibition of Sarin: A Comparative Interaction Study of the Inhibitor and Acetylcholine with a Model Enzyme Cavity. J Phys Chem B. 2006;110:13597–13607. [PubMed]
16. Vasilyev VV. Tetrahedral intermediate formation in the acylation step of acetylcholinesterases. A combined quantum chemical and molecular mechanical model. J Mol Struct (Theochem) 1994;304:129–141.
17. Hurley MM, Wright JB, Lushington GH, White WE. Quantum mechanics and mixed quantum mechanics/molecular mechanics of model nerve agents with acetylcholinesterase. Theor Chem Acc. 2003;109:160–168.
18. Kwasnieski O, Verdier L, Malacria M, Derat E. Fixation of the Two Tabun Isomers in Acetylcholinesterase: A QM/MM Study. J Phys Chem B. 2009;113:10001–10007. [PubMed]
19. Gaussian 03, Revision C.02. Gaussian, Inc; Pittsburgh PA: 2003. For full citation please refer to
20. Gaussian 09, Revision A.01. Gaussian, Inc; Pittsburgh, PA: 2009. For full citation please refer to
21. a) Becke AD. A new mixing of Hartree-Fock and local density-functional theories. J Chem Phys. 1993;98:1372–1377. b) Lee C, Yang W, Parr RG. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys Chem Rev B. 1988;37:785. [PubMed]
22. Hariharan PC, Pople JA. The effect of d-functions on molecular orbital energies for hydrocarbons. Chem Phys Lett. 1972;16:217–219.
23. Tomasi J, Mennucci B, Cammi R. Quantum Mechanical Continuum Solvation Models. Chem Rev. 2005;105:2999–3093. [PubMed]
24. Marenich AV, Cramer CJ, Truhlar DG. Universal solvation model based on solute electron density and a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions. J Phys Chem B. 2009;113:6378–6396. [PubMed]
25. Kryger G, Harel M, Giles K, Toker L, Velan B, Lazar A, Kronman C, Barak D, Ariel N, Shafferman A, Silman I, Sussman JL. Structures of recombinant native and E202Q mutant human acetylcholinesterase complexed with the snake-venom toxin fasciculin-II. Acta Crystallogr, Sect D. 2000;56:1385–1394. [PubMed]
26. Bourne Y, Grassi J, Bougis PE, Marchot P. Conformational flexibility of the acetylcholinesterase tetramer suggested by x-ray crystallography. J Biol Chem. 1999;274:30370–30376. [PubMed]
27. Case DA, Darden TA, Cheatham TE, III, Simmerling CL, Wang J, Duke RE, Luo R, Crowley M, Walker RC, Zhang W, Merz KM, Wang B, Hayik S, Roitberg A, Seabra G, Kolossváry I, Wong KF, Paesani F, Vanicek J, Wu X, Brozell SR, Steinbrecher T, Gohlke H, Yang L, Tan C, Mongan J, Hornak V, Cui G, Mathews DH, Seetin MG, Sagui C, Babin V, Kollman PA. AMBER. University of California; San Francisco: 2008. p. 10.
28. Zhang Y, Kua J, McCammon JA. Role of the Catalytic Triad and Oxyanion Hole in Acetylcholinesterase Catalysis: An ab initio QM/MM Study. J Am Chem Soc. 2002;124:10572–10577. [PubMed]
29. Dolinsky TJ, Nielsen JE, McCammon JA, Baker NA. PDB2PQR: an automated pipeline for the setup of Poisson-Boltzmann electrostatics calculations. Nucleic Acids Research. 2004;32:W665–W667. [PMC free article] [PubMed]
30. Hornberg A, Tunemalm AK, Ekstrom F. Crystal Structures of Acetylcholinesterase in Complex with Organophosphorus Compounds Suggest that the Acyl Pocket Modulates the Aging Reaction by Precluding the Formation of the Trigonal Bipyramidal Transition State. Biochemistry. 2007;46:4815–4825. [PubMed]
31. a) ChemShell, a Computational Chemistry Shell. see b) Sherwood P, de Vries AH, Guest MF, Schreckenbach G, Catlow CRA, French SA, Sokol AA, Bromley ST, Thiel W, Turner AJ, Billeter S, Terstegen F, Thiel S, Kendrick J, Rogers SC, Casci J, Watson M, King F, Karlsen E, Sjøvoll M, Fahmi A, Schäfer A, Lennartz Ch. QUASI: A general purpose implementation of the QM/MM approach and its application to problems in catalysis. J Mol Struct (Theochem) 2003;632:1–23.
32. a) Ahlrichs R, Bär M, Häser M, Horn H, Kölmel C. Electronic structure calculations on workstation computers: The program system Turbomole. Chem Phys Lett. 1989;162:165–169. b) For the current version of TURBOMOLE. see
33. Smith W, Forester TR. DL_POLY_2.0: A general-purpose parallel molecular dynamics simulation package. J Mol Graphics. 1996;14:136–141. [PubMed]
34. MacKerrell AD, Jr, Bashford D, Bellott M, Dunbrack RL, Evanseck JD, Field MJ, Fischer S, Gao J, Guo H, Ha S, Joseph-McCarthy D, Kuchnir L, Kuczera K, Lau FTK, Mattos C, Michnick S, Ngo T, Nguyen DT, Prodhom B, Reiher WE, III, Roux B, Schlenkrich M, Smith JC, Stote R, Straub J, Watanabe M, Wiórkiewicz-Kuczera J, Yin D, Karplus M. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J Phys Chem B. 1998;102:3586–3616. [PubMed]
35. Becke AD. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys Rev A. 1988;38:3098–3100. [PubMed]
36. Weigend F, Ahlrichs R. Balanced basis sets of split valence, triple zeta valence, and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys Chem Chem Phys. 2003;7:3297–3305. [PubMed]