Redistributed charge and dipole schemes for combined quantum mechanical and molecular mechanical calculations

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.

CHARMM-GUI QM/MM Interfacer for a Quantum Mechanical and Molecular Mechanical (QM/MM) Simulation Setup: 1. Semiempirical Methods

Quantum mechanical (QM) treatments, when combined with molecular mechanical (MM) force fields, can effectively handle enzyme-catalyzed reactions without significantly increasing the computational cost. In this context, we present CHARMM-GUI QM/MM Interfacer, a web-based cyberinfrastructure designed to streamline the preparation of various QM/MM simulation inputs with ligand modification. The development of QM/MM Interfacer has been achieved through integration with existing CHARMM-GUI modules, such as PDB Reader and Manipulator, Solution Builder, and Membrane Builder. In addition, new functionalities have been developed to facilitate the one-stop preparation of QM/MM systems and enable interactive and intuitive ligand modifications and QM atom selections. QM/MM Interfacer offers support for a range of semiempirical QM methods, including AM1(+/d), PM3(+/PDDG), MNDO(+/d, +/PDDG), PM6, RM1, and SCC-DFTB, tailored for both AMBER and CHARMM. A nontrivial setup related to ligand modification, link-atom insertion, and charge distribution is automatized through intuitive user interfaces. To illustrate the robustness of QM/MM Interfacer, we conducted QM/MM simulations of three enzyme-substrate systems: dihydrofolate reductase, insulin receptor kinase, and oligosaccharyltransferase. In addition, we have created three tutorial videos about building these systems, which can be found at https://www.charmm-gui.org/demo/qmi. QM/MM Interfacer is expected to be a valuable and accessible web-based tool that simplifies and accelerates the setup process for hybrid QM/MM simulations.

The development of the QM/MM interface and its application for the on-the-fly QM/MM nonadiabatic dynamics in JADE package: Theory, implementation, and applications

Understanding the nonadiabatic dynamics of complex systems is a challenging task in computational photochemistry. Herein, we present an efficient and user-friendly quantum mechanics/molecular mechanics (QM/MM) interface to run on-the-fly nonadiabatic dynamics. Currently, this interface consists of an independent set of codes designed for general-purpose use. Herein, we demonstrate the ability and feasibility of the QM/MM interface by integrating it with our long-term developed JADE package. Tailored to handle nonadiabatic processes in various complex systems, especially condensed phases and protein environments, we delve into the theories, implementations, and applications of on-the-fly QM/MM nonadiabatic dynamics. The QM/MM approach is established within the framework of the additive QM/MM scheme, employing electrostatic embedding, link-atom inclusion, and charge-redistribution schemes to treat the QM/MM boundary. Trajectory surface-hopping dynamics are facilitated using the fewest switches algorithm, encompassing classical and quantum treatments for nuclear and electronic motions, respectively. Finally, we report simulations of nonadiabatic dynamics for two typical systems: azomethane in water and the retinal chromophore PSB3 in a protein environment. Our results not only illustrate the power of the QM/MM program but also reveal the important roles of environmental factors in nonadiabatic processes.

The Influence of Clustered DNA Damage Containing Iz/Oz and OXOdG on the Charge Transfer through the Double Helix: A Theoretical Study

The genome-the source of life and platform of evolution-is continuously exposed to harmful factors, both extra- and intra-cellular. Their activity causes different types of DNA damage, with approximately 80 different types of lesions having been identified so far. In this paper, the influence of a clustered DNA damage site containing imidazolone (Iz) or oxazolone (Oz) and 7,8-dihydro-8-oxo-2'-deoxyguanosine (OXOdG) on the charge transfer through the double helix as well as their electronic properties were investigated. To this end, the structures of oligo-Iz, d[A1Iz2A3OXOG4A5]*d[T5C4T3C2T1], and oligo-Oz, d[A1Oz2A3OXOG4A5]*d[T5C4T3C2T1], were optimized at the M06-2X/6-D95**//M06-2X/sto-3G level of theory in the aqueous phase using the ONIOM methodology; all the discussed energies were obtained at the M06-2X/6-31++G** level of theory. The non-equilibrated and equilibrated solvent-solute interactions were taken into consideration. The following results were found: (A) In all the discussed cases, OXOdG showed a higher predisposition to radical cation formation, and B) the excess electron migration toward Iz and Oz was preferred. However, in the case of oligo-Oz, the electron transfer from Oz2 to complementary C4 was noted during vertical to adiabatic anion relaxation, while for oligo-Iz, it was settled exclusively on the Iz2 moiety. The above was reflected in the charge transfer rate constant, vertical/adiabatic ionization potential, and electron affinity energy values, as well as the charge and spin distribution. It can be postulated that imidazolone moiety formation within the CDL ds-oligo structure and its conversion to oxazolone can significantly influence the charge migration process, depending on the C2 carbon hybridization sp2 or sp3. The above can confuse the single DNA damage recognition and removal processes, cause an increase in mutagenesis, and harm the effectiveness of anticancer therapy.

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.

Proton Coulomb Blockade Effect Involving Covalent Oxygen-Hydrogen Bond Switching

Instead of the canonical Grotthuss mechanism, we show that a knock-on proton transport process is preferred between organic functional groups (e.g., -COOH and -OH) and adjacent water molecules in biological proton channel and synthetic nanopores through comprehensive quantum and classical molecular dynamics simulations. The knock-on process is accomplished by the switching of covalent O─H bonds of the functional group under externally applied electric fields. The proton transport through the synthetic nanopore exhibits nonlinear current-voltage characteristics, suggesting an unprecedented proton Coulomb blockade effect. These findings not only enhance the understanding of proton transport in nanoconfined systems but also pave the way for the design of a variety of proton-based nanofluidic devices.

