Periodic Boundary Conditions in QM/MM Calculations: Implementation and Tests

Assessment of Free Energies From Electrostatic Embedding Density Functional Tight Binding-Based/Molecular Mechanics in Periodic Boundary Conditions

Electrostatic embedding quantum mechanics/molecular mechanics (QM/MM) methods in periodic boundary conditions (PBC) can successfully describe the condensed phase reactivity of a fragment treated at the QM level with an atomistic description of an electrostatic environment treated at the MM level. The computational cost of ab initio QM methods limits the phase space sampling, thus affecting statistical quantities like free energies. Here, we describe the implementation of a PBC-adapted QM/MM model based on the semi-empirical density-functional based tight-binding (DFTB) method within the GAMESS-US quantum package interfaced with Tinker. Further, we take advantage of the free energy methods provided by a newly developed interface with the PLUMED plugin. The versatility of the implementation is illustrated by the prediction of the free energy profile for three different families of reactions in solution. Overall, using the DFTB/MM, it has been possible to obtain results that are at least in a qualitatively agreement with respect to the experimental data or high-level ab initio simulations.

High-performance QM/MM Enhanced Sampling Molecular Dynamics Simulations with GENESIS SPDYN and QSimulate-QM

A new module for quantum mechanical/molecular mechanical (QM/MM) calculations is implemented in a molecular dynamics (MD) program, SPDYN in GENESIS, interfaced with an electronic structure program, QSimulate-QM. The periodic boundary condition (PBC) in QM/MM simulation is incorporated via QM calculation in real space with duplicated MM charges and particle mesh Ewald (PME) calculation with QM and MM charges. A highly optimized code is implemented in QSimulate-QM, particularly for the density functional tight-binding (DFTB) method, where the interaction between the QM and MM regions is computed utilizing multipole expansions. Together with highly parallelized algorithms in SPDYN, the developed program performs MD simulations based on DFTB in the QM size of ∼100 atoms and MM of ∼100,000 atoms with a better performance than 1 ns/day using one computer node. This feature paves the way for QM/MM-MD enhanced sampling simulations. Various enhanced sampling methods in GENESIS, namely, generalized replica exchange solute tempering (gREST), replica-exchange umbrella sampling (REUS), and path sampling with the string method, are demonstrated at the QM/MM level to compute the Ramachandran plot of alanine dipeptide and the potential of mean force (PMF) of a proton transfer reaction in an enzyme.

An algorithm for very high pressure molecular dynamics simulations

We describe a method to run simulations of ground or excited state dynamics under extremely high pressures. The method is based on the introduction of a fictitious ideal gas that exerts the required pressure on a molecular sample and is therefore called XP-GAS (eXtreme Pressure by Gas Atoms in a Sphere). The algorithm is most suitable for approximately spherical clusters of molecules described by quantum chemistry methods, Molecular Mechanics or mixed QM/MM approaches. We compare the results obtained by the algorithm here presented and by the XP-PCM approach, based on a continuum description of the environment. As a test case, we study the conformational dynamics of 1,3-butadiene either as an isolated molecule ("naked" butadiene) or embedded in a cluster of argon atoms, under pressures up to 15 GPa. Overall, our results show that the XP-GAS QM/MM simulation method is in good agreement with the XP-PCM QM/Continuum model (Cammi model) in describing the effect of the pressure on static properties as the equilibrium geometry of butadiene in the ground state. Furthermore, the comparison of XP-GAS simulations with naked butadiene and butadiene in argon shows the importance, for XP-GAS and related methods, of a realistic representation of the medium in modelling pressure effects.

Semirational Design Strategy To Enhance the Thermostability and Catalytic Activity of Cytochrome P450 105D7 for the Degradation of the Pharmaceutically Active Compounds: Diclofenac

