Open-Source, Python-Based Redevelopment of the ChemShell Multiscale QM/MM Environment

"On-the-Fly" Nonadiabatic Dynamics Simulation on the Ultrafast Photoisomerization of a Molecular Photoswitch Iminothioindoxyl: An RMS-CASPT2 Investigation

Iminothioindoxyl (ITI) is a new class of photoswitch that exhibits many excellent properties including well-separated absorption bands in the visible region for both conformers, ultrafast Z to E photoisomerization as well as the millisecond reisomerization at room temperature for the E isomer, and switchable ability in both solids and various solvents. However, the underlying ultrafast photoisomerization mechanism at the atomic level remains unclear. In this work, we have employed a combination of high-level RMS-CASPT2-based static electronic structure calculations and nonadiabatic dynamics simulations to investigate the ultrafast photoisomerization dynamics of ITI. Based on the minimum-energy structures, minimum-energy conical intersections, linear interpolation internal coordinate paths, and nonadiabatic dynamics simulations, the overall photoisomerization scenario of ITI upon excitation is established. Upon excitation around 416 nm, the molecule will be excited to the S2 state considering its close energy to the experimentally measured absorption maximum and larger oscillator strength, from which ultrafast decay of S2 to S1 state can take place efficiently with a time constant of 62 fs. However, the photoisomerization is not likely to complete in the S2 state since the dihedral associated with the Z to E isomerization changes little during the relaxation. Upon relaxing to the S1 state, the molecule will decay to the S0 state ultrafast with a time constant of 232 fs. In contrast, the decay of the S1 state is important for the isomerization considering that the dihedral related to the isomerization of the hopping structures is close to 90°. Therefore, the S1/S0 intersection region should be important for the isomerization of ITI. Arriving at the S0 state, the molecule can either go back to the original Z reactant or isomerize to the E products. At the end of the 500 fs simulation time, the E configuration accounts for nearly 37% of the final structures. Moreover, the photoisomerization mechanism is different from the isomerization mechanism in the ground state; i.e., instead of the inversion mechanism in the ground state, the photoisomerization prefers the rotation mechanism. Our results not only agree well with previous experimental studies but also provide some novel insights that could be helpful for future improvements in the performance of the ITI photoswitches.

Non-Aufbau electronic structure in radical enzymes and control of the highly reactive intermediates

Radicals are highly reactive, short-lived chemical species that normally react indiscriminately with biological materials, and yet, nature has evolved thousands of enzymes that employ radicals to catalyze thermodynamically challenging chemistry. How these enzymes harness highly reactive radical intermediates to steer the catalysis to the correct outcome is a topic of intense investigation. Here, the results of detailed QM/MM calculations on archetype radical B12-enzymes are presented that provide new insights into how these enzymes control the reactivity of radicals within their active sites. The catalytic cycle in B12-enzymes is initiated through the formation of the 5'-deoxyadenosyl (Ado˙) moiety, a primary carbon-centred radical, which must migrate up to 8 Å to reach the target substrate to engage in the next step of the catalytic process, a hydrogen atom abstraction. Our calculations reveal that Ado˙ within the protein environment exhibits an unusual non-Aufbau electronic structure in which the singly occupied molecular orbital is lower in energy than the doubly occupied orbitals, an uncommon phenomenon known as SOMO-HOMO inversion (SHI). We find that the magnitude of SHI in the initially formed Ado˙ is larger compared to when the Ado˙ is near the intended substrate, leading to the former being relatively less reactive. The modulation of the SHI originates from Coulombic interactions of a quantum nature between a negative charge on a conserved glutamate residue and the spin on the Ado˙. Our findings support a novel hypothesis that these enzymes utilize this quantum Coulombic effect as a means of maintaining exquisite control over the chemistry of highly reactive radical intermediates in enzyme active sites.

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.

MiMiC: A high-performance framework for multiscale molecular dynamics simulations

MiMiC is a framework for performing multiscale simulations in which loosely coupled external programs describe individual subsystems at different resolutions and levels of theory. To make it highly efficient and flexible, we adopt an interoperable approach based on a multiple-program multiple-data (MPMD) paradigm, serving as an intermediary responsible for fast data exchange and interactions between the subsystems. The main goal of MiMiC is to avoid interfering with the underlying parallelization of the external programs, including the operability on hybrid architectures (e.g., CPU/GPU), and keep their setup and execution as close as possible to the original. At the moment, MiMiC offers an efficient implementation of electrostatic embedding quantum mechanics/molecular mechanics (QM/MM) that has demonstrated unprecedented parallel scaling in simulations of large biomolecules using CPMD and GROMACS as QM and MM engines, respectively. However, as it is designed for high flexibility with general multiscale models in mind, it can be straightforwardly extended beyond QM/MM. In this article, we illustrate the software design and the features of the framework, which make it a compelling choice for multiscale simulations in the upcoming era of exascale high-performance computing.

Going for Gold(-Standard): Attaining Coupled Cluster Accuracy in Oxide-Supported Nanoclusters