Bacteria at Work - Experimental and Theoretical Studies Reveal the Catalytic Mechanism of Ectoine Synthase

Ectoine synthase (EctC) catalyses the ultimate step of ectoine biosynthesis, a kosmotropic compound produced as compatible solute by many bacteria and some archaea or eukaryotes. EctC is an Fe2+-dependent homodimeric cytoplasmic protein. Using Mössbauer spectroscopy, molecular dynamics simulations and QM/MM calculations, we determined the most likely coordination number and geometry of the Fe2+ ion and proposed a mechanism of the EctC-catalysed reaction. Most notably, we show that apart from the three amino acids binding to the iron ion (Glu57, Tyr84 and His92), one water molecule and one hydroxide ion are required as additional ligands for the reaction to occur. They fill the first coordination sphere of the Fe2+-cofactor and act as critical proton donors and acceptors during the cyclization reaction.

Fluoride Ion Binding and Translocation in the CLCF Fluoride/Proton Antiporter: Molecular Insights from Combined Quantum-Mechanical/Molecular-Mechanical Modeling

CLCF fluoride/proton antiporters move fluoride ions out of bacterial cells, leading to fluoride resistance in these bacteria. However, many details about their operating mechanisms remain unclear. Here, we report a combined quantum-mechanical/molecular-mechanical (QM/MM) study of a CLCF homologue from Enterococci casseliflavus (Eca), in accord with the previously proposed windmill mechanism. Our multiscale modeling sheds light on two critical steps in the transport cycle: (i) the external gating residue E118 pushing a fluoride in the external binding site into the extracellular vestibule and (ii) an incoming fluoride reconquering the external binding site by forcing out E118. Both steps feature competitions for the external binding site between the negatively charged carboxylate of E118 and the fluoride. Remarkably, the displaced E118 by fluoride accepts a proton from the nearby R117, initiating the next transport cycle. We also demonstrate the importance of accurate quantum descriptions of fluoride solvation. Our results provide clues to the mysterious E318 residue near the central binding site, suggesting that the transport activities are unlikely to be disrupted by the glutamate interacting with a well-solvated fluoride at the central binding site. This differs significantly from the structurally similar CLC chloride/proton antiporters, where a fluoride trapped deep in the hydrophobic pore causes the transporter to be locked down. A free-energy barrier of 10-15 kcal/mol was estimated via umbrella sampling for a fluoride ion traveling through the pore to repopulate the external binding site.

Performance of point charge embedding schemes for excited states in molecular organic crystals

Modeling excited state processes in molecular crystals is relevant for several applications. A popular approach for studying excited state molecular crystals is to use cluster models embedded in point charges. In this paper, we compare the performance of several embedding models in predicting excited states and S1-S0 optical gaps for a set of crystals from the X23 molecular crystal database. The performance of atomic charges based on ground or excited states was examined for cluster models, Ewald embedding, and self-consistent approaches. We investigated the impact of various factors, such as the level of theory, basis sets, embedding models, and the level of localization of the excitation. We consider different levels of theory, including time-dependent density functional theory and Tamm-Dancoff approximation (TDA) (DFT functionals: ωB97X-D and PBE0), CC2, complete active space self-consistent field, and CASPT2. We also explore the impact of selection of the QM region, charge leakage, and level of theory for the description of different kinds of excited states. We implemented three schemes based on distance thresholds to overcome overpolarization and charge leakage in molecular crystals. Our findings are compared against experimental data, G0W0-BSE, periodic TDA, and optimally tuned screened range-separated functionals.

MoBioTools: A toolkit to setup quantum mechanics/molecular mechanics calculations

We present a toolkit that allows for the preparation of QM/MM input files from a conformational ensemble of molecular geometries. The package is currently compatible with trajectory and topology files in Amber, CHARMM, GROMACS and NAMD formats, and has the possibility to generate QM/MM input files for Gaussian (09 and 16), Orca (≥4.0), NWChem and (Open)Molcas. The toolkit can be used in command line, so that no programming experience is required, although it presents some features that can also be employed as a python application programming interface. We apply the toolkit in four situations in which different electronic-structure properties of organic molecules in the presence of a solvent or a complex biological environment are computed: the reduction potential of the nucleobases in acetonitrile, an energy decomposition analysis of tyrosine interacting with water, the absorption spectrum of an azobenzene derivative integrated into a voltage-gated ion channel, and the absorption and emission spectra of the luciferine/luciferase complex. These examples show that the toolkit can be employed in a manifold of situations for both the electronic ground state and electronically excited states. It also allows for the automatic correction of the active space in the case of CASSCF calculations on an ensemble of geometries, as it is shown for the azobenzene derivative photoswitch case.

Quantum computing in pharma: A multilayer embedding approach for near future applications

Quantum computers are special purpose machines that are expected to be particularly useful in simulating strongly correlated chemical systems. The quantum computer excels at treating a moderate number of orbitals within an active space in a fully quantum mechanical manner. We present a quantum phase estimation calculation on F₂ in a (2,2) active space on Rigetti's Aspen-11 QPU. While this is a promising start, it also underlines the need for carefully selecting the orbital spaces treated by the quantum computer. In this work, a scheme for selecting such an active space automatically is described and simulated results obtained using both the quantum phase estimation (QPE) and variational quantum eigensolver (VQE) algorithms are presented and combined with a subtractive method to enable accurate description of the environment. The active occupied space is selected from orbitals localized on the chemically relevant fragment of the molecule, while the corresponding virtual space is chosen based on the magnitude of interactions with the occupied space calculated from perturbation theory. This protocol is then applied to two chemical systems of pharmaceutical relevance: the enzyme [Fe] hydrogenase and the photosenzitizer temoporfin. While the sizes of the active spaces currently amenable to a quantum computational treatment are not enough to demonstrate quantum advantage, the procedure outlined here is applicable to any active space size, including those that are outside the reach of classical computation.

