Stepping-Stone CBH: Benchmark and Application of a Multilayered Isodesmic-Based Correction Scheme.

We present a generalization of the connectivity-based hierarchy (CBH) of isodesmic-based correction schemes to a multilayered fragmentation platform for overall cost reduction while retaining high accuracy. The newly developed multilayered CBH approach, called stepping-stone CBH (SSCBH), is benchmarked on a diverse set of 959 medium-sized organic molecules. Applying SSCBH corrections to the PBEh-D3 density functional resulted in an average error of 0.76 kcal/mol for the full test set compared to accurate CCSD(T)-quality enthalpies and an even lower error of 0.44 kcal/mol on a subset containing only acyclic molecules. These results rival the traditional CBH-3 approach at a greatly reduced cost, allowing larger fragment corrections to be made at the MP2 level of theory rather than with G4. Our SSCBH approach will enable more widespread applications of CBH methods to a broader range of organic and biomolecular systems.

Toward Post-Hartree-Fock Accuracy for Protein-Ligand Affinities Using the Molecules-in-Molecules Fragmentation-Based Method.

The complexity and size of large molecular systems, such as protein-ligand complexes, pose computational challenges for accurate post-Hartree-Fock calculations. This study delivers a thorough benchmarking of the Molecules-in-Molecules (MIM) method, presenting a clear and accessible strategy for layer/theory selections in post-Hartree-Fock computations on substantial molecular systems, notably protein-ligand complexes. An approach is articulated, enabling augmented computational efficiency by strategically canceling out common subsystem energy terms between complexes and proteins within the supermolecular equation. Employing DLPNO-based post-Hartree-Fock methods in conjunction with the three-layer MIM method (MIM3), this study demonstrates the achievement of protein-ligand binding energies with remarkable accuracy (errors <1 kcal mol-1), while significantly reducing computational costs. Furthermore, noteworthy correlations between theoretically computed interaction energies and their experimental equivalents were observed, with R2 values of approximately 0.90 and 0.78 for CDK2 and BZT-ITK sets, respectively, thus validating the efficacy of the MIM method in calculating binding energies. By highlighting the crucial role of diffuse or small Pople-style basis sets in the middle layer for reducing energy errors, this work provides valuable insights and practical methodologies for interaction energy computations in large molecular complexes and opens avenues for their application across a diverse range of molecular systems.

Leveraging DFT and Molecular Fragmentation for Chemically Accurate pKa Prediction Using Machine Learning.

We present a quantum mechanical/machine learning (ML) framework based on random forest to accurately predict the pKas of complex organic molecules using inexpensive density functional theory (DFT) calculations. By including physics-based features from low-level DFT calculations and structural features from our connectivity-based hierarchy (CBH) fragmentation protocol, we can correct the systematic error associated with DFT. The generalizability and performance of our model are evaluated on two benchmark sets (SAMPL6 and Novartis). We believe the carefully curated input of physics-based features lessens the model's data dependence and need for complex deep learning architectures, without compromising the accuracy of the test sets. As a point of novelty, our work extends the applicability of CBH, employing it for the generation of viable molecular descriptors for ML.

Solvent-Dependent Folding Behavior of a Helix-Forming Aryl-Triazole Foldamer.

Aromatic foldamers make up a novel class of bioinspired molecules that display helical conformations and have functions that rely on control over their coil-helix folding preferences. While the folding has been extensively examined by experiment, it has rarely been paired with the types of atomic level insights offered by theory. We present the results of all-atom molecular dynamics (MD) simulations to examine the role of solvent polarity on driving the helical folding behavior of the aryl-triazole foldamer. The temperature-dependent enhanced sampling technique, replica-exchange MD simulations, was employed to understand the folding phenomena. The simulation results show that in a low polarity solvent (dichloromethane), the foldamer prefers to stay in the unfolded state. The unfolded state has four dipolar triazoles (5 D) in their favored all-anti geometries and favoring anti-parallel geometries. However, an increase in solvent polarity using acetonitrile resulted in solvophobic collapse, yielding the helically folded form as the predominant state. The folded helix has an all-syn geometry, with triazoles in parallel arrangements. Intermediate conformations with a mixture of syn and anti arrangements of the triazoles are of lower abundance in both the DCM and MeCN solvents. The chiral handedness of the helix observed experimentally is assigned as left-handed by correlation with computed electronic circular dichroism spectra using time-dependent density functional theory.