The structure of oxide-supported metal nanoclusters plays an essential role in their sharply enhanced catalytic activity over that of bulk metals. Simulations provide the atomic-scale resolution needed to understand these systems. However, the sensitive mix of metal-metal and metal-support interactions, which govern their structure, puts stringent requirements on the method used, requiring calculations beyond standard density functional theory (DFT). The method of choice is coupled cluster theory [specifically CCSD(T)], but its computational cost has so far prevented its application to these systems. In this work, we showcase two approaches to make CCSD(T) accuracy readily achievable in oxide-supported nanoclusters. First, we leverage the SKZCAM protocol to provide the first benchmarks of oxide-supported nanoclusters, revealing that it is specifically metal-metal interactions that are challenging to capture with DFT. Second, we propose a CCSD(T) correction (ΔCC) to the metal-metal interaction errors in DFT, reaching accuracy comparable to that of the SKZCAM protocol at significantly lower cost. This approach forges a path toward studying larger systems at reliable accuracy, which we highlight by identifying a ground-state structure in agreement with experiments for Au20 on MgO, a challenging system where DFT models have yielded conflicting predictions.

Plugin-based interoperability and ecosystem management for the MolSSI Driver Interface Project

The MolSSI Driver Interface (MDI) Project is an effort to simplify and standardize the process of enabling tight interoperability between independently developed code bases and is supported by numerous software packages across the domain of chemical physics. It enables a wide variety of use cases, including quantum mechanics/molecular mechanics, advanced sampling, path integral molecular dynamics, machine learning, ab initio molecular dynamics, etc. We describe two major developments within the MDI Project that provide novel solutions to key interoperability challenges. The first of these is the development of the MDI Plugin System, which allows MDI-supporting libraries to be used as highly modular plugins, with MDI enforcing a standardized application programming interface across plugins. Codes can use these plugins without linking against them during their build process, and end-users can select which plugin(s) they wish to use at runtime. The MDI Plugin System features a sophisticated callback system that allows codes to interact with plugins on a highly granular level and represents a significant advancement toward increased modularity among scientific codes. The second major development is MDI Mechanic, an ecosystem management tool that utilizes Docker containerization to simplify the process of developing, validating, maintaining, and deploying MDI-supporting codes. Additionally, MDI Mechanic provides a framework for launching MDI simulations in which each interoperating code is executed within a separate computational environment. This eliminates the need to compile multiple production codes within a single computational environment, reducing opportunities for dependency conflicts and lowering the barrier to entry for users of MDI-enabled codes.

Environment-Driven Variability in Absolute Band Edge Positions and Work Functions of Reduced Ceria

The absolute band edge positions and work function (Φ) are the key electronic properties of metal oxides that determine their performance in electronic devices and photocatalysis. However, experimental measurements of these properties often show notable variations, and the mechanisms underlying these discrepancies remain inadequately understood. In this work, we focus on ceria (CeO2), a material renowned for its outstanding oxygen storage capacity, and combine theoretical and experimental techniques to demonstrate environmental modifications of its ionization potential (IP) and Φ. Under O-deficient conditions, reduced ceria exhibits a decreased IP and Φ with significant sensitivity to defect distributions. In contrast, the IP and Φ are elevated in O-rich conditions due to the formation of surface peroxide species. Surface adsorbates and impurities can further augment these variabilities under realistic conditions. We rationalize the shifts in energy levels by separating the individual contributions from bulk and surface factors, using hybrid quantum mechanical/molecular mechanical (QM/MM) embedded-cluster and periodic density functional theory (DFT) calculations supported by interatomic-potential-based electrostatic analyses. Our results highlight the critical role of on-site electrostatic potentials in determining the absolute energy levels in metal oxides, implying a dynamic evolution of band edges under catalytic conditions.

Multiscale biomolecular simulations in the exascale era

The complexity of biological systems and processes, spanning molecular to macroscopic scales, necessitates the use of multiscale simulations to get a comprehensive understanding. Quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations are crucial for capturing processes beyond the reach of classical MD simulations. The advent of exascale computing offers unprecedented opportunities for scientific exploration, not least within life sciences, where simulations are essential to unravel intricate molecular mechanisms underlying biological processes. However, leveraging the immense computational power of exascale computing requires innovative algorithms and software designs. In this context, we discuss the current status and future prospects of multiscale biomolecular simulations on exascale supercomputers with a focus on QM/MM MD. We highlight our own efforts in developing a versatile and high-performance multiscale simulation framework with the aim of efficient utilization of state-of-the-art supercomputers. We showcase its application in uncovering complex biological mechanisms and its potential for leveraging exascale computing.

Interoperable workflows by exchanging grid-based data between quantum-chemical program packages