Challenges in Protein QM/MM Simulations with Intra-Backbone Link Atoms

Hybrid quantum mechanical/molecular mechanical (QM/MM) simulations fuel discoveries in many fields of science including computational biochemistry and enzymology. Development of more convenient tools leads to an increase in the number of works in which mechanical insights into enzymes' mode of operation are obtained. Most commonly, these tools feature hydrogen-capping (link atom) approach to provide coupling between QM and MM subsystems across a covalent bond. Extensive studies were conducted to provide a solid foundation for the correctness of such an approach when a bond to a nonpolar MM atom is considered. However, not every task may be accomplished this way. Certain scenarios of using QM/MM in computational enzymology encourage or even necessitate the incorporation of backbone atoms into the QM region. Two out of three backbone atoms are polar, and in QM/MM with electrostatic embedding, a neighboring link atom will be hyperpolarized. Several schemes to mitigate this effect were previously proposed alongside a rigorous assessment of quantitative effects on model systems. However, it was not clear whether they may translate into qualitatively different results and how link atom hyperpolarization may manifest itself in a real-life enzymological scenario. Here, we show that the consequences of such an artifact may be severe and may completely overturn the conclusions drawn from the simulations. Our case advocates for the use of charge redistribution schemes whenever intra-backbone QM/MM boundaries are considered. Moreover, we addressed how different boundary types and charge redistribution schemes influence backbone dynamics. We showed that the results are heavily dependent on which boundary MM terms are retained, with charge alteration being of secondary importance. In the worst case, only three intra-backbone boundaries may be used with relative confidence in the adequacy of resulting simulations, irrespective of the hyperpolarization mitigation scheme. Thus, advances in the field are certainly needed to fuel new discoveries. As of now, we believe that issues raised in this work might encourage authors in the field to report what boundaries, boundary MM terms, and charge redistribution schemes they are using, so their results may be correctly interpreted.

Intermolecular Interactions of Nucleoside Antibiotic Tunicamycin with On-Target MraYCB-TUN and Off-Target DPAGT1-TUN in the Active Sites Delineated by Quantum Mechanics/Molecular Mechanics Calculations

Tunicamycin (TUN) is a nucleoside antibiotic with a complex structure comprising uracil, tunicamine sugar, N-acetylglucosamine (GlcNAc), and fatty acyl tail moieties. TUN, known as a canonical inhibitor, blocks vital functions of certain transmembrane protein families, for example, the insect enzyme dolichyl phosphate α-N-acetylglucosaminylphosphotransferase (DPAGT1) of Spodoptera frugiperda and the bacterial enzyme phospho-N-acetylmuramoylpentapeptide translocase (MraY_CB) of Clostridium bolteae. Accurate description of protein-drug interactions has an immense impact on structure-based drug design, while the main challenge is to create proper topology and parameter entries for TUN in modeling protein-TUN interactions given the structural complexity. Starting from DPAGT1-TUN and MraY_CB-TUN crystal structures, we first sketched these structural complexes on the basis of the CHARMM36 force field and optimized each of them using quantum mechanics/molecular mechanics (QM/MM) calculations. By continuing calculations on the active site (QM region) of each optimized structure, we specified the characteristics of intermolecular interactions contributing to the binding of TUN to each active site by quantum theory of atoms in molecules (QTAIM) and natural bond orbital (NBO) analyses at the M06-2X/6-31G** level. The results outlined that TUN insertion into each active site requires multiple weak, moderate, and strong hydrogen bonds accompanying charge-dipole, dipole-dipole, and hydrophobic interactions among different TUN moieties and adjacent residues. The water-mediated interactions also play central roles in situating the uracil and tunicamine moieties of TUN within the DPAGT1 active site as well as in preserving the uracil-binding pocket in the MraY_CB active site. The TUN binds more strongly to DPAGT1 than to MraY_CB. The information garnered here is valuable particularly for better understanding mode of action at the molecular level, as it is conducive to developing next generations of nucleoside antibiotics.

INAQS, a Generic Interface for Nonadiabatic QM/MM Dynamics: Design, Implementation, and Validation for GROMACS/Q-CHEM simulations

The accurate description of large molecular systems in complex environments remains an ongoing challenge for the field of computational chemistry. This problem is even more pronounced for photoinduced processes, as multiple excited electronic states and their corresponding nonadiabatic couplings must be taken into account. Multiscale approaches such as hybrid quantum mechanics/molecular mechanics (QM/MM) offer a balanced compromise between accuracy and computational burden. Here, we introduce an open-source software package (INAQS) for nonadiabatic QM/MM simulations that bridges the sampling capabilities of the GROMACS MD package and the excited-state infrastructure of the Q-CHEM electronic structure software. The interface is simple and can be adapted easily to other MD codes. The code supports a variety of different trajectory-based molecular dynamics, ranging from Born-Oppenheimer to surface hopping dynamics. To illustrate the power of this combination, we simulate electronic absorption spectra, free-energy surfaces along a reaction coordinate, and the excited-state dynamics of 1,3-cyclohexadiene in solution.

How Reproducible Are QM/MM Simulations? Lessons from Computational Studies of the Covalent Inhibition of the SARS-CoV-2 Main Protease by Carmofur

