QM:QM embedding using electronic densities within an ONIOM framework: Energies and analytic gradients

A Vision for the Future of Multiscale Modeling

The rise of modern computer science enabled physical chemistry to make enormous progresses in understanding and harnessing natural and artificial phenomena. Nevertheless, despite the advances achieved over past decades, computational resources are still insufficient to thoroughly simulate extended systems from first principles. Indeed, countless biological, catalytic and photophysical processes require ab initio treatments to be properly described, but the breadth of length and time scales involved makes it practically unfeasible. A way to address these issues is to couple theories and algorithms working at different scales by dividing the system into domains treated at different levels of approximation, ranging from quantum mechanics to classical molecular dynamics, even including continuum electrodynamics. This approach is known as multiscale modeling and its use over the past 60 years has led to remarkable results. Considering the rapid advances in theory, algorithm design, and computing power, we believe multiscale modeling will massively grow into a dominant research methodology in the forthcoming years. Hereby we describe the main approaches developed within its realm, highlighting their achievements and current drawbacks, eventually proposing a plausible direction for future developments considering also the emergence of new computational techniques such as machine learning and quantum computing. We then discuss how advanced multiscale modeling methods could be exploited to address critical scientific challenges, focusing on the simulation of complex light-harvesting processes, such as natural photosynthesis. While doing so, we suggest a cutting-edge computational paradigm consisting in performing simultaneous multiscale calculations on a system allowing the various domains, treated with appropriate accuracy, to move and extend while they properly interact with each other. Although this vision is very ambitious, we believe the quick development of computer science will lead to both massive improvements and widespread use of these techniques, resulting in enormous progresses in physical chemistry and, eventually, in our society.

Development of analytic gradients for the Huzinaga quantum embedding method and its applications to large-scale hybrid and double hybrid DFT forces

The theory of analytic gradients is presented for the projector-based density functional theory (DFT) embedding approach utilizing the Huzinaga-equation. The advantages of the Huzinaga-equation-based formulation are demonstrated. In particular, it is shown that the projector employed does not appear in the Lagrangian, and the potential risk of numerical problems is avoided at the evaluation of the gradients. The efficient implementation of the analytic gradient theory is presented for approaches where hybrid DFT, second-order Møller-Plesset perturbation theory, or double hybrid DFT are embedded in lower-level DFT environments. To demonstrate the applicability of the method and to gain insight into its accuracy, it is applied to equilibrium geometry optimizations, transition state searches, and potential energy surface scans. Our results show that bond lengths and angles converge rapidly with the size of the embedded system. While providing structural parameters close to high-level quality for the embedded atoms, the embedding approach has the potential to relax the coordinates of the environment as well. Our demonstrations on a 171-atom zeolite and a 570-atom protein system show that the Huzinaga-equation-based embedding can accelerate (double) hybrid gradient computations by an order of magnitude with sufficient active regions and enables affordable force evaluations or geometry optimizations for molecules of hundreds of atoms.

Hydrogen Sulfide Ligation in Hemoglobin I of Lucina pectinata─A QM/MM and Local Mode Study

In this study, we investigated the interaction between the H2S ligand and the heme pocket of hemoglobin I (HbI) of Lucina pectinata for the wild-type protein; three known mutations where distal glutamine is replaced by hydrophobic valine (Gln64Val) and hydrophilic histidine in both protonation forms (Gln64Hisϵ and Gln64Hisδ); five known mutations of the so-called phenyl cage, replacing the hydrophobic phenylalanines Phe29 and Phe43 with tyrosine (Tyr), valine (Val), or leucine (Leu); and two additional mutations, Phe68Tyr and Phe68Val, in order to complement previous studies with new insights about the binding mechanism at the molecular level. A particular focus was on the intrinsic strengths of the chemical bonds involved, utilizing local vibrational force constants based on combined quantum mechanical-molecular mechanical calculations. Wild-type protein and mutations clustered into two distinct groups: Group 1 protein systems with a proton acceptor in the distal protein pocket, close to one of the H2S bonds, and Group 2 protein systems without a hydrogen acceptor close by in the active site of the protein. According to our results, the interactions between H2S and HbI of Lucina pectinata involve two important elements, namely, binding of H2S to Fe of the heme group, followed by the proton transfer from the HS bond to the distal residue. The distal residue is additionally stabilized by a second proton transfer from the distal residue to COO- of the propionate group in heme. We could identify the FeS bond as a key player and discovered that the strength of this bond depends on two mutual factors, namely, the strength of the HS bond involved in the proton transfer and the electrostatic field of the protein pocket qualifying the FeS bond as a sensitive probe for monitoring changes in H2S ligation upon protein mutations. We hope our study will inspire and guide future experimental studies, targeting new promising mutations such as Phe68Tyr, Phe68Val, or Phe43Tyr/Phe68Val.