Quantum-chemical subsystem and embedding methods require complex workflows that may involve multiple quantum-chemical program packages. Moreover, such workflows require the exchange of voluminous data that go beyond simple quantities, such as molecular structures and energies. Here, we describe our approach for addressing this interoperability challenge by exchanging electron densities and embedding potentials as grid-based data. We describe the approach that we have implemented to this end in a dedicated code, PyEmbed, currently part of a Python scripting framework. We discuss how it has facilitated the development of quantum-chemical subsystem and embedding methods and highlight several applications that have been enabled by PyEmbed, including wave-function theory (WFT) in density-functional theory (DFT) embedding schemes mixing non-relativistic and relativistic electronic structure methods, real-time time-dependent DFT-in-DFT approaches, the density-based many-body expansion, and workflows including real-space data analysis and visualization. Our approach demonstrates, in particular, the merits of exchanging (complex) grid-based data and, in general, the potential of modular software development in quantum chemistry, which hinges upon libraries that facilitate interoperability.

Systematic QM/MM Study for Predicting 31P NMR Chemical Shifts of Adenosine Nucleotides in Solution and Stages of ATP Hydrolysis in a Protein Environment

NMR (nuclear magnetic resonance) spectroscopy allows for important atomistic insights into the structure and dynamics of biological macromolecules; however, reliable assignments of experimental spectra are often difficult. Herein, quantum mechanical/molecular mechanical (QM/MM) calculations can provide crucial support. A major problem for the simulations is that experimental NMR signals are time-averaged over much longer time scales, and since computed chemical shifts are highly sensitive to local changes in the electronic and structural environment, sufficiently large averages over representative structural ensembles are essential. This entails high computational demands for reliable simulations. For NMR measurements in biological systems, a nucleus of major interest is 31P since it is both highly present (e.g., in nucleic acids) and easily observable. The focus of our present study is to develop a robust and computationally cost-efficient framework for simulating 31P NMR chemical shifts of nucleotides. We apply this scheme to study the different stages of the ATP hydrolysis reaction catalyzed by p97. Our methodology is based on MM molecular dynamics (MM-MD) sampling, followed by QM/MM structure optimizations and NMR calculations. Overall, our study is one of the most comprehensive QM-based 31P studies in a protein environment and the first to provide computed NMR chemical shifts for multiple nucleotide states in a protein environment. This study sheds light on a process that is challenging to probe experimentally and aims to bridge the gap between measured and calculated NMR spectroscopic properties.

Level-Shifted Embedded Cluster Method for Modeling the Chemistry of Metal Oxides

The embedded cluster method has been used extensively in the study of the chemical and physical properties of metal oxides. This method has been a popular tool due to its relatively high accuracy and low computational cost. An even more promising option may entail integrating the embedded cluster method with the combined quantum mechanical and molecular mechanical (QM/MM) approach, thereby enabling further consideration of interactions within the entire system for superior results. We aim to accurately model the chemistry of metal oxides using this combined scheme. Here, using the prototypical MgO(100) surface as a test system, with Mg9O14 as the cluster in the quantum mechanical region, we show that the embedded cluster with untailored boundary effective core potentials (ECPs) can have frontier orbital energy levels that substantially deviate from the quantum mechanical reference results. This occurs even when Mg9O9, which retains the stoichiometry of MgO, is used as the cluster in the quantum mechanical region. As a result, the chemical properties of the embedded cluster models differ from those of the quantum mechanical reference model. To address this issue, we propose a new variant of the embedded cluster method called the level-shifted embedded cluster (LSEC) method, which allows the energy levels to be shifted to match the reference levels by tuning the boundary ECPs. Our validation calculations on the adsorption of various adsorbates with different properties on the MgO(100) surface show that the overall performance of QM/MM with the LSEC method is excellent for the adsorption energies, geometries, and charge properties. The excellent performance holds for both the nonstoichiometric and stoichiometric clusters (i.e., Mg9O14 and Mg9O9, respectively), demonstrating the robustness of the LSEC method. We expect that the LSEC method can be combined with QM/MM or used separately for future chemical studies of metal oxides and other ionically bonded systems.

Improved Sampling of Adaptive Path Collective Variables by Stabilized Extended-System Dynamics

Because of the complicated multistep nature of many biocatalytic reactions, an a priori definition of reaction coordinates is difficult. Therefore, we apply enhanced sampling algorithms along with adaptive path collective variables (PCVs), which converge to the minimum free energy path (MFEP) during the simulation. We show how PCVs can be combined with the highly efficient well-tempered metadynamics extended-system adaptive biasing force (WTM-eABF) hybrid sampling algorithm, offering dramatically increased sampling efficiency due to its fast adaptation to path updates. For this purpose, we address discontinuities of PCVs that can arise due to path shortcutting or path updates with a novel stabilization algorithm for extended-system methods. In addition, we show how the convergence of simulations can be further accelerated by utilizing the multistate Bennett's acceptance ratio (MBAR) estimator. These methods are applied to the first step of the enzymatic reaction mechanism of pseudouridine synthases, where the ability of path WTM-eABF to efficiently explore intricate molecular transitions is demonstrated.

Seamless Multilayer─A Novel Total Energy Partition Scheme for Embedded and Hybrid Calculations