This work explores the level of transparency in reporting the details of computational protocols that is required for practical reproducibility of quantum mechanics/molecular mechanics (QM/MM) simulations. Using the reaction of an essential SARS-CoV-2 enzyme (the main protease) with a covalent inhibitor (carmofur) as a test case of chemical reactions in biomolecules, we carried out QM/MM calculations to determine the structures and energies of the reactants, the product, and the transition state/intermediate using analogous QM/MM models implemented in two software packages, NWChem and Q-Chem. Our main benchmarking goal was to reproduce the key energetics computed with the two packages. Our results indicate that quantitative agreement (within the numerical thresholds used in calculations) is difficult to achieve. We show that rather minor details of QM/MM simulations must be reported in order to ensure the reproducibility of the results and offer suggestions toward developing practical guidelines for reporting the results of biosimulations.

Role of the CarH photoreceptor protein environment in the modulation of cobalamin photochemistry

The photochemistry of cobalamins has recently been found to have biological importance, with the discovery of bacterial photoreceptor proteins, such as CarH and AerR. CarH and AerR, are involved in the light regulation of carotenoid biosynthesis and bacteriochlorophyll biosynthesis, respectively, in bacteria. Experimental transient absorption spectroscopic studies have indicated unusual photochemical behavior of 5'-deoxy-5'-adenosylcobalamin (AdoCbl) in CarH, with excited-state charge separation between cobalt and adenosyl and possible heterolytic cleavage of the Co-adenosyl bond, as opposed to the homolytic cleavage observed in aqueous solution and in many AdoCbl-based enzymes. We employ molecular dynamics and hybrid quantum mechanical/molecular mechanical calculations to obtain a microscopic understanding of the modulation of the excited electronic states of AdoCbl by the CarH protein environment, in contrast to aqueous solution and AdoCbl-based enzymes. Our results indicate a progressive stabilization of the electronic states involving charge transfer (CT) from cobalt/corrin to adenine on changing the environment from gas phase to water to solvated CarH. The solvent exposure of the adenosyl ligand in CarH, the π-stacking interaction between a tryptophan and the adenine moiety, and the hydrogen-bonding interaction between a glutamate and the lower axial ligand of cobalt are found to contribute to the stabilization of the states involving CT to adenine. The combination of these three factors, the latter two of which can be experimentally tested via mutagenesis studies, is absent in an aqueous solvent environment and in AdoCbl-based enzymes. The favored CT from metal and/or corrin to adenine in CarH may promote heterolytic cleavage of the cobalt-adenosyl bond proposed by experimental studies. Overall, this work provides novel, to our knowledge, physical insights into the mechanism of CarH function and directions for future experimental investigations. The fundamental understanding of the mechanism of CarH functioning will serve the development of optogenetic tools based on the new class of B₁₂-dependent photoreceptors.

Spin-projected QM/MM Free Energy Simulations for Oxidation Reaction of Guanine in B-DNA by Singlet Oxygen

Guanine is the most susceptible base to oxidation damage induced by reactive oxygen species including singlet oxygen (¹ O₂ , ¹ Δ_g ). We clarify whether the first step of guanine oxidation in B-DNA proceeds via either a zwitterionic or a diradical intermediate. The free energy profiles are calculated by means of a combined quantum mechanical and molecular mechanical (QM/MM) method coupled with the adaptive biasing force (ABF) method. To describe the open-shell electronic structure of ¹ O₂ correctly, the broken-symmetry spin-unrestricted density functional theory (BS-UDFT) with an approximate spin projection (AP) correction is applied to the QM region. We find that the effect of spin contamination on the activation and reaction free energies is up to ∼8 kcal mol^-1 , which is too large to be neglected. The QM(AP-ULC-BLYP)/MM-based free energy calculations also reveal that the reaction proceeds through a diradical transition state, followed by a conversion to a zwitterionic intermediate. Our computed activation energy of 5.2 kcal mol^-1 matches experimentally observed range (0∼6 kcal mol^-1 ).

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 continuum models provide an effective electrostatic embedding model for fragment-based chemical shift prediction in challenging systems

Ab initio nuclear magnetic resonance chemical shift prediction provides an important tool for interpreting and assigning experimental spectra, but it becomes computationally prohibitive in large systems. The computational costs can be reduced considerably by fragmentation of the large system into a series of contributions from many smaller subsystems. However, the presence of charged functional groups and the need to partition the system across covalent bonds create complications in biomolecules that typically require the use of large fragments and careful descriptions of the electrostatic environment. The present work shows how a model that combines chemical shielding contributions from non-overlapping monomer and dimer fragments embedded in a polarizable continuum model provides a simple, easy-to-implement, and computationally inexpensive approach for predicting chemical shifts in complex systems. The model's performance proves rather insensitive to the continuum dielectric constant, making the selection of the optimal embedding dielectric less critical. The PCM-embedded fragment model is demonstrated to perform well across systems ranging from molecular crystals to proteins.

Electronic Structure Modeling of Metal-Organic Frameworks