Vibrational Analysis of Constrained Molecular Systems.

Vibrational spectroscopy, including infrared (IR), Raman spectroscopy, and vibrational circular dichroism, is instrumental in determining the structure and composition of molecules. These techniques are highly sensitive to molecular conformations. However, full molecular optimization, necessary for theoretical vibrational spectra, can lead to unintended conformational changes, especially in large biomolecules like polypeptides. To address this, dihedral angle constraints can be imposed during optimization to preserve the molecule's native conformation. Constraint-optimized molecular geometries, not being true stationary points in the full configurational space, pose challenges for traditional vibrational analysis. We address this by considering such geometries as subspace minima, reformulating vibrational analysis to incorporate constraints. Normal modes and spectra consistent with these constraints are obtained by projecting the force constant matrix onto a space orthogonal to the constrained coordinates. This method, illustrated by the example of enkephalin, yields 3N - 6 - m nonzero frequencies after constraint projection, demonstrating its applicability to biomolecules with flexible conformations. Our approach offers a comprehensive mathematical framework to compute vibrational spectra of molecules with conformationally flexible subunits under environmental constraints.

Fragment-based models for dissociation of strong acids in water: Electrostatic embedding minimizes the dependence on the fragmentation schemes.

Fragmentation methods such as MIM (Molecules-in-Molecules) provide a route to accurately model large systems and have been successful in predicting their structures, energies, and spectroscopic properties. However, their use is often limited to systems at equilibrium due to the inherent complications in the choice of fragments in systems away from equilibrium. Furthermore, the presence of charges resulting from any heterolytic bond breaking may increase the fragmentation error. We have previously suggested EE-MIM (Electrostatically Embedded Molecules-In-Molecules) as a method to mitigate the errors resulting from the missing long-range interactions in molecular clusters in equilibrium. Here, we show that the same method can be applied to improve the performance of MIM to solve the longstanding problem of dependency of the fragmentation energy error on the choice of the fragmentation scheme. We chose four widely used acid dissociation reactions (HCl, HClO4, HNO3, and H2SO4) as test cases due to their importance in chemical processes and complex reaction potential energy surfaces. Electrostatic embedding improves the performance at both one and two-layer MIM as shown by lower EE-MIM1 and EE-MIM2 errors. The EE-MIM errors are also demonstrated to be less dependent on the choice of the fragmentation scheme by analyzing the variation in fragmentation energy at the points with more than one possible fragmentation scheme (points where the fragmentation scheme changes). EE-MIM2 with M06-2X as the low-level resulted in a variation of less than 1 kcal/mol for all the cases and 1 kJ/mol for all but three cases, rendering our method fragmentation scheme-independent for acid dissociation processes.

Structure Dependence of CO2 Reduction Electrocatalyzed by Metal-Nanographene Complexes: A Computational Study.

Improving the energy efficiency of electrocatalytic reduction of CO2 requires tuning of redox properties of electrocatalysts to match redox potentials of the substrate. Recently, we introduced nanographenes as ligands for metal complexes for such purposes by taking advantage of size-dependent properties of the conjugated systems. Here, we use computations to investigate the structure dependence of the electrocatalysis at Re(diimine)(CO)3Cl complexes with nanographene ligands that contain a polycyclic aromatic hydrocarbon moiety through a pyrazinyl linkage. We show that the reduction potentials of the complexes depend not only on conjugation size but also on shape and geometry of the ligands, revealing another parameter in tuning the redox properties of the electrocatalysts. In addition, our work reveals a compromise between reduction potentials and activation of this class of electrocatalysts.

