Figures
Abstract
The frozen domain effective fragment molecular orbital method is extended to allow for the treatment of a single fragment at the MP2 level of theory. The approach is applied to the conversion of chorismate to prephenate by Chorismate Mutase, where the substrate is treated at the MP2 level of theory while the rest of the system is treated at the RHF level. MP2 geometry optimization is found to lower the barrier by up to 3.5 kcal/mol compared to RHF optimzations and ONIOM energy refinement and leads to a smoother convergence with respect to the basis set for the reaction profile. For double zeta basis sets the increase in CPU time relative to RHF is roughly a factor of two.
Citation: Christensen AS, Steinmann C, Fedorov DG, Jensen JH (2014) Hybrid RHF/MP2 Geometry Optimizations with the Effective Fragment Molecular Orbital Method. PLoS ONE 9(2): e88800. https://doi.org/10.1371/journal.pone.0088800
Editor: Freddie Salsbury Jr., Wake Forest University, United States of America
Received: October 28, 2013; Accepted: January 15, 2014; Published: February 18, 2014
Copyright: © 2014 Christensen et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: ASC is funded by the Novo Nordisk STAR PhD program. DGF is supported by the Next Generation Super Computing Project, Nanoscience Program (NEXT, Japan) and Computational Materials Science Initiative (CMSI, Japan). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Combined quantum mechanical/molecular mechanical (QM/MM) and fragment-based quantum mechanical methods [1]–[13], recently reviewed [13], [14], are becoming increasingly popular for large molecular systems. In the fragment molecular orbital method (FMO) [15]–[17] one does fragment calculations in the presence of the embedding potential of all the other fragments, whereas in the recently developed effective fragment molecular orbital method (EFMO) [18], [19] fragment polarizabilities are used instead to approximate the many-body polarization.
For fast geometry optimizations, FMO with the frozen domain and dimers (FDD) [20] has been proposed and EFMO/FDD has been used to map the reaction path of the conversion of chorismate to prephanate in Chorismate Mutase at the RHF level for geometry optimization combined with ONIOM for energy refinement [21]. Chorismate Mutase has also been studied extensively by many groups. Particularly, the group of Mulholland has invested considerable amount of resources to evaluate the catalytic effect of Chorismate Mutase [22]–[28]. Other related QM/MM work on Chorismate Mutase includes FMO energetics refinement by Ishida et al. [29] and the work of Claeyssens et al.[30] who used linear scaling coupled cluster methods to obtain the reaction barrier on structures optimized using a QM/MM approach with density functional theory used to describe the QM region. This study specifically underlines the importance of energy corrections at a correlated level of theory, which in turn requires reliable optimization of the molecular structure. Our recent study [21] emphasizes that in addition to a high quality energy evaluation on the reaction complex, a conformational sampling of the reaction complex geometry is needed in order to obtain a reliable energy barrier, since the reaction barrier can fluctuate by up to 15 kcal/mol between geometry optimizations on different starting conformations.
Our previous methodology was to estimate the reaction barrier in Chorismate Mutase using an EFMO-RHF geometry optimization with an ONIOM MP2 energy correction [21]. It was clear, however, that the RHF based optimization did not always lead to a reliable MP2 correction. For example, the MP2 energy did not converge in a smooth manner with respect to the basis set size. One likely explanation is, that it is in general not optimal to deal with reaction complexes for which the structure is calculated using an uncorrelated wave function method such as RHF.
In this work, we have created a method to obtain a correlated (MP2 level) reaction complex geometry using the EFMO method on a large system, and show a calculation of the reaction barrier in Chorismate Mutase as an example.
We extend EFMO/FDD to enable treatment of only one fragment at the MP2 level and show that it is a good compromise between efficiency and accuracy. Note that the effects of conformational sampling are not investigated in this paper.
This paper is organized as follows: First we present the EFMO method and our extension to the EFMO energy and gradient. Second we compare our method to similar ONIOM calclations on the reaction barrier of the conversion of chorismate to prephanate in Chorismate Mutase.
Theory
The basics of EFMO can be summarized as follows. The system is divided into fragments and we use the adaptive frozen orbital technique (AFO) [31] to treat fragment boundaries by freezing the molecular orbitals corresponding to detached covalent bonds. Ab initio calculations of fragments are carried out without embedding, and the total polarization is evaluated using fragment polarizabilities. In the next step, ab initio calculations of dimers are carried out to account for two-body quantum effects such as the charge transfer between fragment pairs within a cut-off distance, . The total energy in the two-body EFMO expansion is then:
(1)
Here is the quantum mechanical gas-phase energy of each monomer fragment,
is the quantum mechanical two-body polarization energy between two fragments,
is the classical two-body polarization energy between two fragments, and
is the classical polarization energy of the system.
In the frozen domain method (FD) [20], the geometry of the molecular system is optimized only for a smaller subsystem called the active domain, while the atoms in the rest of the system are fixed.
For a given molecular system, we define two domains (“frozen”) and
(“active”). Domain
is defined as all atoms having a frozen geometry and domain
is defined as all atoms whose positions are optimized. Each domain is further divided into a number of molecular fragments. In the frozen domain and dimers methods (FDD) [20], the domain with frozen geometry is further divided as a polarizable domain with frozen geometry,
and a domain with frozen geometry and fragment electron densities that are not updated after they have been calculated the first time. The EFMO energy [21] is then given by:
(2)where
and
are the internal energies of domains
and
, respectively,
is the interaction between domains
and
,
is the interaction between domains
and
and
is the classical total polarization energy of the whole system. In our EFMO-RHF:MP2 extension, we evaluate the internal energies of domain
and
at the RHF level. Furthermore, we specify a single fragment
(“high level”) from the active domain to be treated at the MP2 level of theory (see Fig. 1 for a schematic overview). The total EFMO-RHF:MP2 energy is then given as
(3)where
is the MP2 correlation energy of fragment
.
denotes the frozen domain (green);
denotes the polarizable domain (blue);
denotes the active domain (red);
denotes fragment
, for which the MP2 energy and gradients are evaluated (yellow).
The corresponding EFMO energy gradients of each domain in the FDD approximation:(4)
(6)This gives the following EFMO-RHF:MP2 energy gradients:
(7)
(9)Where
contains the gradient of the MP2 correlation energy for fragment H
A.
Methods
All calculations were carried out in a development version of GAMESS [32] where FMO and EFMO are implemented [33].
Starting structures for Chorismate Mutase were obtained from Steinmann et al. [21] who prepared the structures following Claeyssens et al. [28]. The preparation can be summarized as follows: The experimental structure of Chorismate Mutase was obtained from the Protein Data Bank (PDB code: 2CHT) and protonated using PDB2PQR at pH 7. The inhibitors were manually replaced with Chorismate in the reactant state. The complexes were simulated in GROMACS with the CHARMM27 force field at 300K. The structure was then prepared for fragment based calculations in FragIt. [34] All residues with an atom within a distance of 2.0 Å from any atom in chorismate were assigned to the A (active) domain. All atoms in the prephanate/chorismate reaction complex were assigned to the H fragment. See Fig. 2. The total system consists of 313 fragments, divided as 213 fragments in the frozen domain, 92 fragment in the polarizable
domain, and 8 fragments in active
domain of which one fragment (the reaction complex) is treated at a higher level, i.e. in the
domain.
(A) shows the atoms and bonds of the active layer, , with the MP2 fragment,
, highlighted in yellow. (B) additionally shows the polarizable, but with frozen geometry, buffer layer
(in blue) surrounding the active layer. (C) additionally shows the
layer in which both geometry and densities are frozen.
The adiabatic mapping was carried out using the presented EFMO-RHF:MP2 gradient with 6-31G(d) basis set on all atoms. Two additional runs were also carried out, in these cases with the cc-pVDZ or cc-pVTZ on chorismate and 6-31G(d) on remaining atoms. The EFMO-RHF/6-31G(d):MP2//cc-pVTZ reaction path was obtained starting from the converged structures in the EFMO-RHF/6-31G(d):MP2//cc-pVDZ reaction path.
The RESDIM keyword was set to 1.5 and the optimization convergence criterion was set to Hartree/Bohr. Each step of the reaction path was obtained by imposing harmonic constraints on
and
with a force constant of 500 kcal/Å. The FDD approximation was enabled by setting MODFD = 3 in all calculations. GAMESS input files to calculate the reaction path at the EFMO-RHF/6-31G(d):MP2/cc-pVDZ level of theory in File S1.
Timings for the optimization procedure were carried out on 80 Intel Xeon X5550 CPU cores distributed across 10 nodes and the Generalized Distributed Data Interface (GDDI) was used to run the code in parallel [35].
Results
Transition State Structure
We define the reaction coordinate similarly to Claeyssens et al.[28] as the difference in bond length between the breaking O2-C1 bond and the forming C4–C3 bond in chorismate, i.e. (see Fig. 3). The reaction coordinate of the transition state was found to be −0.17 Å using the 6-31G(d) basis set on the MP2 fragment and −0.43 Å for both the cc-pVTZ and cc-pVDZ basis set reaction paths. The convergence with respect to basis set is in good, quantitative agreement with the coordinates obtained by Claeyssens et al. [30] using a QM/MM approach, treating the reaction complex at the LCCSD(T0) level of theory (−0.4 Å). In comparison, the corresponding MP2:RHF ONIOM calculations by Steinmann et al. [21] resulted in transition state reaction coordinates of 0.13, −0.36, and 0.13 Å with the cc-pVDZ, cc-pVTZ and cc-pVQZ basis sets used in the MP2 calculation, respectively.
The atoms describing the reaction coordinate are marked with numbers one through four [21].
The reaction coordinate found using MP2 to optimize the reaction complex substantially improves obtained the transition state structure compared to our MP2:RHF ONIOM approach and is in good agreement with a high-level calculation [30].
Reaction Barrier
Electronic energy barriers and reaction coordinates for the transition state are given in Table 1 and Fig. 4. We find the electronic energy barrier at the EFMO-RHF/6-31G(d):MP2/6-31G(d) level of theory to be 20.95 kcal/mol. Increasing the size of the basis set on the MP2 fragment decreases the barrier to 19.21 kcal/mol with the cc-pVDZ basis set and 18.34 kcal/mol with the cc-pVTZ basis set.
The three reation paths are calculated using the FDD/EFMO-RHF:MP2 approach with three different basis sets on the reaction complex in the MP2 layer. The 6-31G(d) basis set was used for the RHF layer in all three cases.
In comparison, the corresponding MP2:RHF ONIOM calculations by Steinmann et al. resulted in barriers of 22.24, 19.75, and 21.79 kcal/mol with the cc-pVDZ, cc-pVTZ and cc-pVQZ basis sets, respectively. In contrast to the ONIOM approach, we find that for increasing basis set sizes, the electronic energy barrier is systematically reduced.
The experimental enthalpy barrier has been measured to be 12.7 kcal/mol.[28], [36]. The large difference between the calculated reaction barrier and the experimentally measured barrier is likely caused by lack of conformational sampling and the relatively small size of the active geometry region, . We have previously shown [21] that the reaction barrier can fluctuate by up to 15 kcal/mol between different conformational samples. A more accurate estimation of the reaction barrier (compared to the experimental value) using this approach would thus likely require averaging the barrier over a large number of conformational samples with a larger active region.
Reaction Energy
The energy difference between the product and reactant state is found to be −3.2 kcal/mol using the 6-31G(d) basis set on chorismate. Increasing the basis set to cc-pVDZ and cc-pVTZ on chorismate decreased the reaction energy to −6.83 kcal/mol and −6.17 kcal/mol, respectively. The ONIOM approach by Steinmann et al. found the reaction energy to be between −5.48 kcal/mol to −0.82 kcal/mol. However, in the ONIOM approach increasing the basis set from cc-pVTZ on chorismate increased the reaction energy from −5.48 kcal/mol to −1.17 kcal/mol. We find that all three basis sets are in close agreement, and only a 0.7 kcal/mol difference between the cc-pVDZ and cc-pVTZ reaction paths.
The reaction energy calculated using the presented method is around −6 kcal/mol when the cc-pVDZ or cc-pVTZ are used in the MP2 calculation, which contrasts our earlier calculations where the reaction energy systematically increases as the basis set size is increased (see Table 1).
As we discuss in the previous subsection, a more accurate estimation would likely require averaging this values over a large number of conformational samples and possibly a larger active region.
Timings
Running on 80 cores distributed on 10 compute nodes and using the default compute node load balancing scheme, the average time for a geometry optimization step was 760s at the EFMO-RHF/6-31G(d) level of theory[21]. For the EFMO-RHF/6-31G(d):MP2/6-31G(d) calculation, this time increased to 1526 s per step. Increasing the basis set on the MP2 part of the system to cc-pVDZ and cc-pVTZ increased the time to 1967 s and 18845 s, respectively (see Table 2). The large increase in calculation time from cc-pVDZ to cc-pVTZ was found to be due to sub-optimal load balancing in GDDI during the MP2 part of the calculation. Subsequently, one optimization using the cc-pVTZ was carried out, in which the calculation of the MP2 fragment energy and gradient was distributed across all 10 nodes. This reduced the average gradient step time from 18845 s to 10911 s. In other words, the slower calculation used 10 GDDI groups in the second (MP2) layer, whereas the faster one had 1 group, during the monomer step. The latter run is more efficient because the MP2 fragment was calculated by all 10 nodes, whereas in the former only by 1 node.
A calculation of the reaction at the EFMO-RHF/6-31G(d):MP2/cc-pVDZ level is thus about 2.5 times more expensive than the same calculation at the EFMO-RHF/6-31G(d) level of theory. But as we have shown, the calculated reaction coordinate is essentially the same as that found using a coupled cluster approach[30] when applying the EFMO-RHF/6-31G(d):MP2/cc-pVDZ level of theory.
Conclusion
We have implemented an scheme for optimizing a reaction complex using a correlated method in the EFMO/FDD approximation.[21] Our method is computationally efficient when a moderately sized basis sets is used on the correlated fragment. At the EFMO-RHF/6-31G(d):MP2/cc-pVDZ level of theory, the method is about 2.5 times slower than the same calculation at the EFMO-RHF/6-31G(d) level. However, the use of a correlated method (MP2) in the optimization step substantially improves the calculated transition state compared to similar uncorrelated optimization with a subsequent MP2 ONIOM energy correction.
The modest increase in computational cost compared to an similar uncorrelated calculation makes the presented method very attractive for cases where electron correlation is essential for a correct and reliable geometry optimization. The method is a special case within the FMO or EFMO approximations and thus requires no further approximations, such as those in the ONIOM method, and is carried out in a single calculation in the GAMESS program.
The method is thus a general method to obtain geometry optimized correlated structures inside large molecular systems when using FMO or EFMO. For example the EFMO method has been used to estimate hydrolysis barriers for the enzyme Bacillus circulans xylanase [37].
Our EFMO-RHF:MP2 approach does not achieve chemical accuracy in predicting enthalpy barrier of the conversion of chorismate to prephanate in Chorismate Mutase, but as we show in our previous work this is likely due to the lack of structural sampling [21].
In conclusion, we have demonstrated that our method serves as a rigorous and viable alternative to the widely used ONIOM approach. Source code to add the method to GAMESS can be found at: https://github.com/andersx/efmo-rhf-mp2.
Supporting Information
File S1.
GAMESS input files to calculate the reaction path at the EFMO-RHF/6-31G(d):MP2/cc-pVDZ level of theory.
https://doi.org/10.1371/journal.pone.0088800.s001
(BZ2)
Author Contributions
Conceived and designed the experiments: ASC CSS JHJ. Performed the experiments: ASC CSS DGF. Analyzed the data: ASC CSS JHJ. Contributed reagents/materials/analysis tools: ASC CSS DGF JHJ. Wrote the paper: ASC JHJ.
References
- 1. Otto P, Ladik J (1975) Investigation of the interaction between molecules at medium distances: I. SCF LCAO MO supermolecule, perturbational and mutually consistent calculations for two interacting HF and CH2O molecules. Chem Phys 8: 192–200.
- 2. Gao J (1997) Toward a Molecular Orbital Derived Empirical Potential for Liquid Simulations. J? Phys Chem? B 101: 657–663.
- 3. Korchowiec J, Gu FL, Aoki Y (2005) Elongation method at restricted open-shell Hartree€Fock level of theory. Int? J? Quantum Chem 105: 875–882.
- 4. Chen XH, Zhang JZH (2004) MFCC-downhill simplex method for molecular structure optimization. J? Theor Comput Chem 3: 277–289.
- 5. Gordon MS, Mullin JM, Pruitt SR, Roskop LB, Slipchenko LV, et al. (2009) Accurate methods for large molecular systems. J? Phys Chem? B 113: 9646–9663.
- 6. Flick JC, Kosenkov D, Hohenstein EG, Sherrill CD, Slipchenko LV (2012) Accurate Prediction of Noncovalent Interaction Energies with the Effective Fragment Potential Method: Comparison of Energy Components to Symmetry-Adapted Perturbation Theory for the S22 Test Set. J? Chem Theory Comput 8: 2835–2843.
- 7. Kobayashi M, Nakai H (2012) How does it become possible to treat delocalized and/or open-shell systems in fragmentation-based linear-scaling electronic structure calculations? The case of the divide-and-conquer method. Phys Chem Chem Phys 14: 7629–7639.
- 8. Huang L, Massa L (2012) Electron density study of the anti-Alzheimer€s disease drug donepezil from conventional x-ray data and invariom database application. Future Medicinal Chemistry 4: 1479–1494.
- 9. Söderhjelm P, Kongsted J, Ryde U (2010) Ligand affinities estimated by quantum chemical calculations. J? Chem Theory Comput 6: 1726–1737.
- 10. Sahu N, Yeole SD, Gadre SR (2013) Appraisal of Molecular Tailoring Approach for Large Clusters. J? Chem Phys 138: 104101.
- 11. Frank A, Möller HM, Exner TE (2012) Toward the Quantum Chemical Calculation of NMR Chemical Shifts of Proteins 2: Level of Theory, Basis Set, and Solvents Model Dependence. J? Chem Theory Comput 8: 1480–1492.
- 12. Kurbanov EK, Leverentz HR, Truhlar DG, Amin EA (2012) Electrostatically Embedded Many-Body Expansion for Neutral and Charged Metalloenzyme Model Systems. J? Chem Theory Comput 8: 1–5.
- 13. Gordon MS, Fedorov DG, Pruitt SR, Slipchenko LV (2012) Fragmentation methods: a route to accurate calculations on large systems. Chem Rev 112: 632–672.
- 14. Senn HM, Thiel W (2009) QM/MM Methods for Biomolecular Systems. Angew Chem Int Ed 48: 1198–1229.
- 15. Kitaura K, Ikeo E, Asada T, Nakano T, Uebayasi M (1999) Fragment molecular orbital method: an approximate computational method for large molecules. Chem Phys Lett 313: 701–706.
- 16. Fedorov DG, Kitaura K (2007) Extending the Power of Quantum Chemistry to Large Systems with the Fragment Molecular Orbital Method. J? Phys Chem? A 111: 6904–6914.
- 17. Fedorov DG, Nagata T, Kitaura K (2012) Exploring chemistry with the fragment molecular orbital method. Phys Chem Chem Phys 14: 7562–7577.
- 18. Steinmann C, Fedorov DG, Jensen JH (2010) Effective fragment molecular orbital method: A merger of the effective fragment potential and fragment molecular orbital methods. J? Phys Chem? A 114: 8705–8712.
- 19. Steinmann C, Fedorov DG, Jensen J H (2012) The effective fragment molecular orbital method for fragments connected by covalent bonds. PLoS ONE 7(7): e41117.
- 20. Fedorov DG, Alexeev Y, Kitaura K (2011) Geometry optimization of the active site of a large system with the fragment molecular orbital method. J? Phys Chem Lett 2: 282–288.
- 21. Steinmann C, Fedorov DG, Jensen JH (2013) Mapping Enzymatic Catalysis using the Effective Fragment Molecular Orbital Method: Towards all ab initio Biochemistry. PLoS ONE 8(4): e60602.
- 22. Lyne PD, Mulholland AJ, Richards WG (1995) Insights into chorismate mutase catalysis from a combined QM/MM simulation of the enzyme reaction. J? Am Chem Soc 117: 11345–11350.
- 23. Ranaghan KE, Ridder L, Szefczyk B, Sokalski WA, Hermann JC, et al. (2003) Insights into enzyme catalysis from QM/MM modelling: transition state stabilization in chorismate mutase. Mol Phys 101: 2695–2714.
- 24. Szefczyk B, Mulholland AJ, Ranaghan KE, Sokalski WA (2004) Differential transition-state stabilization in enzyme catalysis: Quantum chemical analysis of interactions in the chorismate mutase reaction and prediction of the optimal catalytic field. J? Am Chem Soc 126: 16148–16159.
- 25. Ranaghan KE, Mulholland AJ (2004) Conformational effects in enzyme catalysis: QM/MM free energy calculation of the ‘NAC’ contribution in chorismate mutase. Chem Commun 10: 1238–1239.
- 26. Claeyssens F, Ranaghan KE, Manby FR, Harvey JN, Mulholland AJ (2005) Multiple high-level QM/MM reaction paths demonstrate transition-state stabilization in chorismate mutase: correlation of barrier height with transition-state stabilization. Chem Commun 40: 5068–5070.
- 27. Szefczyk B, Claeyssens F, Mulholland AJ, Sokalski WA (2007) Quantum chemical analysis of reaction paths in chorismate mutase: Conformational effects and electrostatic stabilization. Int? J? Quantum Chem 107: 2274–2285.
- 28. Claeyssens F, Ranaghan KE, Lawan N, Macrae SJ, Manby FR, et al. (2011) Analysis of chorismate mutase catalysis by QM/MM modelling of enzyme-catalysed and uncatalysed reactions. Org Biomol Chem 9: 1578–1590.
- 29. Ishida T, Fedorov DG, Kitaura K (2006) All electron quantum chemical calculation of the entire enzyme system confirms a collective catalytic device in the chorismate mutase reaction. J? Phys Chem? B 110: 1457–1463.
- 30. Claeyssens F, Harvey JN, Manby FR, Mata RA, Mulholland AJ, et al. (2006) High Accuracy Computation of Reaction Barriers in Enzymes. Angew Chem Int Ed 45: 6856–6859.
- 31. Fedorov DG, Jensen JH, Deka RC, Kitaura K (2008) Covalent bond fragmentation suitable to describe solids in the fragment molecular orbital method J Phys Chem A. 112: 11808–11816.
- 32. Schmidt MW, Baldridge KK, Boatz JA, Elbert ST, Gordon MS, et al. (1993) General Atomic and Molecular Electronic Structure System. J? Comput Chem 14: 1347–1363.
- 33. Fedorov DG, Kitaura K (2004) The importance of three-body terms in the fragment molecular orbital method. J? Chem Phys 120: 6832–6840.
- 34. Steinmann C, Ibsen MW, Hansen AS, Jensen JH (2012) FragIt: A Tool to Prepare Input Files for Fragment Based Quantum Chemical Calculations. PLoS ONE 7(9): e44480.
- 35. Fedorov D, Olson R, Kitaura K, Gordon M, Koseki S (2004) A new hierarchical parallelization scheme: Generalized distributed data interface (gddi), and an application to the fragment molecular orbital method (fmo). J? Comput Chem 25: 872–880.
- 36. Kast P, Asif-Ullah M, Jiang N, Hilvert D (1996) Exploring the active site of chorismate mutase by combinatorial mutagenesis and selection: The importance of electrostatic catalysis. Proc Natl Acad Sci 93: 5043–5048.
- 37. Hediger MR, Steinmann C, De Vico L, Jensen JH (2013) A computational method for the systematic screening of reaction barriers in enzymes: searching for Bacillus circulans xylanase mutants with greater activity towards a synthetic substrate. PeerJ 1: e111.