EE-ONIOM-CT Method to Efficiently Account for the Missing Interactions in ONIOM: Energies and Analytic Gradients

Hybrid methods such as ONIOM (QM:QM) are widely used for the study of local processes in large systems. However, the intrinsic need for system partitioning often leads to a less-than-desirable performance for some important chemical processes. This is due to the missing interactions in the chemically important model region (i.e., active site) of the high-level theory. The missing interactions can be categorized into two classes, viz. charge transfer (i.e., charge redistribution) and long-range electrostatic interactions. Our group presented two entirely different methods to treat these deficiencies individually. ONIOM-CT and ONIOM-EE methods have been demonstrated to improve the performance of ONIOM by incorporating charge transfer and missing electrostatic interactions, respectively. In general, the inclusion of the missing interactions separately in two different calculations may not be sufficient to reach a high accuracy. Thus, it is highly desirable to develop a method to correct both deficiencies simultaneously. In this work, we aim to connect the methods ONIOM-CT and ONIOM-EE for a more comprehensive treatment. A "stepwise" model was found to be necessary for a robust performance. This model employs a stepwise procedure by first satisfying the ONIOM-CT condition for charge balance before accounting for the electrostatic interactions from the rest of the system perturbatively. This has the advantage of easy interpretation due to the clear separation of the two effects. We demonstrate the performance of our method using embedding charges determined from a Mulliken population analysis. An efficient analytic gradient expression for this method is derived and implemented by requiring three sets of z-vector self-consistent equations. The performance of our method is assessed against full system calculations in high-level theory for a set of three proton transfer reactions representing different degrees of electrostatic embedding.

A density fitting scheme for the fast evaluation of molecular electrostatic potential

Molecular electrostatic potential (MEP) is a significant and crucial physical quantity that can be applied to a large number of scenarios, such as the prediction of nucleophilic or electrophilic attacks, fitting atomic charges, σ-hole, and so forth. The computational cost for the MEP has an O(N² ) scaling with the increase of atoms, which is intractable and laborious for macromolecules. Herein, a density fitting molecular electrostatic potential (DF-MEP) is used to reduce the computational costs for the macromolecular MEP. It is found that the accuracy of DF-MEP is almost identical to the conventional molecular electrostatic potential (Conv-MEP), while the computational costs can be reduced to an O(N) scaling, for example, the computational time of 699,200 grids for the Trp-cage molecule (304 atoms) only takes 16.6 s at the B3LYP-D3(BJ)/def2-SVP level of theory with 16 CPU cores compared with 3060.2 s for the Conv-MEP method.

Computational prediction of the supramolecular self-assembling properties of organic molecules: the role of conformational flexibility of amide moieties

Two families of organic molecules with different backbones have been considered. The first family is based on a macrolactam-like unit that is constrained in a particular conformation. The second family is composed by a substituted central phenyl that allows a larger mobility for its substituents. They have however a common feature, three amide moieties (within the cycle for the macrolactam-like molecule and as substituents for the phenyl) that permit hydrogen bonding when molecules are stacked. In this study we propose a computational protocol to unravel the ability of the different families to self-assemble into organic nanotubes. Starting from the monomer and going towards larger assemblies like dimers, trimers, and pentamers we applied the different protocols to rationalize the behavior of the different assemblies. Both structures and thermodynamics were investigated to give a complete picture of the process. Thanks to the combination of a quantum mechanics approach and molecular dynamics simulations along with the use of tailored tools (non covalent interaction visualization) and techniques (umbrella sampling), we have been able to differentiate the two families and highlight the best candidate for self-assembling purposes.