In this paper, we provide general formulation of a multilayer approach, covering both additive and subtractive quantum mechanics/molecular mechanics (QM/MM) as special cases. After that, we suggest a novel definition of QM/MM total energy based on the consideration of a system divided into three layers. In a simplified form, it is E = E Q M ( 1 + 2 ) - E Q M ( 2 ) + E M M ( 2 + 3 ) , where layers 1, 2, and 3 represent inner QM, outer QM, and classical MM regions, respectively. The novel formulation is also not limited by only QM/MM combination of methods─in fact, any computational methods can be combined in a hybrid calculation. In this paper, we call the new approach seamless multilayer. Test calculations performed for silica and boric oxide show that the new approach requires no QM/MM interface parameterization as well as no or very simple correction terms for boundary atoms. This can greatly facilitate QM/MM studies of covalent inorganic solids. However, test calculations of α-Al2O3 show that for ionic compounds, the new method requires some additional development.

Bulk and Surface Contributions to Ionisation Potentials of Metal Oxides

Determining the absolute band edge positions in solid materials is crucial for optimising their performance in wide-ranging applications including photocatalysis and electronic devices. However, obtaining absolute energies is challenging, as seen in CeO2 , where experimental measurements show substantial discrepancies in the ionisation potential (IP). Here, we have combined several theoretical approaches, from classical electrostatics to quantum mechanics, to elucidate the bulk and surface contributions to the IP of metal oxides. We have determined a theoretical bulk contribution to the IP of stoichiometric CeO2 of only 5.38 eV, while surface orientation results in intrinsic IP variations ranging from 4.2 eV to 8.2 eV. Highly tuneable IPs were also found in TiO2 , ZrO2 , and HfO2 , in which surface polarisation plays a pivotal role in long-range energy level shifting. Our analysis, in addition to rationalising the observed range of experimental results, provides a firm basis for future interpretations of experimental and computational studies of oxide band structures.

Computational Study of the Solid-State Incorporation of Sn(II) Acetate into Zeolite β

Sn-doped zeolites are potent Lewis acid catalysts for important reactions in the context of green and sustainable chemistry; however, their synthesis can have long reaction times and harsh chemical requirements, presenting an obstacle to scale-up and industrial application. To incorporate Sn into the β zeolite framework, solid-state incorporation (SSI) has recently been demonstrated as a fast and solvent-free synthetic method, with no impairment to the high activity and selectivity associated with Sn-β for its catalytic applications. Here, we report an ab initio computational study that combines periodic density functional theory with high-level embedded-cluster quantum/molecular mechanical (QM/MM) to elucidate the mechanistic steps in the synthetic process. Initially, once the Sn(II) acetate precursor coordinates to the β framework, acetic acid forms via a facile hydrogen transfer from the β framework onto the monodentate acetate ligand, with low kinetic barriers for subsequent dissociation of the ligand from the framework-bound Sn. Ketonization of the dissociated acetic acid can occur over the Lewis acidic Sn(II) site to produce CO2 and acetone with a low kinetic barrier (1.03 eV) compared to a gas-phase process (3.84 eV), helping to explain product distributions in good accordance with experimental analysis. Furthermore, we consider the oxidation of the Sn(II) species to form the Sn(IV) active site in the material by O2- and H2O-mediated mechanisms. The kinetic barrier for oxidation via H2 release is 3.26 eV, while the H2O-mediated dehydrogenation process has a minimum barrier of 1.38 eV, which indicates the possible role of residual H2O in the experimental observations of SSI synthesis. However, we find that dehydrogenation is facilitated more significantly by the presence of dioxygen (O2), introduced in the compressed air gas feed, via a two-step process oxidation process that forms H2O2 as an intermediate and has greatly reduced kinetic barriers of 0.25 and 0.26 eV. The results provide insight into how Sn insertion into β occurs during SSI and demonstrate the possible mechanism of top-down synthetic procedures for metal insertion into zeolites.

Characterization and structure-based protein engineering of a regiospecific saponin acetyltransferase from Astragalus membranaceus

Acetylation contributes to the bioactivity of numerous medicinally important natural products. However, little is known about the acetylation on sugar moieties. Here we report a saponin acetyltransferase from Astragalus membranaceus. AmAT7-3 is discovered through a stepwise gene mining approach and characterized as the xylose C3'/C4'-O-acetyltransferse of astragaloside IV (1). To elucidate its catalytic mechanism, complex crystal structures of AmAT7-3/1 and AmAT7-3A310G/1 are obtained, which reveal a large active pocket decided by a specific sequence AADAG. Combining with QM/MM computation, the regiospecificity of AmAT7-3 is determined by sugar positioning modulated by surrounding amino acids including #A310 and #L290. Furthermore, a small mutant library is built using semi-rational design, where variants A310G and A310W are found to catalyze specific C3'-O and C4'-O acetylation, respectively. AmAT7-3 and its variants are also employed to acetylate other bioactive saponins. This work expands the understanding of saponin acetyltransferases, and provide efficient catalytic tools for saponin acetylation.

Multiscale QM/MM modelling of catalytic systems with ChemShell