MIM-ML: A Novel Quantum Chemical Fragment-Based Random Forest Model for Accurate Prediction of NMR Chemical Shifts of Nucleic Acids.

We developed a random forest machine learning (ML) model for the prediction of 1H and 13C NMR chemical shifts of nucleic acids. Our ML model is trained entirely on reproducing computed chemical shifts obtained previously on 10 nucleic acids using a Molecules-in-Molecules (MIM) fragment-based density functional theory (DFT) protocol including microsolvation effects. Our ML model includes structural descriptors as well as electronic descriptors from an inexpensive low-level semiempirical calculation (GFN2-xTB) and trained on a relatively small number of DFT chemical shifts (2080 1H chemical shifts and 1780 13C chemical shifts on the 10 nucleic acids). The ML model is then used to make chemical shift predictions on 8 new nucleic acids ranging in size from 600 to 900 atoms and compared directly to experimental data. Though no experimental data was used in the training, the performance of our model is excellent (mean absolute deviation of 0.34 ppm for 1H chemical shifts and 2.52 ppm for 13C chemical shifts for the test set), despite having some nonstandard structures. A simple analysis suggests that both structural and electronic descriptors are critical for achieving reliable predictions. This is the first attempt to combine ML from fragment-based DFT calculations to predict experimental chemical shifts accurately, making the MIM-ML model a valuable tool for NMR predictions of nucleic acids.

Investigating the Stacking Interactions Responsible for Driving 3D Self-Association of Tricarb Macrocycles.

We have investigated the noncovalent forces that play a crucial role in the three-dimensional (3D) self-association of the tricarb macrocycle (composed of alternating triazoles and carbazoles) to understand the multilayer stacks observed through electron microscopy. To explore this idea quantitatively, we have investigated a stacked dimer model of tricarb, where we consider homochiral as well as heterochiral forms of the dimer. We have computed the rotational potential energy surface of the dimer by conducting an angle-dependent scan between the two macrocycles using different levels of theory including the RI-MP2 ab initio method. We observe that dimers oriented at an angle of 60° exhibit the highest stability, while a secondary minimum is observed at an angle of 30°. While density functional theory (DFT) describes the behavior of both minima very close to that obtained with RI-MP2, semiempirical and MM models appear to obtain only a shoulder instead of the second minimum. To further understand the underlying interactions responsible for stabilizing the self-assembly of the macrocycles, we employed energy decomposition analysis (EDA) using SAPT0. This quantitative assessment allowed us to identify the major contributing noncovalent interactions, including electrostatic, exchange-repulsion, dispersion, and induction. Finally, we expanded our study to evaluate the accuracy of the MIM (molecules-in-molecules) fragmentation methodology in capturing the crucial π-stacking interactions. Our benchmarking results using the MIM method accurately replicated the angle-dependent PES results. This shows the efficacy of MIM in predicting the noncovalent interactions responsible for the construction of 3D and other higher-order nanoarchitectures for tricarb and related compounds.

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.

Triplet States of Cyanostar and Its Anion Complexes.

The design of advanced optical materials based on triplet states requires knowledge of the triplet energies of the molecular building blocks. To this end, we report the triplet energy of cyanostar (CS) macrocycles, which are the key structure-directing units of small-molecule ionic isolation lattices (SMILES) that have emerged as programmable optical materials. Cyanostar is a cyclic pentamer of covalently linked cyanostilbene units that form π-stacked dimers when binding anions as 2:1 complexes. The triplet energies, ET, of the parent cyanostar and its 2:1 complex around PF6- are measured to be 1.96 and 2.02 eV, respectively, using phosphorescence quenching studies at room temperature. The similarity of these triplet energies suggests that anion complexation leaves the triplet energy relatively unchanged. Similar energies (2.0 and 1.98 eV, respectively) were also obtained from phosphorescence spectra of the iodinated form, I-CS, and of complexes formed with PF6- and IO4- recorded at 85 K in an organic glass. Thus, measures of the triplet energies likely reflect geometries close to those of the ground state either directly by triplet energy transfer to the ground state or indirectly by using frozen media to inhibit relaxation. Density functional theory (DFT) and time-dependent DFT were undertaken on a cyanostar analogue, CSH, to examine the triplet state. The triplet excitation localizes on a single olefin whether in the single cyanostar or its π-stacked dimer. Restriction of the geometrical changes by forming either a dimer of macrocycles, (CSH)2, or a complex, (CSH)2·PF6-, reduces the relaxation resulting in an adiabatic energy of the triplet state of 2.0 eV. This structural constraint is also expected for solid-state SMILES materials. The obtained T1 energy of 2.0 eV is a key guide line for the design of SMILES materials for the manipulation of triplet excitons by triplet state engineering in the future.