Pharmaceutically active compounds are an important category of emerging pollutants, and their biological transformation processes in the environment are crucial for understanding and evaluating the migration, transformation, and environmental fate of emerging pollutants. The cytochrome P450 105 enzyme family has been proven to play an important role in the degradation of exogenous environmental pollutants. However, its thermostability and catalytic activity still need to be improved to better adapt to complex environmental conditions. This work elucidates the key mechanisms and important residues of the degradation reaction through multiple computational strategies, establishes a mutation library, and obtains 21 single-point mutation designs. Experimental verification showed that 16 single mutants had enhanced thermostability, with the R89F and L197Y mutants showing the highest increases in thermostability at 135 and 119% relative to the wild-type enzyme, respectively. Additionally, as a result of the higher specific activity of D390Q, it was selected for combination mutagenesis, ultimately resulting in three combination mutants (R89F/L197Y, R89F/D390Q, and R89F/L197Y/D390Q) with enhanced thermostability and catalytic activity. This study provides a modification approach for constructing efficient enzyme variants through semirational design and can contribute to the development of control technologies for emerging pollutants.

PyDFT-QMMM: A modular, extensible software framework for DFT-based QM/MM molecular dynamics

PyDFT-QMMM is a Python-based package for performing hybrid quantum mechanics/molecular mechanics (QM/MM) simulations at the density functional level of theory. The program is designed to treat short-range and long-range interactions through user-specified combinations of electrostatic and mechanical embedding procedures within periodic simulation domains, providing necessary interfaces to external quantum chemistry and molecular dynamics software. To enable direct embedding of long-range electrostatics in periodic systems, we have derived and implemented force terms for our previously described QM/MM/PME approach [Pederson and McDaniel, J. Chem. Phys. 156, 174105 (2022)]. Communication with external software packages Psi4 and OpenMM is facilitated through Python application programming interfaces (APIs). The core library contains basic utilities for running QM/MM molecular dynamics simulations, and plug-in entry-points are provided for users to implement custom energy/force calculation and integration routines, within an extensible architecture. The user interacts with PyDFT-QMMM primarily through its Python API, allowing for complex workflow development with Python scripting, for example, interfacing with PLUMED for free energy simulations. We provide benchmarks of forces and energy conservation for the QM/MM/PME and alternative QM/MM electrostatic embedding approaches. We further demonstrate a simple example use case for water solute in a water solvent system, for which radial distribution functions are computed from 100 ps QM/MM simulations; in this example, we highlight how the solvation structure is sensitive to different basis-set choices due to under- or over-polarization of the QM water molecule's electron density.

Analytic Gradients for the Electrostatic Embedding QM/MM Model in Periodic Boundary Conditions Using Particle-Mesh Ewald Sums and Electrostatic Potential Fitted Charge Operators

Long-range electrostatic effects are fundamental for describing chemical reactivity in the condensed phase. Here, we present the methodology of an efficient quantum mechanical/molecular mechanical (QM/MM) model in periodic boundary conditions (PBC) compatible with QM/MM boundaries at chemical bonds. The method combines electrostatic potential fitted charge operators and electrostatic potentials derived from the smooth particle-mesh Ewald (PME) sum approach. The total energy and its analytic first derivatives with respect to QM, MM, and lattice vectors allow QM/MM molecular dynamics (MD) in the most common thermodynamic ensembles. We demonstrate the robustness of the method by performing a QM/MM MD equilibration of methanol in water. We simulate the cis/trans isomerization free-energy profiles in water of proline amino acid and a proline-containing oligopeptide, showing a correct description of the reaction barrier. Our PBC-compatible QM/MM model can efficiently be used to study the chemical reactivity in the condensed phase and enzymatic catalysis.

Cysteine Radical and Glutamate Collaboratively Enable C-H Bond Activation and C-N Bond Cleavage in a Glycyl Radical Enzyme HplG