Aggregation-Induced Emission: A Challenge for Computational Chemistry Taking TPA-BMO as an Example*

A multi-environment computational approach is proposed to study the modulation of the emission behavior of the triphenylamine (Z)-4-benzylidene-2-methyloxazol-5(4H)-one (TPA-BMO) molecule [Tang et al., J. Phys. Chem. C 119, 21875 (2015)]. We aim at (1) proposing a realistic description of the molecule in several environments (solution, aggregate, polymer matrix), (2) modelling its absorption and emission properties, and (3) providing a qualitative understanding of the experimental observations by highlighting the photophysical phenomena leading to the emission modulation. To this purpose, we rely on (TD-)DFT calculations and classical Molecular Dynamics simulations, but also on the hybrid ONIOM QM/QM' approach and the in situ chemical polymerization methodology. In low-polar solvents, the investigation of the potential energy surfaces and the modulation of the emission quantum yield can be attributed to possible photophysical energy dissipation caused by low-frequency vibrational modes. In the aggregate and in the polymer matrix, the emission modulation can be qualitatively interpreted in terms of the possible restriction of the intramolecular vibrations. For these two systems, our study highlights that a careful modelling of the environment is far from trivial but is fundamental to model the optical properties of the fluorophore.

A self-consistent coulomb bath model using density fitting

A self-consistent Coulomb bath model is presented to provide an accurate and efficient way of performing calculations for interfragment electrostatic and polarization interactions. In this method, a condensed-phase system is partitioned into molecular fragment blocks. Each fragment is embedded in the Coulomb bath due to other fragments. Importantly, the present Coulomb bath is represented using a density fitting method in which the electron densities of molecular fragments are fitted using an atom-centered auxiliary basis set of Gaussian type. The Coulomb bath is incorporated into an effective Hamiltonian for each fragment, with which the electron density is optimized through an iterative double self-consistent field (DSCF) procedure to realize the mutual many-body polarization effects. In this work, the accuracy of interfragment interaction energies enumerated using the Coulomb bath is tested, showing a good agreement with the exact results from an energy decomposition analysis. The qualitative features of many-body polarization effects are visualized by electron density difference plots. It is also shown that the present DSCF method can yield fast and robust convergence with near-linear scaling in performance with increase in system size.

Point charge embedding for ONIOM excited states calculations

Hybrid quantum mechanical methods can assist in the interpretation and prediction of the electronic spectra of large molecular structures. In this work, we study the performance of the ONIOM (Our own N-layered Integrated molecular Orbital molecular Mechanics) hybrid method for the calculation of transition energies and oscillator strengths by embedding the core region in a field of fixed point charges. These charges introduce polarization effects from the substituent groups to the core region. We test various charge definitions, with particular attention to the issue of overpolarization near the boundary between layers. To minimize this issue, we fit the charges on the electrostatic potential of the entire structure in the presence of the link atoms used to cap dangling bonds. We propose two constrained fitting strategies: one that produces an average set of charges common to both model system calculations, EE(L1), and one that produces two separate sets of embedding charges, EE(L2). The results from our tests show that indeed electronic embedding with constrained-fitted charges tends to improve the performance of ONIOM compared to non-embedded calculations. However, the EE(L2) charges work best for transition energies, and the EE(L1) charges work best for oscillator strengths. This may be an indication that fixed point charges do not have enough flexibility to adapt to each system, and other effects (e.g., polarization of the embedding field) may be necessary.

Multi-state extrapolation of UV/Vis absorption spectra with QM/QM hybrid methods