Tautomerization of H+KPGG: Entropic Consequences of Strong Hydrogen-Bond Networks in Peptides.

Ion mobility spectrometry-mass spectrometry and quantum chemical calculations are used to determine the structures and stabilities of the singly protonated peptide H+KPGG. The two peaks making up the IMS distribution are shown to be tautomers differing by the location of the extra proton on either the lysine side chain or the N-terminus. The lysine-protonated tautomer is strongly preferred entropically while being disfavored in terms of the electronic energy and enthalpy. This relationship is shown, through comparison of all low-lying conformers of both tautomers, to be related to the strong hydrogen-bond network of the N-terminally protonated tautomer. A general relationship is demonstrated wherein stronger cross-peptide hydrogen-bond networks result in entropically disfavored conformers. Further effects of the H+KPGG hydrogen-bond network are probed by computationally examining singly and doubly methylated analogues. These results demonstrate the importance of the entropic consequences of hydrogen bonds to peptide stability as well as techniques for perturbing the hydrogen-bond network and folding preferences of peptides via minimal chemical modification.

Approaching Coupled Cluster Accuracy with Density Functional Theory Using the Generalized Connectivity-Based Hierarchy.

This Perspective reviews connectivity-based hierarchy (CBH), a systematic hierarchy of error-cancellation schemes developed in our group with the goal of achieving chemical accuracy using inexpensive computational techniques ("coupled cluster accuracy with DFT"). The hierarchy is a generalization of Pople's isodesmic bond separation scheme that is based only on the structure and connectivity and is applicable to any organic and biomolecule consisting of covalent bonds. It is formulated as a series of rungs involving increasing levels of error cancellation on progressively larger fragments of the parent molecule. The method and our implementation are discussed briefly. Examples are given for the applications of CBH involving (1) energies of complex organic rearrangement reactions, (2) bond energies of biofuel molecules, (3) redox potentials in solution, (4) pKa predictions in the aqueous medium, and (5) theoretical thermochemistry combining CBH with machine learning. They clearly show that near-chemical accuracy (1-2 kcal/mol) is achieved for a variety of applications with DFT methods irrespective of the underlying density functional used. They demonstrate conclusively that seemingly disparate results, often seen with different density functionals in many chemical applications, are due to an accumulation of systematic errors in the smaller local molecular fragments that can be easily corrected with higher-level calculations on those small units. This enables the method to achieve the accuracy of the high level of theory (e.g., coupled cluster) while the cost remains that of DFT. The advantages and limitations of the method are discussed along with areas of ongoing developments.

Interpretable Graph-Network-Based Machine Learning Models via Molecular Fragmentation.

Chemists have long benefitted from the ability to understand and interpret the predictions of computational models. With the current shift to more complex deep learning models, in many situations that utility is lost. In this work, we expand on our previously work on computational thermochemistry and propose an interpretable graph network, FragGraph(nodes), that provides decomposed predictions into fragment-wise contributions. We demonstrate the usefulness of our model in predicting a correction to density functional theory (DFT)-calculated atomization energies using Δ-learning. Our model predicts G4(MP2)-quality thermochemistry with an accuracy of <1 kJ mol-1 for the GDB9 dataset. Besides the high accuracy of our predictions, we observe trends in the fragment corrections which quantitatively describe the deficiencies of B3LYP. Node-wise predictions significantly outperform our previous model predictions from a global state vector. This effect is most pronounced as we explore the generality by predicting on more diverse test sets indicating node-wise predictions are less sensitive to extending machine learning models to larger molecules.