Hydroxyprolines are abundant in nature and widely utilized by many living organisms. Isomerization of trans-4-hydroxy-d-proline (t4D-HP) to generate 2-amino-4-ketopentanoate has been found to need a glycyl radical enzyme HplG, which catalyzes the cleavage of the C-N bond, while dehydration of trans-4-hydroxy-l-proline involves a homologous enzyme of HplG. Herein, molecular dynamics simulations and quantum mechanics/molecular mechanics (QM/MM) calculations are employed to understand the reaction mechanism of HplG. Two possible reaction pathways of HplG have been explored to decipher the origin of its chemoselectivity. The QM/MM calculations reveal that the isomerization proceeds via an initial hydrogen shift from the Cγ site of t4D-HP to a catalytic cysteine radical, followed by cleavage of the Cδ-N bond in t4D-HP to form a radical intermediate that captures a hydrogen atom from the cysteine. Activation of the Cδ-H bond in t4D-HP to bring about dehydration of t4D-HP possesses an extremely high energy barrier, thus rendering the dehydration pathway implausible in HplG. On the basis of the current calculations, conserved residue Glu429 plays a pivotal role in the isomerization pathway: the hydrogen bonding between it and t4D-HP weakens the hydroxyalkyl Cγ-Hγ bond, and it acts as a proton acceptor to trigger the cleavage of the C-N bond in t4D-HP. Our current QM/MM calculations rationalize the origin of the experimentally observed chemoselectivity of HplG and propose an H-bond-assisted bond activation strategy in radical-containing enzymes. These findings have general implications on radical-mediated enzymatic catalysis and expand our understanding of how nature wisely and selectively activates the C-H bond to modulate catalytic selectivity.

Polarizable Embedding without Artificial Boundary Polarization

We present a fully self-consistent polarizable embedding (PE) model that does not suffer from unphysical boundary polarization. This is achieved through the use of the minimum-image convention (MIC) in the induced electrostatics. It is a simple yet effective approach that includes a more physically accurate description of the polarization throughout the molecular system. Using PE with MIC (PE-MIC), we shed new light on the limitations of commonly employed cutoff models, such as the droplet model, when used in PE calculations. Specifically, we investigate the effects of the unphysical polarization at the outer boundary by comparing induced dipoles and the associated electrostatic potentials, as well as some optical properties of solute-solvent and biomolecular systems. We show that the magnitude of the inaccuracies caused by the unphysical polarization depends on multiple parameters: the nature of the quantum subsystem and of the environment, the cutoff model and distance, and the calculated property.

An efficient electrostatic embedding QM/MM method using periodic boundary conditions based on particle-mesh Ewald sums and electrostatic potential fitted charge operators

Hybrid quantum mechanics/molecular mechanics (QM/MM) models are successful at describing the properties and reactivity of biological macromolecules. Combining ab initio QM/MM methods and periodic boundary conditions (PBC) is currently the optimal approach for modeling chemical processes in an infinite environment, but frequently, these models are too time-consuming for general applicability to biological systems in a solution. Here, we define a simple and efficient electrostatic embedding QM/MM model in PBC, combining the benefits of electrostatic potential fitted atomic charges and particle-mesh Ewald sums, which can efficiently treat systems of an arbitrary size at a reasonable computational cost. To illustrate this, we apply our scheme to extract the lowest singlet excitation energies from a model for Arabidopsis thaliana cryptochrome 1 containing circa 93 000 atoms, accurately reproducing the experimental absorption maximum.

Automatized protocol and interface to simulate QM/MM time-resolved transient absorption at TD-DFT level with COBRAMM