Hybrid quantum mechanical/molecular mechanical (QM/MM) methods are a powerful computational tool for the investigation of all forms of catalysis, as they allow for an accurate description of reactions occurring at catalytic sites in the context of a complicated electrostatic environment. The scriptable computational chemistry environment ChemShell is a leading software package for QM/MM calculations, providing a flexible, high performance framework for modelling both biomolecular and materials catalysis. We present an overview of recent applications of ChemShell to problems in catalysis and review new functionality introduced into the redeveloped Python-based version of ChemShell to support catalytic modelling. These include a fully guided workflow for biomolecular QM/MM modelling, starting from an experimental structure, a periodic QM/MM embedding scheme to support modelling of metallic materials, and a comprehensive set of tutorials for biomolecular and materials modelling.

Stability of Single Gold Atoms on Defective and Doped Diamond Surfaces

Polycrystalline boron-doped diamond (BDD) is widely used as a working electrode material in electrochemistry, and its properties, such as its stability, make it an appealing support material for nanostructures in electrocatalytic applications. Recent experiments have shown that electrodeposition can lead to the creation of stable small nanoclusters and even single gold adatoms on the BDD surfaces. We investigate the adsorption energy and kinetic stability of single gold atoms adsorbed onto an atomistic model of BDD surfaces by using density functional theory. The surface model is constructed using hybrid quantum mechanics/molecular mechanics embedding techniques and is based on an oxygen-terminated diamond (110) surface. We use the hybrid quantum mechanics/molecular mechanics method to assess the ability of different density functional approximations to predict the adsorption structure, energy, and barrier for diffusion on pristine and defective surfaces. We find that surface defects (vacancies and surface dopants) strongly anchor adatoms on vacancy sites. We further investigated the thermal stability of gold adatoms, which reveals high barriers associated with lateral diffusion away from the vacancy site. The result provides an explanation for the high stability of experimentally imaged single gold adatoms on BDD and a starting point to investigate the early stages of nucleation during metal surface deposition.

Understanding the Excited-State Relaxation Mechanisms of Xanthophyll Lutein by Multi-configurational Electronic Structure Calculations

The contradictory behaviors in light harvesting and non-photochemical quenching make xanthophyll lutein the most attractive functional molecule in photosynthesis. Despite several theoretical simulations on the spectral properties and excited-state dynamics, the atomic-level photophysical mechanisms need to be further studied and established, especially for an accurate description of geometric and electronic structures of conical intersections for the lowest several electronic states of lutein. In the present work, semiempirical OM2/MRCI and multi-configurational restricted active space self-consistent field methods were performed to optimize the minima and conical intersections in and between the 1Ag-, 2Ag-, 1Bu+, and 1Bu- states. Meanwhile, the relative energies were refined by MS-CASPT2(10,8)/6-31G*, which can reproduce correct electronic state properties as those in the spectroscopic experiments. Based on the above calculation results, we proposed a possible excited-state relaxation mechanism for lutein from its initially populated 1Bu+ state. Once excited to the optically bright 1Bu+ state, the system will propagate along the key reaction coordinate, i.e., the stretching vibration of the conjugated carbon chain. During this period of time, the 1Bu- state will participate in and forms a resonance state between the 1Bu- and 1Bu+ states. Later, the system will rapidly hop to the 2Ag- state via the 1Bu+/2Ag- conical intersection. Finally, the lutein molecule will survive in the 2Ag- state for a relatively long time before it internally converts to the ground state directly or via a twisted S1/S0 conical intersection. Notably, though the photophysical picture may be very different in solvents and proteins, the current theoretical study proposed a promising calculation protocol and also provided many valuable mechanistic insights for lutein and similar carotenoids.

Computational infrared and Raman spectra by hybrid QM/MM techniques: a study on molecular and catalytic material systems

Vibrational spectroscopy is one of the most well-established and important techniques for characterizing chemical systems. To aid the interpretation of experimental infrared and Raman spectra, we report on recent theoretical developments in the ChemShell computational chemistry environment for modelling vibrational signatures. The hybrid quantum mechanical and molecular mechanical approach is employed, using density functional theory for the electronic structure calculations and classical forcefields for the environment. Computational vibrational intensities at chemical active sites are reported using electrostatic and fully polarizable embedding environments to achieve more realistic vibrational signatures for materials and molecular systems, including solvated molecules, proteins, zeolites and metal oxide surfaces, providing useful insight into the effect of the chemical environment on the signatures obtained from experiment. This work has been enabled by the efficient task-farming parallelism implemented in ChemShell for high-performance computing platforms.  This article is part of a discussion meeting issue 'Supercomputing simulations of advanced materials'.

Non-covalent interactions that tune the reactivities of the flavins in bifurcating electron transferring flavoprotein