Owing to their molecular building blocks, yet highly crystalline nature, metal-organic frameworks (MOFs) sit at the interface between molecule and material. Their diverse structures and compositions enable them to be useful materials as catalysts in heterogeneous reactions, electrical conductors in energy storage and transfer applications, chromophores in photoenabled chemical transformations, and beyond. In all cases, density functional theory (DFT) and higher-level methods for electronic structure determination provide valuable quantitative information about the electronic properties that underpin the functions of these frameworks. However, there are only two general modeling approaches in conventional electronic structure software packages: those that treat materials as extended, periodic solids, and those that treat materials as discrete molecules. Each approach has features and benefits; both have been widely employed to understand the emergent chemistry that arises from the formation of the metal-organic interface. This Review canvases these approaches to date, with emphasis placed on the application of electronic structure theory to explore reactivity and electron transfer using periodic, molecular, and embedded models. This includes (i) computational chemistry considerations such as how functional, k-grid, and other model variables are selected to enable insights into MOF properties, (ii) extended solid models that treat MOFs as materials rather than molecules, (iii) the mechanics of cluster extraction and subsequent chemistry enabled by these molecular models, (iv) catalytic studies using both solids and clusters thereof, and (v) embedded, mixed-method approaches, which simulate a fraction of the material using one level of theory and the remainder of the material using another dissimilar theoretical implementation.

Clustered DNA Damage: Electronic Properties and Their Influence on Charge Transfer. 7,8-Dihydro-8-Oxo-2'-Deoxyguaosine Versus 5',8-Cyclo-2'-Deoxyadenosines: A Theoretical Approach

Approximately 3 × 10¹⁷ DNA damage events take place per hour in the human body. Within clustered DNA lesions, they pose a serious problem for repair proteins, especially for iron–sulfur glycosylases (MutyH), which can recognize them by the electron-transfer process. It has been found that the presence of both 5′,8-cyclo-2′-deoxyadenosine (cdA) diastereomers in the ds-DNA structure, as part of a clustered lesion, can influence vertical radical cation distribution within the proximal part of the double helix, i.e., d[~oxoGcAoxoG~] (7,8-dihydro-8-oxo-2′-deoxyguaosine - ^oxodG). Here, the influence of cdA, “the simplest tandem lesion”, on the charge transfer through ds-DNA was taken into theoretical consideration at the M062x/6-31+G** level of theory in the aqueous phase. It was shown that the presence of (5′S)- or (5′R)-cdA leads to a slowdown in the hole transfer by one order of magnitude between the neighboring dG→^oxodG in comparison to “native” ds-DNA. Therefore, it can be concluded that such clustered lesions can lead to defective damage recognition with a subsequent slowing down of the DNA repair process, giving rise to an increase in mutations. As a result, the unrepaired, ^oxodG: dA base pair prior to genetic information replication can finally result in GC → TA or AT→CG transversion. This type of mutation is commonly observed in human cancer cells. Moreover, because local multiple damage sites (LMSD) are effectively produced as a result of ionization factors, the presented data in this article might be useful in developing a new scheme of radiotherapy treatment against the background of DNA repair efficiency.

The Influence of (5'R)- and (5'S)-5',8-Cyclo-2'-Deoxyadenosine on UDG and hAPE1 Activity. Tandem Lesions are the Base Excision Repair System's Nightmare

DNA lesions are formed continuously in each living cell as a result of environmental factors, ionisation radiation, metabolic processes, etc. Most lesions are removed from the genome by the base excision repair system (BER). The activation of the BER protein cascade starts with DNA damage recognition by glycosylases. Uracil-DNA glycosylase (UDG) is one of the most evolutionary preserved glycosylases which remove the frequently occurring 2′-deoxyuridine from single (ss) and double-stranded (ds) oligonucleotides. Conversely, the unique tandem lesions (5′R)- and (5′S)-5′,8-cyclo-2′-deoxyadenosine (cdA) are not suitable substrates for BER machinery and are released from the genome by the nucleotide excision repair (NER) system. However, the cyclopurines appearing in a clustered DNA damage structure can influence the BER process of other lesions like dU. In this article, UDG inhibition by 5′S- and 5′R-cdA is shown and discussed in an experimental and theoretical manner. This phenomenon was observed when a tandem lesion appears in single or double-stranded oligonucleotides next to dU, on its 3′-end side. The cdA shift to the 5′-end side of dU in ss-DNA stops this effect in both cdA diastereomers. Surprisingly, in the case of ds-DNA, 5′S-cdA completely blocks uracil excision by UDG. Conversely, 5′R-cdA allows glycosylase for uracil removal, but the subsequently formed apurinic/apyrimidinic (AP) site is not suitable for human AP-site endonuclease 1 (hAPE1) activity. In conclusion, the appearance of the discussed tandem lesion in the structure of single or double-stranded DNA can stop the entire base repair process at its beginning, which due to UDG and hAPE1 inhibition can lead to mutagenesis. On the other hand, the presented results can cast some light on the UDG or hAPE1 inhibitors being used as a potential treatment.

Multilink F* Method for Combined Quantum Mechanical and Molecular Mechanical Calculations of Complex Systems

Combined quantum mechanical and molecular mechanical (QM/MM) studies on catalysis in metal-organic frameworks (MOFs) are relatively undeveloped in contrast to the wide use of QM/MM for enzyme catalysis. One reason is that the currently available methods for treating QM-MM boundaries are not fully compatible with the combination of features in MOFs, namely, their high connectivity, their polar bonds (e.g., metal-oxygen bonds), and their potential boundary atoms with high partial atomic charges. The treatment of polar bonds can be improved by using tuned link atoms, but both the widely used H link atom method and the F* link atom method provide limited options in placing the QM-MM boundary in MOFs and other covalently bonded solids, which seriously reduces the efficiency of QM/MM calculations. Here, we propose a generalized version of the F* link atom method with greater flexibility for the placement of the QM-MM boundary in MOFs and with a practical scheme for tuning. The new method, called the multilink F* method, allows a large part of an inorganic node of a MOF to be partitioned into the MM subsystem to increase the efficiency. Our validation calculations on dimerization of ethylene to 1-butene by a nickel catalyst supported on a MOF show that the overall performance of QM/MM calculations with the multilink F* method is excellent for energies, geometries, and partial atomic charges.