We present a series of new implementations that we recently introduced in COBRAMM, the open-source academic software developed in our group. The goal of these implementations is to offer an automatized workflow and interface to simulate time-resolved transient absorption (TA) spectra of medium-to-big chromophore embedded in a complex environment. Therefore, the excited states absorption and the stimulated emission are simulated along nonadiabatic dynamics performed with trajectory surface hopping. The possibility of treating systems from medium to big size is given by the use of time-dependent density functional theory (TD-DFT) and the presence of the environment is taken into account employing a hybrid quantum mechanics/molecular mechanics (QM/MM) scheme. The full implementation includes a series of auxiliary scripts to properly setup the QM/MM system, the calculation of the wavefunction overlap along the dynamics for the propagation, the evaluation of the transition dipole moment at linear response TD-DFT level, and scripts to setup, run and analyze the TA from an ensemble of trajectories. Altogether, we believe that our implementation will open the door to the easily simulate the time-resolved TA of systems so far computationally inaccessible.

DFT-based QM/MM with Particle-Mesh Ewald for Direct, Long-Range Electrostatic Embedding

We present a DFT-based, QM/MM implementation with long-range electrostatic embedding achieved by direct real-space integration of the particle mesh Ewald (PME) computed electrostatic potential. The key transformation is the interpolation of the electrostatic potential from the PME grid to the DFT quadrature grid, from which integrals are easily evaluated utilizing standard DFT machinery. We provide benchmarks of the numerical accuracy with choice of grid size and real-space corrections, and demonstrate that good convergence is achieved while introducing nominal computational overhead. Furthermore, the approach requires only small modification to existing software packages, as is demonstrated with our implementation in the OpenMM and Psi4 software. After presenting convergence benchmarks, we evaluate the importance of long-range electrostatic embedding in three solute/solvent systems modeled with QM/MM. Water and BMIM/BF4 ionic liquid were considered as ``simple' and ``complex' solvents respectively, with water and p-phenylenediamine (PPD) solute molecules treated at QM level of theory. While electrostatic embedding with standard real-space truncation may introduce negligible error for simple systems such as water solute in water solvent, errors become more significant when QM/MM is applied to complex solvents such as ionic liquids. An extreme example is the electrostatic embedding energy for oxidized PPD in BMIM/BF4 for which real-space truncation produces severe error even at 2-3 nm cutoff distances. This latter example illustrates that utilization of QM/MM to compute redox potentials within concentrated electrolytes/ionic media requires carefully chosen long-range electrostatic embedding algorithms, with our presented algorithm providing a general and robust approach.

First Hyperpolarizability of Water in Bulk Liquid Phase: Long-Range Electrostatic Effects Included via the Second Hyperpolarizability

The molecular first hyperpolarizability β contributes to second-order optical non-linear signals collected from molecular liquids. For the Second Harmonic Generation (SHG) response, the first hyperpolarizability β (2ω,ω,ω) often depends on...

Energy, Structures, and Response Properties with a Fully Coupled QM/AMOEBA/ddCOSMO Implementation

We present the implementation of a fully coupled polarizable QM/MM/continuum model based on the AMOEBA polarizable force field and the domain decomposition implementation of the conductor-like screening model. Energies, response properties, and analytical gradients with respect to both QM and MM nuclear positions are available, and a generic, atomistic cavity can be employed. The model is linear scaling in memory requirements and computational cost with respect to the number of classical atoms and is therefore suited to model large, complex systems. Using three variants of the green-fluorescent protein, we investigate the overall computational cost of such calculations and the effect of the continuum model on the convergence of the computed properties with respect to the size of the embedding. We also demonstrate the fundamental role of polarization effects by comparing polarizable and nonpolarizable embeddings to fully QM ones.

An enhanced sampling QM/AMOEBA approach: The case of the excited state intramolecular proton transfer in solvated 3-hydroxyflavone