Bifurcating electron transferring flavoproteins (Bf-ETFs) tune chemically identical flavins to two contrasting roles. To understand how, we used hybrid quantum mechanical molecular mechanical calculations to characterize non-covalent interactions applied to each flavin by the protein. Our computations replicated the differences between the reactivities of the flavins: the electron transferring flavin (^ETflavin) was calculated to stabilize anionic semiquinone (ASQ) as needed to execute its single-electron transfers, whereas the Bf flavin (^Bfflavin) was found to disfavor the ASQ state more than does free flavin and to be less susceptible to reduction. The stability of ^ETflavin ASQ was attributed in part to H-bond donation to the flavin O2 from a nearby His side chain, via comparison of models employing different tautomers of His. This H-bond between O2 and the ET site was uniquely strong in the ASQ state, whereas reduction of ^ETflavin to the anionic hydroquinone (AHQ) was associated with side chain reorientation, backbone displacement and reorganization of its H-bond network including a Tyr from the other domain and subunit of the ETF. The Bf site was less responsive overall, but formation of the ^Bfflavin AHQ allowed a nearby Arg side chain to adopt an alternative rotamer that can H-bond to the ^Bfflavin O4. This would stabilize the anionic ^Bfflavin and rationalize effects of mutation at this position. Thus, our computations provide insights on states and conformations that have not been possible to characterize experimentally, offering explanations for observed residue conservation and raising possibilities that can now be tested.

The application of QM/MM simulations in heterogeneous catalysis

The QM/MM simulation method is provenly efficient for the simulation of biological systems, where an interplay of extensive environment and delicate local interactions drives a process of interest through a funnel on a complex energy landscape. Recent advances in quantum chemistry and force-field methods present opportunities for the adoption of QM/MM to simulate heterogeneous catalytic processes, and their related systems, where similar intricacies exist on the energy landscape. Herein, the fundamental theoretical considerations for performing QM/MM simulations, and the practical considerations for setting up QM/MM simulations of catalytic systems, are introduced; then, areas of heterogeneous catalysis are explored where QM/MM methods have been most fruitfully applied. The discussion includes simulations performed for adsorption processes in solvent at metallic interfaces, reaction mechanisms within zeolitic systems, nanoparticles, and defect chemistry within ionic solids. We conclude with a perspective on the current state of the field and areas where future opportunities for development and application exist.

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.

Toward a Consistent Prediction of Defect Chemistry in CeO2

Polarizable shell-model potentials are widely used for atomic-scale modeling of charged defects in solids using the Mott-Littleton approach and hybrid Quantum Mechanical/Molecular Mechanical (QM/MM) embedded-cluster techniques. However, at the pure MM level of theory, the calculated defect energetics may not satisfy the requirement of quantitative predictions and are limited to only certain charged states. Here, we proposed a novel interatomic potential development scheme that unifies the predictions of all relevant charged defects in CeO2 based on the Mott-Littleton approach and QM/MM electronic-structure calculations. The predicted formation energies of oxygen vacancies accompanied by different excess electron localization patterns at the MM level of theory reach the accuracy of density functional theory (DFT) calculations using hybrid functionals. The new potential also accurately reproduces a wide range of physical properties of CeO2, showing excellent agreement with experimental and other computational studies. These findings provide opportunities for accurate large-scale modeling of the partial reduction and nonstoichiometry in CeO2, as well as a prototype for developing robust interatomic potentials for other defective crystals.

Influence of Solvent on Selective Catalytic Reduction of Nitrogen Oxides with Ammonia over Cu-CHA Zeolite

The copper-exchanged zeolite Cu-CHA has received considerable attention in recent years, owing to its application in the selective catalytic reduction (SCR) of NO_x species. Here, we study the NH₃-SCR reaction mechanism on Cu-CHA using the hybrid quantum mechanical/molecular mechanical (QM/MM) technique and investigate the effects of solvent on the reactivity of active Cu species. To this end, a comparison is made between water- and ammonia-solvated and bare Cu species. The results show the promoting effect of solvent on the oxidation component of the NH₃-SCR cycle since the formation of important nitrate species is found to be energetically more favorable on the solvated Cu sites than in the absence of solvent molecules. Conversely, both solvent molecules are predicted to inhibit the reduction component of the NH₃-SCR cycle. Diffuse reflectance infrared fourier-transform spectroscopy (DRIFTS) experiments exploiting (concentration) modulation excitation spectroscopy (MES) and phase-sensitive detection (PSD) identified spectroscopic signatures of Cu-nitrate and Cu-nitrosamine (H₂NNO), important species which had not been previously observed experimentally. This is further supported by the QM/MM-calculated harmonic vibrational analysis. Additional insights are provided into the reactivity of solvated active sites and the formation of key intermediates including their formation energies and vibrational spectroscopic signatures, allowing the development of a detailed understanding of the reaction mechanism. We demonstrate the role of solvated active sites and their influence on the energetics of important species that must be explicitly considered for an accurate understanding of NH₃-SCR kinetics.

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.

General embedded cluster protocol for accurate modeling of oxygen vacancies in metal-oxides