Cooperativity in Photofoldamer Chloride Double Helices Turned On with Sequences and Solvents, Around with Guests, and Off with Light.

Photofoldamers are sequence-defined receptors capable of switching guest binding on and off. When two foldamer strands wrap around the guest into 2:1 double helical complexes, cooperativity emerges, and with it comes the possibility to switch cooperativity with light and other stimuli. We use lessons from nonswitchable sequence isomers of aryl-triazole foldamers to guide how to vary the sequence location of azobenzenes from the end (FEND) to the interior (FIN) and report their impact on the cooperative formation of 2:1 complexes with Cl-. This sequence change produces a 125-fold increase from anti-cooperative (α = 0.008) for FEND to non-cooperative with FIN (α = 1.0). Density functional theory (DFT) studies show greater H-bonding and a more relaxed double helix for FIN. The solvent and guest complement the synthetic designs. Use of acetonitrile to enhance solvophobicity further enhances cooperativity in FIN (α = 126) but lowers the difference in cooperativity between sequences. Surprisingly, the impact of the sequence on cooperativity is inverted when the guest size is increased from Cl- (3.4 Å) to BF4- (4.1 Å). While photoconversion of interior azobenzenes was poor, the cis-cis isomer forms 1:1 complexes around chloride consistent with switching cooperativity. The effect of the guest, solvent, and light on the double-helix cooperativity depends on the sequence.

Quantitative Prediction of Vertical Ionization Potentials from DFT via a Graph-Network-Based Delta Machine Learning Model Incorporating Electronic Descriptors.

While accurate wave function theories like CCSD(T) are capable of modeling molecular chemical processes, the associated steep computational scaling renders them intractable for treating large systems or extensive databases. In contrast, density functional theory (DFT) is much more computationally feasible yet often fails to quantitatively describe electronic changes in chemical processes. Herein, we report an efficient delta machine learning (ΔML) model that builds on the Connectivity-Based Hierarchy (CBH) scheme─an error correction approach based on systematic molecular fragmentation protocols─and achieves coupled cluster accuracy on vertical ionization potentials by correcting for deficiencies in DFT. The present study integrates concepts from molecular fragmentation, systematic error cancellation, and machine learning. First, we show that by using an electron population difference map, ionization sites within a molecule may be readily identified, and CBH correction schemes for ionization processes may be automated. As a central feature of our work, we employ a graph-based QM/ML model, which embeds atom-centered features describing CBH fragments into a computational graph to further increase accuracy for the prediction of vertical ionization potentials. In addition, we show that the incorporation of electronic descriptors from DFT, namely electron population difference features, improves model performance well beyond chemical accuracy (1 kcal/mol) to approach benchmark accuracy. While the raw DFT results are strongly dependent on the underlying functional used, for our best models, the performance is robust and much less dependent on the functional.

Accurate and Cost-Effective NMR Chemical Shift Predictions for Nucleic Acids Using a Molecules-in-Molecules Fragmentation-Based Method.