We present an extension of the polarizable quantum mechanical (QM)/AMOEBA approach to enhanced sampling techniques. This is achieved by connecting the enhanced sampling PLUMED library to the machinery based on the interface of Gaussian and Tinker to perform QM/AMOEBA molecular dynamics. As an application, we study the excited state intramolecular proton transfer of 3-hydroxyflavone in two solvents: methanol and methylcyclohexane. By using a combination of molecular dynamics and umbrella sampling, we find an ultrafast component of the transfer, which is common to the two solvents, and a much slower component, which is active in the protic solvent only. The mechanisms of the two components are explained in terms of intramolecular vibrational redistribution and intermolecular hydrogen-bonding, respectively. Ground and excited state free energies along an effective reaction coordinate are finally obtained allowing for a detailed analysis of the solvent mediated mechanism.

Exploring the Minimum-Energy Pathways and Free-Energy Profiles of Enzymatic Reactions with QM/MM Calculations

Understanding molecular mechanisms of enzymatic reactions is of vital importance in biochemistry and biophysics. Here, we introduce new functions of hybrid quantum mechanical/molecular mechanical (QM/MM) calculations in the GENESIS program to compute the minimum-energy pathways (MEPs) and free-energy profiles of enzymatic reactions. For this purpose, an interface in GENESIS is developed to utilize a highly parallel electronic structure program, QSimulate-QM (https://qsimulate.com), calling it as a shared library from GENESIS. Second, algorithms to search the MEP are implemented, combining the string method (E et al. J. Chem. Phys. 2007, 126, 164103) with the energy minimization of the buffer MM region. The method implemented in GENESIS is applied to an enzyme, triosephosphate isomerase, which converts dihyroxyacetone phosphate to glyceraldehyde 3-phosphate in four proton-transfer processes. QM/MM-molecular dynamics simulations show performances of greater than 1 ns/day with the density functional tight binding (DFTB), and 10-30 ps/day with the hybrid density functional theory, B3LYP-D3. These performances allow us to compute not only MEP but also the potential of mean force (PMF) of the enzymatic reactions using the QM/MM calculations. The barrier height obtained as 13 kcal mol-1 with B3LYP-D3 in the QM/MM calculation is in agreement with the experimental results. The impact of conformational sampling in PMF calculations and the level of electronic structure calculations (DFTB vs B3LYP-D3) suggests reliable computational protocols for enzymatic reactions without high computational costs.

Computational investigations of allostery in aromatic amino acid biosynthetic enzymes

Allostery, in which binding of ligands to remote sites causes a functional change in the active sites, is a fascinating phenomenon observed in enzymes. Allostery can occur either with or without significant conformational changes in the enzymes, and the molecular basis of its mechanism can be difficult to decipher using only experimental techniques. Computational tools for analyzing enzyme sequences, structures, and dynamics can provide insights into the allosteric mechanism at the atomic level. Combining computational and experimental methods offers a powerful strategy for the study of enzyme allostery. The aromatic amino acid biosynthesis pathway is essential in microorganisms and plants. Multiple enzymes involved in this pathway are sensitive to feedback regulation by pathway end products and are known to use allostery to control their activities. To date, four enzymes in the aromatic amino acid biosynthesis pathway have been computationally investigated for their allosteric mechanisms, including 3-deoxy-d-arabino-heptulosonate 7-phosphate synthase, anthranilate synthase, chorismate mutase, and tryptophan synthase. Here we review the computational studies and findings on the allosteric mechanisms of these four enzymes. Results from these studies demonstrate the capability of computational tools and encourage future computational investigations of allostery in other enzymes of this pathway.

A simplified charge projection scheme for long-range electrostatics in ab initio QM/MM calculations

In a previous work [Pan et al., Molecules 23, 2500 (2018)], a charge projection scheme was reported, where outer molecular mechanical (MM) charges [>10 Å from the quantum mechanical (QM) region] were projected onto the electrostatic potential (ESP) grid of the QM region to accurately and efficiently capture long-range electrostatics in ab initio QM/MM calculations. Here, a further simplification to the model is proposed, where the outer MM charges are projected onto inner MM atom positions (instead of ESP grid positions). This enables a representation of the long-range MM electrostatic potential via augmentary charges (AC) on inner MM atoms. Combined with the long-range electrostatic correction function from Cisneros et al. [J. Chem. Phys. 143, 044103 (2015)] to smoothly switch between inner and outer MM regions, this new QM/MM-AC electrostatic model yields accurate and continuous ab initio QM/MM electrostatic energies with a 10 Å cutoff between inner and outer MM regions. This model enables efficient QM/MM cluster calculations with a large number of MM atoms as well as QM/MM calculations with periodic boundary conditions.

Polarizable embedding QM/MM: the future gold standard for complex (bio)systems?

We provide a perspective of the induced dipole formulation of polarizable QM/MM, showing how efficient implementations will enable their application to the modeling of dynamics, spectroscopy, and reactivity in complex biosystems.

A benchmark for the size of the QM system required for accurate hybrid QM/MM calculations on the metal site of the protein copper, zinc superoxide dismutase

The protein superoxide dismutase 1 (SOD1) is a copper and zinc-binding protein that has been implicated in the pathogenesis of amyotrophic lateral sclerosis (ALS). The Zn(II) binding to SOD1 is critical for the stability of the protein, and has been by itself implicated in ALS pathogenesis. Hence, the quantum mechanical (QM) study of the Zn(II)-site of SOD1 is relevant for understanding ALS. The hybrid QM-molecular mechanics (QM/MM) approach commonly employed for the QM study of proteins is highly dependent on the size of the sub-system treated quantum-mechanically. The size of the QM system also determines the computational feasibility of a given method. In the present work, we compare optimized geometries for the metal site and Zn(II) dissociation energies obtained with a QM/MM methodology employing different sizes for the QM sub-system. We find that geometries converge rapidly to RMSDs of around 0.3 Å, and fails to converge further, while a QM system of 480 atoms was required for converging the Zn(II) interaction energy of SOD1 to within 5 kcal*mol^-1, and a 611-atoms QM system for a 1 kcal*mol^-1 convergence with respect to our reference, 1280 QM-atoms system. Graphical Abstract The size of the QM system is critical for both the accuracy and the computational cost of a QM/MM calculation. We have identified a optimum balance for the study of the active site of the coppper, zinc superoxide dismutase.

Multistep Explicit Solvation Protocol for Calculation of Redox Potentials

The calculation of molecular redox potentials in aqueous solution presents a challenge to quantum chemistry due to the need to calculate charged, open-shell species experiencing large solvent effects. Traditionally, redox potentials are calculated via the use of density functional theory and continuum solvation methods, but such protocols have been found to often suffer from large errors, particularly in the case of aqueous solution. While explicit solvation models hold promise of higher accuracy to describe solvent effects in general, their complicated use and lack of well-defined, reliable protocols has hindered their adoption. In this study, we present an explicit-solvation-based approach for the calculation of molecular redox potentials. We combine the use of affordable semiempirical QM/MM molecular dynamics (making use of the recently proposed GFN-xTB method by Grimme et al.) for both redox states and use the linear response approximation to relate vertical ionization energies to the adiabatic redox potential. Simulation length, averaging over snapshots, and accounting for bulk and polarization effects are systematically evaluated using phenol as a working example. We find that it is crucial to reliably account for bulk solvation effects in these calculations, as well as polarization effects which we divide up into short-range and long-range contributions. The short-range polarization contribution is accounted for via QM-region expansion, while the long-range contribution is accounted for via Drude-polarizable QM/MM. Our multistep protocol has been coded to be used in a fully automatic way in a local version of Chemshell. It has been evaluated on a test set of oxidation potentials of organic molecules and found to give gas-solution redox shifts with a mean absolute error of 0.13 eV with respect to experiment, compared to mean absolute errors of 0.26 and 0.21 eV with CPCM and SMD continuum models, respectively.

Representation of the QM Subsystem for Long-Range Electrostatic Interaction in Non-Periodic Ab Initio QM/MM Calculations

In QM/MM calculations, it is essential to handle electrostatic interactions between the QM and MM subsystems accurately and efficiently. To achieve maximal efficiency, it is convenient to adopt a hybrid scheme, where the QM electron density is used explicitly in the evaluation of short-range QM/MM electrostatic interactions, while a multipolar representation for the QM electron density is employed to account for the long-range QM/MM electrostatic interactions. In order to avoid energy discontinuity at the cutoffs, which separate the short- and long-range QM/MM electrostatic interactions, a switching function should be utilized to ensure a smooth potential energy surface. In this study, we benchmarked the accuracy of such hybrid embedding schemes for QM/MM electrostatic interactions using different multipolar representations, switching functions and cutoff distances. For test systems (neutral and anionic oxyluciferin in MM (aqueous and enzyme) environments), the best accuracy was acquired with a combination of QM electrostatic potential (ESP) charges and dipoles and two switching functions (long-range electrostatic corrections (LREC) and Switch) in the treatment of long-range QM/MM electrostatics. It allowed us to apply a 10Å distance cutoff and still obtain QM/MM electrostatics/polarization energies within 0.1 kcal/mol and time-dependent density functional theory (TDDFT)/MM vertical excitation energies within 10-3 eV from theoretical reference values.

Efficient Computation of Free Energy Surfaces of Diels⁻Alder Reactions in Explicit Solvent at Ab Initio QM/MM Level

For Diels–Alder (DA) reactions in solution, an accurate and converged free energy (FE) surface at ab initio (ai) quantum mechanical/molecular mechanical (QM/MM) level is imperative for the understanding of reaction mechanism. However, this computation is still far too expensive. In a previous work, we proposed a new method termed MBAR+wTP, with which the computation of the ai FE profile can be accelerated by several orders of magnitude via a three-step procedure: (I) an umbrella sampling (US) using a semi-empirical (SE) QM/MM Hamiltonian is performed; (II) the FE profile is generated using the Multistate Bennett Acceptance Ratio (MBAR) analysis; and (III) a weighted Thermodynamic Perturbation (wTP) from the SE Hamiltonian to the ai Hamiltonian is performed to obtain the ai QM/MM FE profile using weight factors from the MBAR analysis. In this work, this method is extended to the calculations of two-dimensional FE surfaces of two Diels–Alder reactions of cyclopentadiene with either acrylonitrile or 1-4-naphthoquinone at ai QM/MM level. The accurate activation free energies at the ai QM/MM level, which are much closer to the experimental measurements than those calculated by other methods, indicate that this MBAR+wTP method can be applied in the studies of complex reactions in condensed phase with much-enhanced efficiency.

Dioxygen binding to Fe-MOF-74: microscopic insights from periodic QM/MM calculations

Accurate MOF-FF parameter sets were determined for the ferrous and ferric forms of an iron-based metal–organic framework (MOF) called Fe-MOF-74. For this purpose, density functional theory (DFT) calculations were applied to truncated cluster models of Fe-MOF-74, and the DFT-calculated geometries and energy derivatives were used for the force-field parameterization. The resultant parameter sets performed remarkably well in reproducing the experimentally determined structure of the MOF. We also performed periodic quantum mechanics (QM) / molecular mechanics (MM) calculations employing a subtractive scheme called ONIOM, with the optimized MOF-FF parameters used for the MM calculations, in an attempt to evaluate the binding energies between O2 and several Fe-MOF-74 variants. The calculated binding energy for Fe-MOF-74 agreed very well with the experimental value, and QM/MM geometry optimization calculations confirmed that the O2-bound complex has a side-on geometry. Our calculations also predicted that, when the two neighboring iron ions around the O2-binding site are replaced with other metal ions (Mg2+, Ni2+, Zn2+, Co2+, or Mn2+), there are noticeable variations in the binding energy, indicating that these substituted metal ions affect the O2 binding indirectly.