The O vacancy (Ov) formation energy, EOv, is an important property of a metal-oxide, governing its performance in applications such as fuel cells or heterogeneous catalysis. These defects are routinely studied with density functional theory (DFT). However, it is well-recognized that standard DFT formulations (e.g., the generalized gradient approximation) are insufficient for modeling the Ov, requiring higher levels of theory. The embedded cluster method offers a promising approach to compute EOv accurately, giving access to all electronic structure methods. Central to this approach is the construction of quantum(-mechanically treated) clusters placed within suitable embedding environments. Unfortunately, current approaches to constructing the quantum clusters either require large system sizes, preventing application of high-level methods, or require significant manual input, preventing investigations of multiple systems simultaneously. In this work, we present a systematic and general quantum cluster design protocol that can determine small converged quantum clusters for studying the Ov in metal-oxides with accurate methods, such as local coupled cluster with single, double, and perturbative triple excitations. We apply this protocol to study the Ov in the bulk and surface planes of rutile TiO2 and rock salt MgO, producing the first accurate and well-converged determinations of EOv with this method. These reference values are used to benchmark exchange–correlation functionals in DFT, and we find that all the studied functionals underestimate EOv, with the average error decreasing along the rungs of Jacob’s ladder. This protocol is automatable for high-throughput calculations and can be generalized to study other point defects or adsorbates.

Quantum Chemistry Common Driver and Databases (QCDB) and Quantum Chemistry Engine (QCEngine): Automation and interoperability among computational chemistry programs

Community efforts in the computational molecular sciences (CMS) are evolving toward modular, open, and interoperable interfaces that work with existing community codes to provide more functionality and composability than could be achieved with a single program. The Quantum Chemistry Common Driver and Databases (QCDB) project provides such capability through an application programming interface (API) that facilitates interoperability across multiple quantum chemistry software packages. In tandem with the Molecular Sciences Software Institute and their Quantum Chemistry Archive ecosystem, the unique functionalities of several CMS programs are integrated, including CFOUR, GAMESS, NWChem, OpenMM, Psi4, Qcore, TeraChem, and Turbomole, to provide common computational functions, i.e., energy, gradient, and Hessian computations as well as molecular properties such as atomic charges and vibrational frequency analysis. Both standard users and power users benefit from adopting these APIs as they lower the language barrier of input styles and enable a standard layout of variables and data. These designs allow end-to-end interoperable programming of complex computations and provide best practices options by default.

Flyby reaction trajectories: Chemical dynamics under extrinsic force

Dynamic effects are an important determinant of chemical reactivity and selectivity, but the deliberate manipulation of atomic motions during a chemical transformation is not straightforward. Here, we demonstrate that extrinsic force exerted upon cyclobutanes by stretching pendant polymer chains influences product selectivity through force-imparted nonstatistical dynamic effects on the stepwise ring-opening reaction. The high product stereoselectivity is quantified by carbon-13 labeling and shown to depend on external force, reactant stereochemistry, and intermediate stability. Computational modeling and simulations show that, besides altering energy barriers, the mechanical force activates reactive intramolecular motions nonstatistically, setting up "flyby trajectories" that advance directly to product without isomerization excursions. A mechanistic model incorporating nonstatistical dynamic effects accounts for isomer-dependent mechanochemical stereoselectivity.

Computer-aided understanding and engineering of enzymatic selectivity

Enzymes offering chemo-, regio-, and stereoselectivity enable the asymmetric synthesis of high-value chiral molecules. Unfortunately, the drawback that naturally occurring enzymes are often inefficient or have undesired selectivity toward non-native substrates hinders the broadening of biocatalytic applications. To match the demands of specific selectivity in asymmetric synthesis, biochemists have implemented various computer-aided strategies in understanding and engineering enzymatic selectivity, diversifying the available repository of artificial enzymes. Here, given that the entire asymmetric catalytic cycle, involving precise interactions within the active pocket and substrate transport in the enzyme channel, could affect the enzymatic efficiency and selectivity, we presented a comprehensive overview of the computer-aided workflow for enzymatic selectivity. This review includes a mechanistic understanding of enzymatic selectivity based on quantum mechanical calculations, rational design of enzymatic selectivity guided by enzyme-substrate interactions, and enzymatic selectivity regulation via enzyme channel engineering. Finally, we discussed the computational paradigm for designing enzyme selectivity in silico to facilitate the advancement of asymmetric biosynthesis.

QMCube (QM3 ): An all-purpose suite for multiscale QM/MM calculations

QMCube (QM³ ) is a suite written in the Python programming language, initially focused on multiscale QM/MM simulations of biological systems, but open enough to address other kinds of problems. It allows the user to combine highly efficient QM and MM programs, providing unified access to a wide range of computational methods. The suite also supplies additional modules with extra functionalities. These modules facilitate common tasks such as performing the setup of the models or process the data generated during the simulations. The design of QM³ has been carried out considering the least number of external dependencies (only an algebra library, already included in the distribution), which makes it extremely portable. Also, the modular structure of the suite should help to expand and develop new computational methods.

Two-photon conversion of a bacterial phytochrome