We have developed, implemented, and assessed an efficient protocol for the prediction of NMR chemical shifts of large nucleic acids using our molecules-in-molecules (MIM) fragment-based quantum chemical approach. To assess the performance of our approach, MIM-NMR calculations are calibrated on a test set of three nucleic acids, where the structure is derived from solution-phase NMR studies. For DNA systems with multiple conformers, the one-layer MIM method with trimer fragments (MIM1trimer) is benchmarked to get the lowest energy structure, with an average error of only 0.80 kcal/mol with respect to unfragmented full molecule calculations. The MIMI-NMRdimer calibration with respect to unfragmented full molecule calculations shows a mean absolute deviation (MAD) of 0.06 and 0.11 ppm, respectively, for 1H and 13C nuclei, but the performance with respect to experimental NMR chemical shifts is comparable to the more expensive MIM1-NMR and MIM2-NMR methods with trimer subsystems. To compare with the experimental chemical shifts, a standard protocol is derived using DNA systems with Protein Data Bank (PDB) IDs 1SY8, 1K2K, and 1KR8. The effect of structural minimizations is employed using a hybrid mechanics/semiempirical approach and used for computations in solution with implicit and explicit-implicit solvation models in our MIM1-NMRdimer methodology. To demonstrate the applicability of our protocol, we tested it on seven nucleic acids, including structures with nonstandard residues, heteroatom substitutions (F and B atoms), and side chain mutations with a size ranging from ∼300 to 1100 atoms. The major improvement for predicted MIM1-NMRdimer calculations is obtained from structural minimizations and implicit solvation effects. A significant improvement with the explicit-implicit solvation model is observed only for two smaller nucleic acid systems (1KR8 and 7NBK), where the expensive first solvation shell is replaced by the microsolvation model, in which a single water molecule is added for each solvent-exposed amino and imino protons, along with the implicit solvation. Overall, our target accuracy of ∼0.2-0.3 ppm for 1H and ∼2-3 ppm for 13C has been achieved for large nucleic acids. The proposed MIM-NMR approach is accurate and cost-effective (linear scaling with system size), and it can aid in the structural assignments of a wide range of complex biomolecules.

Proton-Coupled, Low-Energy Pathway for Electrocatalytic CO2 Reduction at Re(Diimine) Complexes with a Conjugated Pyrazinyl Moiety.

Large conjugated carbon framework has been incorporated as the diimine ligand for Re(α-diimine)(CO)3Cl complexes with a pyrazinyl linkage, either to increase energy efficiency or to turn them into heterogeneous catalysts for selective electrocatalytic CO2 reduction. However, there exists a nonmonotonic dependence of CO2 reduction overpotential on the conjugation size of the ligands. Understanding its origin could facilitate heterogenization of molecular catalysts with improved energy efficiency. Here, we show that the conjugated pyrazinyl moiety plays a crucial role in catalysis by enabling a proton-coupled, lower-energy pathway for CO2 reduction. With ligands of moderate size, the pathway leads to previously unknown intermediates and decreases CO2 reduction overpotential. Because the pathway hinges on the basicity of the pyrazinyl nitrogen, we propose that it imposes a limit on the conjugation size of the ligand for the pathway to be effective.

ONIOM Method with Charge Transfer Corrections (ONIOM-CT): Analytic Gradients and Benchmarking.

Hybrid methods such as ONIOM that treat different regions of a large molecule using different methods are widely used to investigate chemical reactions in a variety of materials and biological systems. However, there are inherent sources of significant errors due to the standard treatment of the boundary between the regions using hydrogen link atoms. In particular, an unbalanced charge distribution in the chemically important model region is a potential source of such problems. We have previously suggested ONIOM-CT (ONIOM with charge transfer corrections) which addresses this issue by applying a potential in the form of point charges to obtain a desired charge redistribution. The metric for charge redistribution relies on the type of population analysis used to obtain the charges. ONIOM-CT has been implemented using Mulliken and Löwdin population analyses and has been shown to improve computed reaction energies for illustrative chemical reactions. In this work, we derive and implement the analytic gradients for ONIOM-CT that requires solving two sets of coupled-perturbed self-consistent equations, one each for the model system and the full system. However, both are needed only at the low level of theory, allowing for an efficient formulation and implementation for both Mulliken and Löwdin population analyses. Benchmarking and illustrative geometry optimizations have been carried out for a previously studied set of reactions involving a single link atom between regions. Additionally, we have generalized our method for the treatment of model systems involving multiple link atoms to enable applications for a broader set of problems. The generalized methods are illustrated for both charge models. Furthermore, we have studied a set of three proton transfer reactions and demonstrate that significant improvement is achieved by ONIOM-CT over ONIOM using both Mulliken and Löwdin population analyses.