Surface Hopping within an Exciton Picture. An Electrostatic Embedding Scheme

We report the development and the implementation of an exciton approach that allows ab initio nonadiabatic dynamics simulations of electronic excitation energy transfer in multichromophoric systems. For the dynamics, a trajectory-based strategy is used within the surface hopping formulation. The approach features a consistent hybrid formulation that allows the construction of potential energy surfaces and gradients by combining quantum mechanics and molecular mechanics within an electrostatic embedding scheme. As an application, the study of a molecular dyad consisting of a covalently bound BODIPY moiety and a tetrathiophene group is presented using time-dependent density functional theory (TDDFT). The results obtained with the exciton model are compared to previously performed full TDDFT dynamics of the same system. Our results show excellent agreement with the full TDDFT results, indicating that the couplings that lead to excitation energy transfer (EET) are dominated by Coulomb interaction terms and that charge-transfer states are not necessary to properly describe the nonadiabatic dynamics of the system. The exciton model also reveals ultrafast coherent oscillations of the excitation between the two units in the dyad, which occur during the first 50 fs.

The AT Interstrand Cross-Link: Structure, Electronic Properties, and Influence on Charge Transfer in dsDNA

The interaction of chemical and physical agents with genetic material can lead to almost 80 different DNA damage formations. The targeted intentional DNA damage by radiotherapy or chemotherapy is a front-line anticancer therapy. An interstrand cross-link can result from ionization radiation or specific chemical agents, such as trans-/cisplatin activity. Here, the influence of the adenine and thymidine (AT) interstrand linkage, the covalent bond between the adenine N6 and thymidine C5 methylene group, on the isolated base pair as well as double-stranded DNA (dsDNA) was taken into quantum mechanical/molecular mechanical (QM/MM) consideration at the m062x/6-31+G*:UFF level of theory in the aqueous phase. All the results presented in this article, for the first time, show that an AT-interstrand cross-link (ICL) changes the positive and negative charge migration process due to a higher activation energy forced by the cross-link’s presence. However, the final radical cation destination in cross-linked DNA is left in the same place as in a native double-stranded-deoxyoligonucleotide. Additionally, the direction of the radical anion transfer was found to be opposite to that of native dsDNA. Therefore, it can be postulated that the appearance of the AT-ICL does not disturb the hole migration in the double helix, with subsequent effective changes in the electron migration process.

Parametrization of Combined Quantum Mechanical and Molecular Mechanical Methods: Bond-Tuned Link Atoms

Combined quantum mechanical and molecular mechanical (QM/MM) methods are the most powerful available methods for high-level treatments of subsystems of very large systems. The treatment of the QM−MM boundary strongly affects the accuracy of QM/MM calculations. For QM/MM calculations having covalent bonds cut by the QM−MM boundary, it has been proposed previously to use a scheme with system-specific tuned fluorine link atoms. Here, we propose a broadly parametrized scheme where the parameters of the tuned F link atoms depend only on the type of bond being cut. In the proposed new scheme, the F link atom is tuned for systems with a certain type of cut bond at the QM−MM boundary instead of for a specific target system, and the resulting link atoms are call bond-tuned link atoms. In principle, the bond-tuned link atoms can be as convenient as the popular H link atoms, and they are especially well adapted for high-throughput and accurate QM/MM calculations. Here, we present the parameters for several kinds of cut bonds along with a set of validation calculations that confirm that the proposed bond-tuned link-atom scheme can be as accurate as the system-specific tuned F link-atom scheme.

NAMD goes quantum: an integrative suite for hybrid simulations