In nature, sensory photoreceptors underlie diverse spatiotemporally precise and generally reversible biological responses to light. Photoreceptors also serve as genetically encoded agents in optogenetics to control by light organismal state and behavior. Phytochromes represent a superfamily of photoreceptors that transition between states absorbing red light (Pr) and far-red light (Pfr), thus expanding the spectral range of optogenetics to the near-infrared range. Although light of these colors exhibits superior penetration of soft tissue, the transmission through bone and skull is poor. To overcome this fundamental challenge, we explore the activation of a bacterial phytochrome by a femtosecond laser emitting in the 1 μm wavelength range. Quantum chemical calculations predict that bacterial phytochromes possess substantial two-photon absorption cross sections. In line with this notion, we demonstrate that the photoreversible Pr ↔ Pfr conversion is driven by two-photon absorption at wavelengths between 1170 and 1450 nm. The Pfr yield was highest for wavelengths between 1170 and 1280 nm and rapidly plummeted beyond 1300 nm. By combining two-photon activation with bacterial phytochromes, we lay the foundation for enhanced spatial resolution in optogenetics and unprecedented penetration through bone, skull, and soft tissue.

f90wrap: an automated tool for constructing deep Python interfaces to modern Fortran codes

f90wrap is a tool to automatically generate Python extension modules which interface to Fortran libraries that makes use of derived types. It builds on the capabilities of the popular f2py utility by generating a simpler Fortran 90 interface to the original Fortran code which is then suitable for wrapping with f2py, together with a higher-level Pythonic wrapper that makes the existance of an additional layer transparent to the final user. f90wrap has been used to wrap a number of large software packages of relevance to the condensed matter physics community, including the QUIP molecular dynamics code and the CASTEP density functional theory code.

Balanced polarizable Drude force field parameters for molecular anions: phosphates, sulfates, sulfamates, and oxides

Polarizable force fields are emerging as a more accurate alternative to additive force fields in terms of modeling and simulations of a variety of chemicals including biomolecules. Explicit treatment of induced polarization in charged species such as phosphates and sulfates offers the potential for achieving an improved atomistic understanding of the physical forces driving their interactions with their environments. To help achieve this, in this study we present balanced Drude polarizable force field parameters for molecular ions including phosphates, sulfates, sulfamates, and oxides. Better balance was primarily achieved in the relative values of minimum interaction energies and distances of the anionic model compounds with water at the Drude and quantum mechanical (QM) model chemistries. Parametrization involved reoptimizing available parameters as well as extending the force field to new molecules with the goal of achieving self-consistency with respect to the Lennard-Jones and electrostatic parameters targeting QM and experimental hydration free energies. The resulting force field parameters achieve consistent treatment across the studied anions, facilitating more balanced simulations of biomolecules and small organic molecules in the context of the classical Drude polarizable force field. Graphical abstract.

NWChem: Past, present, and future

Specialized computational chemistry packages have permanently reshaped the landscape of chemical and materials science by providing tools to support and guide experimental efforts and for the prediction of atomistic and electronic properties. In this regard, electronic structure packages have played a special role by using first-principle-driven methodologies to model complex chemical and materials processes. Over the past few decades, the rapid development of computing technologies and the tremendous increase in computational power have offered a unique chance to study complex transformations using sophisticated and predictive many-body techniques that describe correlated behavior of electrons in molecular and condensed phase systems at different levels of theory. In enabling these simulations, novel parallel algorithms have been able to take advantage of computational resources to address the polynomial scaling of electronic structure methods. In this paper, we briefly review the NWChem computational chemistry suite, including its history, design principles, parallel tools, current capabilities, outreach, and outlook.

fromage: A library for the study of molecular crystal excited states at the aggregate scale

The study of photoexcitations in molecular aggregates faces the twofold problem of the increased computational cost associated with excited states and the complexity of the interactions among the constituent monomers. A mechanistic investigation of these processes requires the analysis of the intermolecular interactions, the effect of the environment, and 3D arrangements or crystal packing on the excited states. A considerable number of techniques have been tailored to navigate these obstacles; however, they are usually restricted to in‐house codes and thus require a disproportionate effort to adopt by researchers approaching the field. Herein, we present the FRamewOrk for Molecular AGgregate Excitations (fromage), which implements a collection of such techniques in a Python library complemented with ready‐to‐use scripts. The program structure is presented and the principal features available to the user are described: geometrical analysis, exciton characterization, and a variety of ONIOM schemes. Each is illustrated by examples of diverse organic molecules in condensed phase settings. The program is available at https://github.com/Crespo-Otero-group/fromage.

Advances in Sustainable Catalysis: A Computational Perspective

The enormous challenge of moving our societies to a more sustainable future offers several exciting opportunities for computational chemists. The first principles approach to "catalysis by design" will enable new and much greener chemical routes to produce vital fuels and fine chemicals. This prospective outlines a wide variety of case studies to underscore how the use of theoretical techniques, from QM/MM to unrestricted DFT and periodic boundary conditions, can be applied to biocatalysis and to both homogeneous and heterogenous catalysts of all sizes and morphologies to provide invaluable insights into the reaction mechanisms they catalyze.