In this work, we present a simple approach to simulate absorption spectra from hybrid QM/QM calculations. The goal is to obtain reliable spectra for compounds that are too large to be treated efficiently at a high level of theory. The present approach is based on the extrapolation of the entire absorption spectrum obtained by individual subcalculations. Our program locates the main spectral features in each subcalculation, e.g., band peaks and shoulders, and fits them to Gaussian functions. Each Gaussian is then extrapolated with a formula similar to that of ONIOM (Our own N-layered Integrated molecular Orbital molecular Mechanics). However, information about individual excitations is not necessary so that difficult state-matching across subcalculations is avoided. This multi-state extrapolation thus requires relatively low implementation effort while affording maximum flexibility in the choice of methods to be combined in the hybrid approach. The test calculations show the efficacy and robustness of this methodology in reproducing the spectrum computed for the entire molecule at a high level of theory.

Quantum Mechanics/Molecular Mechanics Modeling of Enzymatic Processes: Caveats and Breakthroughs

Nature has developed large groups of enzymatic catalysts with the aim to transfer substrates into useful products, which enables biosystems to perform all their natural functions. As such, all biochemical processes in our body (we drink, we eat, we breath, we sleep, etc.) are governed by enzymes. One of the problems associated with research on biocatalysts is that they react so fast that details of their reaction mechanisms cannot be obtained with experimental work. In recent years, major advances in computational hardware and software have been made and now large (bio)chemical systems can be studied using accurate computational techniques. One such technique is the quantum mechanics/molecular mechanics (QM/MM) technique, which has gained major momentum in recent years. Unfortunately, it is not a black-box method that is easily applied, but requires careful set-up procedures. In this work we give an overview on the technical difficulties and caveats of QM/MM and discuss work-protocols developed in our groups for running successful QM/MM calculations.

MoD-QM/MM Structural Refinement Method: Characterization of Hydrogen Bonding in the Oxytricha nova G-Quadruplex

A generalization of the Moving-Domain Quantum Mechanics/Molecular Mechanics (MoD-QM/MM) hybrid method [Gascon, J. A.; Leung, S. S. F.; Batista, E. R.; Batista, V. S. J. Chem. Theory Comput. 2006, 2, 175-186] is introduced to provide a self-consistent computational protocol for structural refinement of extended systems. The method partitions the system into molecular domains that are iteratively optimized as quantum mechanical (QM) layers embedded in their surrounding molecular environment to obtain an ab initio quality description of the geometry and the molecular electrostatic potential of the extended system composed of those constituent fragments. The resulting methodology is benchmarked as applied to model systems that allow for full QM optimization as well as through refinement of the hydrogen bonding geometry in Oxytricha nova guanine quadruplex for which several studies have been reported, including the X-ray structure and NMR data. Calculations of (1)H NMR chemical shifts based on the gauge independent atomic orbital (GIAO) method and direct comparisons with experiments show that solvated MoD-QM/MM structures, sampled from explicit solvent molecular dynamics simulations, allow for NMR simulations in much improved agreement with experimental data than models based on the X-ray structure or those optimized using classical molecular mechanics force fields.

Multilevel X-Pol: a fragment-based method with mixed quantum mechanical representations of different fragments

The explicit polarization (X-Pol) method is a fragment-based quantum mechanical model, in which a macromolecular system or other large or complex system in solution is partitioned into monomeric fragments. The present study extends the original X-Pol method, where all fragments are treated using the same electronic structure theory, to multilevel representations, called multilevel X-Pol, in which different electronic structure methods are used to describe different fragments. The multilevel X-Pol method has been implemented into a locally modified version of Gaussian 09. A key ingredient that is used to couple interfragment electrostatic interactions at different levels of theory is the use of the response density for the post-self-consistent-field energy. (The response density is also called the generalized density.) The method is useful for treating fragments in a small region of the system such as a solute molecule or the substrate and amino acids in the active site of an enzyme with a high-level theory, and the fragments in the rest of the system by a lower-level and computationally more efficient method. Multilevel X-Pol is illustrated here by applications to hydrogen bonding complexes in which one fragment is treated with the hybrid M06 density functional, Møller-Plesset perturbation theory, or coupled cluster theory, and the other fragments are treated by Hartree-Fock theory or the B3LYP or M06 hybrid density functionals.