Hybrid methods that combine quantum mechanics (QM) and molecular mechanics (MM) can be applied to studies of reaction mechanisms in locations ranging from active sites of small enzymes to multiple sites in large bioenergetic complexes. By combining the widely used molecular dynamics and visualization programs NAMD and VMD with the quantum chemistry packages ORCA and MOPAC, we created an integrated, comprehensive, customizable, and easy-to-use suite (http://www.ks.uiuc.edu/Research/qmmm). Through the QwikMD interface, setup, execution, visualization, and analysis are streamlined for all levels of expertise.

Chloride Ion Transport by the E. coli CLC Cl-/H+ Antiporter: A Combined Quantum-Mechanical and Molecular-Mechanical Study

We performed steered molecular dynamics (SMD) and umbrella sampling simulations of Cl^- ion migration through the transmembrane domain of a prototypical E. coli CLC Cl^-/H⁺ antiporter by employing combined quantum-mechanical (QM) and molecular-mechanical (MM) calculations. The SMD simulations revealed interesting conformational changes of the protein. While no large-amplitude motions of the protein were observed during pore opening, the side chain rotation of the protonated external gating residue Glu148 was found to be critical for full access of the channel entrance by Cl^-. Moving the anion into the external binding site (S_ext) induced small-amplitude shifting of the protein backbone at the N-terminal end of helix F. As Cl^- traveled through the pore, rigid-body swinging motions of helix R separated it from helix D. Helix R returned to its original position once Cl^- exited the channel. Population analysis based on polarized wavefunction from QM/MM calculations discovered significant (up to 20%) charge loss for Cl^- along the ion translocation pathway inside the pore. The delocalized charge was redistributed onto the pore residues, especially the functional groups containing π bonds (e.g., the Tyr445 side chain), while the charges of the H atoms coordinating Cl^- changed almost negligibly. Potentials of mean force computed from umbrella sampling at the QM/MM and MM levels both displayed barriers at the same locations near the pore entrance and exit. However, the QM/MM PMF showed higher barriers (~10 kcal/mol) than the MM PMF (~2 kcal/mol). Binding energy calculations indicated that the interactions between Cl^- and certain pore residues were overestimated by the semi-empirical PM3 Hamiltonian and underestimated by the CHARMM36 force fields, both of which were employed in the umbrella sampling simulations. In particular, CHARMM36 underestimated binding interactions for the functional groups containing π bonds, missing the stabilizations of the Cl^- ion due to electron delocalization. The results suggested that it is important to explore these quantum effects for accurate descriptions of the Cl^- transport.

Combined quantum mechanical and molecular mechanical method for metal–organic frameworks: proton topologies of NU-1000

A force field is presented for NU-1000 and tested by electronically embedded QM/MM calculations, which yield accurate structures and relative energies for various proton topologies.

Molecular Mechanism of Flavin Photoprotection by Archaeal Dodecin: Photoinduced Electron Transfer and Mg2+-Promoted Proton Transfer

Photoinduced biochemical reactions are ubiquitously governed by derivatives of flavin, which is a key player in a manifold of cellular redox reactions. The photoreactivity of flavins is also one of their greatest disadvantages as the molecules are sensitive to photodegradation. To prevent this unfavorable reaction, UV-light-exposed archaea bacteria, such as Halobacterium salinarum, manage the task of protecting flavin derivatives by dodecin, a protein which stores flavins and efficiently photoprotects them. In this study, we shed light on the photoprotection mechanism, i.e., the excited state quenching mechanism by dodecin using computational methodology. Molecular dynamics (MD) simulations unraveled the hydrogen bond network in the flavin binding pocket as a starting point for proton transfer upon preceding electron transfer. Using high-level ab initio quantum chemical methods, different proton transfer channels have been investigated and an energetically feasible Mg²⁺-promoted channel has been identified fully explaining previous experimental observations. This is the first extensive theoretical study of archaeal dodecin, furthering the understanding of its photocycle and manipulation.

Quantum mechanical force fields for condensed phase molecular simulations

Molecular simulations are powerful tools for providing atomic-level details into complex chemical and physical processes that occur in the condensed phase. For strongly interacting systems where quantum many-body effects are known to play an important role, density-functional methods are often used to provide the model with the potential energy used to drive dynamics. These methods, however, suffer from two major drawbacks. First, they are often too computationally intensive to practically apply to large systems over long time scales, limiting their scope of application. Second, there remain challenges for these models to obtain the necessary level of accuracy for weak non-bonded interactions to obtain quantitative accuracy for a wide range of condensed phase properties. Quantum mechanical force fields (QMFFs) provide a potential solution to both of these limitations. In this review, we address recent advances in the development of QMFFs for condensed phase simulations. In particular, we examine the development of QMFF models using both approximate and ab initio density-functional models, the treatment of short-ranged non-bonded and long-ranged electrostatic interactions, and stability issues in molecular dynamics calculations. Example calculations are provided for crystalline systems, liquid water, and ionic liquids. We conclude with a perspective for emerging challenges and future research directions.

A quantum-mechanics molecular-mechanics scheme for extended systems

We introduce and discuss a hybrid quantum-mechanics molecular-mechanics (QM-MM) approach for Car-Parrinello DFT simulations with pseudopotentials and planewaves basis, designed for the treatment of periodic systems. In this implementation the MM atoms are considered as additional QM ions having fractional charges of either sign, which provides conceptual and computational simplicity by exploiting the machinery already existing in planewave codes to deal with electrostatics in periodic boundary conditions. With this strategy, both the QM and MM regions are contained in the same supercell, which determines the periodicity for the whole system. Thus, while this method is not meant to compete with non-periodic QM-MM schemes able to handle extremely large but finite MM regions, it is shown that for periodic systems of a few hundred atoms, our approach provides substantial savings in computational times by treating classically a fraction of the particles. The performance and accuracy of the method is assessed through the study of energetic, structural, and dynamical aspects of the water dimer and of the aqueous bulk phase. Finally, the QM-MM scheme is applied to the computation of the vibrational spectra of water layers adsorbed at the TiO2 anatase (1 0 1) solid-liquid interface. This investigation suggests that the inclusion of a second monolayer of H2O molecules is sufficient to induce on the first adsorbed layer, a vibrational dynamics similar to that taking place in the presence of an aqueous environment. The present QM-MM scheme appears as a very interesting tool to efficiently perform molecular dynamics simulations of complex condensed matter systems, from solutions to nanoconfined fluids to different kind of interfaces.

The MOD-QM/MM Method: Applications to Studies of Photosystem II and DNA G-Quadruplexes

Quantum mechanics/molecular mechanics (QM/MM) hybrid methods are currently the most powerful computational tools for studies of structure/function relations and catalytic sites embedded in macrobiomolecules (eg, proteins and nucleic acids). QM/MM methodologies are highly efficient since they implement quantum chemistry methods for modeling only the portion of the system involving bond-breaking/forming processes (QM layer), as influenced by the surrounding molecular environment described in terms of molecular mechanics force fields (MM layer). Some of the limitations of QM/MM methods when polarization effects are not explicitly considered include the approximate treatment of electrostatic interactions between QM and MM layers. Here, we review recent advances in the development of computational protocols that allow for rigorous modeling of electrostatic interactions in biomacromolecules and structural refinement, beyond the common limitations of QM/MM hybrid methods. We focus on photosystem II (PSII) with emphasis on the description of the oxygen-evolving complex (OEC) and its high-resolution extended X-ray absorption fine structure spectra (EXAFS) in conjunction with Monte Carlo structural refinement. Furthermore, we review QM/MM structural refinement studies of DNA G4 quadruplexes with embedded monovalent cations and direct comparisons to NMR data.

Simple Link Atom Saccharide Hybrid (SLASH) Treatment for Glycosidic Bonds at the QM/MM Boundary

We investigated link atom approaches for treating the polar C-O bond with particular reference to the glycosidic bond found in complex carbohydrates. We show that cutting this bond after the oxygen in the QM region and saturating the QM system with a hydrogen link atom leads to greater conformational and configurational accuracy at the boundary compared with cutting the bond before oxygen and saturating the QM system with a halogen link atom to represent the oxygen. Furthermore, we find that balancing the MM atom charges and redistributing the boundary atom charges at the QM/MM boundary minimizes the effect of the link atom, both energetically and structurally. This is illustrated via a series of calculations on a set of carbohydrate and carbohydrate-like model compounds. Finally, we confirm the validity of our model by performing molecular dynamics simulations for a typical disaccharide model compound in water. Our postsimulation conformational and configurational analyses show that the oxygen-to-water hydrogen pair distribution functions and the Φ,Ψ distributions at the glycosidic boundary between the quantum and classical regions compare favorably with results obtained from complete QM and complete MM treatments of the saccharide.

Combined Electrostatically Embedded Multiconfiguration Molecular Mechanics and Molecular Mechanical Method: Application to Molecular Dynamics Simulation of a Chemical Reaction in Aqueous Solution with Hybrid Density Functional Theory

We here combine the electrostatically embedded multiconfiguration molecular mechanics (EE-MCMM) method for generating global potential energy surfaces in the presence of an electrostatic potential with molecular mechanics (MM). The resulting EE-MCMM/MM method is illustrated by applying it to carry out a molecular dynamics simulation for the symmetric bimolecular reaction Cl(-) + CH3Cl' → ClCH3 + Cl'(-) in aqueous solution with hybrid density functional theory as the quantum mechanical level. The potential of mean force is calculated, and the free energy barrier is found to be 25.3 kcal/mol, which is in good agreement with previous work. The advantage of the combined EE-MCMM and MM method is that the number of quantum mechanical calculations required for the active subsystem is very small compared to straight direct dynamics.

The pDynamo Program for Molecular Simulations using Hybrid Quantum Chemical and Molecular Mechanical Potentials

The pDynamo program has been developed for the simulation of molecular systems using hybrid quantum chemical (QC) and molecular mechanical (MM) potentials. pDynamo is written in a mixture of the computer languages Python and C and is a successor to the previous version of Dynamo, now denoted fDynamo, that was written in Fortran 90 (J. Comput. Chem. 2000, 21, 1088). The current version of Dynamo has a similar range of functionality to the older one but extends it in some significant ways, including the addition of a density functional theory QC capability. This paper gives a general description of pDynamo and outlines some of the advantages and disadvantages that have been encountered in switching computer languages. Some technical aspects of the implementation of pDynamo's algorithms are also discussed and illustrated with the results of example calculations. pDynamo is available on the Web at the address http://www.pdynamo.org and is released under the CeCILL license which is equivalent to the GNU general public license but conforms to the principles of French law.

Combining Quantum Mechanics Methods with Molecular Mechanics Methods in ONIOM

The purpose of this paper is 2-fold. First, we present several extensions to the ONIOM(QM:MM) scheme. In its original formulation, the electrostatic interaction between the regions is included at the classical level. Here we present the extension to electronic embedding. We show how the behavior of ONIOM with electronic embedding can be more stable than QM/MM with electronic embedding. We also investigate the link atom correction, which is implicit in ONIOM but not in QM/MM. Second, we demonstrate some of the practical aspects of ONIOM(QM:MM) calculations. Specifically, we show that the potential surface can be discontinuous when there is bond breaking and forming closer than three bonds from the MM region.

Theoretical Investigation of the Geometries and UV-vis Spectra of Poly(l-glutamic acid) Featuring a Photochromic Azobenzene Side Chain

The geometries and UV-vis spectra of azobenzene dyes grafted as a side chain on poly(l-glutamic acid) have been investigated using a combination of quantum mechanics/molecular mechanics (QM/MM) and time-dependent density functional theory (TD-DFT) methods at the TD-PBE0/6-311+G(d,p)//B3LYP/6-311G(d,p):Amber ff99 level of theory. The influence of the secondary structure of the polypeptide on the electronic properties of both the trans and cis conformations of azobenzene dyes has been studied. It turns out that the grafted dyes exhibit a red-shift of the π → π* absorption energies mainly due to the auxochromic shift induced by the peptidic group used to link the chromophoric unit to the polypeptide and that specific interactions between the glutamic side chain and the azobenzene moiety lead to a large blue-shift of the n → π* transition.

Excited States and Photochemistry of Chromophores in the Photoactive Proteins Explored by the Combined Quantum Mechanical and Molecular Mechanical Calculations

A photoactive protein usually contains a unique chromophore that is responsible for the initial photoresponse and functions of the photoactive protein are determined by the interaction between the chromophore and its protein surroundings. The combined quantum mechanical and molecular mechanical (QM/MM) approach is demonstrated to be a very useful tool for exploring structures and functions of a photoactive protein with the chromophore and its protein surroundings treated by the QM and MM methods, respectively. In this review, we summarize the basic formulas of the QM/MM approach and emphasize its applications to excited states and photoreactions of chromophores in rhodopsin protein, photoactive yellow protein, and green fluorescent protein.