Quantum Mechanics/molecular Mechanics Research Articles

Quantum mechanics/molecular mechanics (QM/MM) methods in drug design: a comprehensive review of development and applications

The integration of quantum mechanics (QM) and molecular mechanics (MM) methods, known as QM/MM, has emerged as a powerful computational approach in drug design. This hybrid technique combines the accuracy of QM calculations for the reactive region with the computational efficiency of MM methods for the surrounding environment. QM/MM methods have proven invaluable in studying chemical reactions, exploring enzyme mechanisms, and investigating ligand-protein interactions, all of which are crucial for rational drug design. This review provides a comprehensive overview of the QM/MM methodology, its theoretical foundations, and its applications in various aspects of drug design. We discussed the key components, such as QM and MM region partitioning, link atom schemes, and boundary treatments. Additionally, we highlight recent advancements and challenges in QM/MM methods, including polarizable force fields, implicit solvation models, and enhanced sampling techniques. Furthermore, we present illustrative examples showcasing the successful application of QM/MM methods in lead optimization, virtual screening, and the elucidation of biochemical mechanisms relevant to drug targets. Finally, we provide perspectives on future developments and the potential impact of QM/MM methods on the drug discovery pipeline.

The Interaction Mechanism Between C14-Polyacetylene Compounds and the Rat TRPA1 Receptor: An In Silico Study.

Polyacetylene (PA) compounds, as natural products, exhibit remarkable properties and distinctive chemical activities. Three structurally similar C14-PA compounds-Echinophorin D, Echinophorin B, and Echinophorin A-extracted from plants demonstrate varying biological activities on the Transient Receptor Potential Channel A1 (TRPA1) protein, which belongs to the TRP (Transient Receptor Potential) family. In the current study, we investigated the binding modes of these three PA compounds with TRPA1 using molecular dynamics (MD), molecular docking, binding free energy calculations, and quantum mechanics/molecular mechanics (QM/MM) methods. Initially, a putative binding site (site-II) in TRPA1 was identified for these compounds; Echinophorin B was found to stabilize the upward A-loop of TRPA1, which is critical for its activation. Furthermore, the binding affinity calculations of PA compounds through molecular fragment decomposition indicate that the arrangement of two triple bonds and one double bond in C14-PA compounds is vital for regulating TRPA1 bioactivity. Additionally, the lipophilic and electronic properties of the three molecules were analyzed in relation to binding affinity, establishing a correlation between TRPA1 activity and these molecular properties.

Assessing the Partial Hessian Approximation in QM/MM-Based Vibrational Analysis.

The partial Hessian approximation is often used in vibrational analysis of quantum mechanics/molecular mechanics (QM/MM) systems because calculating the full Hessian matrix is computationally impractical. This approach aligns with the core concept of QM/MM, which focuses on the QM subsystem. Thus, using the partial Hessian approximation implies that the main interest is in the local vibrational modes of the QM subsystem. Here, we investigate the accuracy and applicability of the partial Hessian vibrational analysis (PHVA) approach as it is typically used within QM/MM, i.e., only the Hessian belonging to the QM subsystem is computed. We focus on solute-solvent systems with small, rigid solutes. To separate two of the major sources of errors, we perform two separate analyses. First, we study the effects of the partial Hessian approximation on local normal modes, harmonic frequencies, and harmonic IR and Raman intensities by comparing them to those obtained using full Hessians, where both partial and full Hessians are calculated at the QM level. Then, we quantify the errors introduced by QM/MM used with the PHVA by comparing normal modes, frequencies, and intensities obtained using partial Hessians calculated using a QM/MM-type embedding approach to those obtained using partial Hessians calculated at the QM level. Another aspect of the PHVA is the appearance of normal modes resembling the translation and rotation of the QM subsystem. These pseudotranslational and pseudorotational modes should be removed as they are collective vibrations of the atoms in the QM subsystem relative to a frozen MM subsystem and, thus, not well-described. We show that projecting out translation and rotation, usually done for systems in isolation, can adversely affect other normal modes. Instead, the pseudotranslational and pseudorotational modes can be identified and removed.

Dissociation Energies via Embedding Techniques.

Due to the large number of interactions, evaluating interaction energies for large or periodic systems results in time-consuming calculations. Prime examples are liquids, adsorbates, and molecular crystals. Thus, there is a high demand for a cheap but still accurate method to determine interaction energies and gradients. One approach to counteract the computational cost is to fragment a large cluster into smaller subsystems, sometimes called many-body expansion, with the fragments being molecules or parts thereof. These subsystems can then be embedded into larger entities, representing the bigger system. In this work, we test several subsystem approaches and explore their limits and behaviors, determined by calculations of trimer interaction energies. The methods presented here encompass mechanical embedding, point charges, polarizable embedding, polarizable density embedding, and density embedding. We evaluate nonembedded fragmentation, QM/MM (quantum mechanics/molecular mechanics), and QM/QM (quantum mechanics/quantum mechanics) embedding theories. Finally, we make use of symmetry-adapted perturbation theory utilizing density functional theory for the monomers to interpret the results. Depending on the strength of the interaction, different embedding methods and schemes prove favorable to accurately describe a system. The embedding approaches presented here are able to decrease the interaction energy errors with respect to full system calculations by a factor of up to 20 in comparison to simple/unembedded approaches, leading to errors below 0.1 kJ/mol.

Constrained nuclear-electronic orbital QM/MM approach for simulating complex systems with quantum nuclear delocalization effects incorporated

The hybrid quantum mechanics/molecular mechanics (QM/MM) approach, which combines the accuracy of QM methods with the efficiency of MM methods, is widely used in the study of complex systems. However, past QM/MM implementations often neglect or face challenges in addressing nuclear quantum effects, despite their crucial role in many key chemical and biological processes. Recently, our group developed the constrained nuclear-electronic orbital (CNEO) theory, a cost-efficient approach that accurately addresses nuclear quantum effects, especially quantum nuclear delocalization effects. In this work, we integrate CNEO with the QM/MM approach through the electrostatic embedding scheme and apply the resulting CNEO QM/MM to two hydrogen-bonded complexes. We find that both solvation effects and nuclear quantum effects significantly impact hydrogen bond structures and dynamics. Notably, in the glutamic acid–glutamate complex, which mimics a common low barrier hydrogen bond in biological systems, CNEO QM/MM accurately predicts nearly equal proton sharing between the two residues. With an accurate description of both quantum nuclear delocalization effects and environmental effects, CNEO QM/MM is a promising new approach for simulating complex chemical and biological systems.

A Dual-Focus Workflow for Simultaneously Engineering High Activity and Thermal Stability in Methyl Parathion Hydrolase.

Industrial fermentation applications typically require enzymes that exhibit high stability and activity at high temperatures. However, efforts to simultaneously improve these properties are usually limited by a trade-off between stability and activity. This report describes a computational strategy to enhance both activity and thermal stability of the mesophilic organophosphate-degrading enzyme, methyl parathion hydrolase (MPH). To predict hotspot mutation sites, we assembled a library of features associated with the target properties for each residue and then prioritized candidate sites by hierarchical clustering. Subsequent in silico screening with multiple algorithms to simulate selective pressures yielded a subset of 23 candidate mutations. Iterative parallel screening of mutations that improved thermal stability and activity yielded, MPHase-m5b, which exhibited 13.3 °C higher Tm and 4.2 times higher catalytic activity than wild-type (WT) MPH over a wide temperature range. Systematic analysis of crystal structures, molecular dynamics (MD) simulations, and quantum mechanics/molecular mechanics (QM/MM) calculations revealed a wider entrance to the active site that increased substrate access with an extensive network of interactions outside the active site that reinforced αβ/βα sandwich architecture to improve thermal stability. This study thus provides an advanced, rational design framework to improve efficiency in engineering highly active, thermostable biocatalysts for industrial applications.

A Dual‐Focus Workflow for Simultaneously Engineering High Activity and Thermal Stability in Methyl Parathion Hydrolase

AbstractIndustrial fermentation applications typically require enzymes that exhibit high stability and activity at high temperatures. However, efforts to simultaneously improve these properties are usually limited by a trade‐off between stability and activity. This report describes a computational strategy to enhance both activity and thermal stability of the mesophilic organophosphate‐degrading enzyme, methyl parathion hydrolase (MPH). To predict hotspot mutation sites, we assembled a library of features associated with the target properties for each residue and then prioritized candidate sites by hierarchical clustering. Subsequent in silico screening with multiple algorithms to simulate selective pressures yielded a subset of 23 candidate mutations. Iterative parallel screening of mutations that improved thermal stability and activity yielded, MPHase‐m5b, which exhibited 13.3 °C higher Tm and 4.2 times higher catalytic activity than wild‐type (WT) MPH over a wide temperature range. Systematic analysis of crystal structures, molecular dynamics (MD) simulations, and quantum mechanics/molecular mechanics (QM/MM) calculations revealed a wider entrance to the active site that increased substrate access with an extensive network of interactions outside the active site that reinforced αβ/βα sandwich architecture to improve thermal stability. This study thus provides an advanced, rational design framework to improve efficiency in engineering highly active, thermostable biocatalysts for industrial applications.

End Group Effects on Anion Binding in Tetraglycine Peptide: A Computational Study.

The importance of anions in various processes has led to a search for molecules that can effectively recognize and interact with these anions. This study explores how the tetraglycine (Gly)4 peptide in its zwitterionic, neutral, and terminally capped forms acts as a receptor for H2PO4- and HSO4- anions within the framework of supramolecular host-guest chemistry. Using molecular dynamics (MD) simulations, we obtained the conformations of the receptor-anion complexes. Density functional theory (DFT), quantifies the complexes' interaction energies in both gas and solvent phases. Proton transfer within the zwitterionic complex with H2PO4- anion alters peptide charge distribution, affecting its conformation and binding site arrangement, as analysed by quantum mechanics/molecular mechanics (QM/MM) methods. Symmetry-adapted perturbation theory (SAPT) and noncovalent interactions analysis highlight the role of electrostatic interactions in these receptor-anion complexes. It emphasizes the key interactions such as N-H···O and O-H···O=C between the peptide backbone and anions and elucidates the molecular recognition mechanism driven by crucial noncovalent interactions. The termination of the peptide's end groups modulates anion binding sites from the backbone to the charged N-terminal, resulting in distinct binding sites. Our findings provide insights for designing peptides tailored to function as anion receptors in diverse supramolecular chemistry applications.

An insight into the solitonic features of the nonlinear generalized higher-order Schrödinger equation using the solver method

For many nonlinear applications described by the dynamics of nonlinear Schrödinger equation with higher-order terms (HONLSE) such as nonlinear optics, space plasma physics molecular biology, astrophysics, quantum mechanics, superfluid, fluid mechanics, and fiber optics communications, a unique closed-form solution have been obtained using energy equation. In addition, some new solitary solutions HONLSE have been obtained via the unified solver method. The resultant solutions behave as breathers, super solitons, envelope breathers, blow up, localized super waves, periodical super shock, train super solitons, and shock structures. The modulations of Kerr nonlinear, chromatic dispersive, and wave packet drift parameters on the wave characteristics of the obtained solutions have been investigated. It was reported that the model parameters affect the amplitude, steepness, and width of the resultant structures. The provided solution can be used as a box solver for a variety of natural science systems described by distinct nonlinear equations.

Influence of Wobbling Tryptophan and Mutations on PET Degradation Explored by QM/MM Free Energy Calculations.

Plastic-degrading enzymes, particularly poly(ethylene terephthalate) (PET) hydrolases, have garnered significant attention in recent years as potential eco-friendly solutions for recycling plastic waste. However, understanding of their PET-degrading activity and influencing factors remains incomplete, impeding the development of uniform approaches for enhancing PET hydrolases for industrial applications. A key aspect of PET hydrolase engineering is optimizing the PET-hydrolysis reaction by lowering the associated free energy barrier. However, inconsistent findings have complicated these efforts. Therefore, our goal is to elucidate various aspects of enzymatic PET degradation by means of quantum mechanics/molecular mechanics (QM/MM) reaction simulations and analysis, focusing on the initial reaction step, acylation, in two thermophilic PET hydrolases, LCC and PES-H1, along with their highly active variants, LCCIG and PES-H1FY. Our findings highlight the impact of semiempirical QM methods on proton transfer energies, affecting the distinction between a two-step reaction involving a metastable tetrahedral intermediate and a one-step reaction. Moreover, we uncovered a concerted conformational change involving the orientation of the PET benzene ring, altering its interaction with the side-chain of the "wobbling" tryptophan from T-stacking to parallel π-π interactions, a phenomenon overlooked in prior research. Our study thus enhances the understanding of the acylation mechanism of PET hydrolases, in particular by characterizing it for the first time for the promising PES-H1FY using QM/MM simulations. It also provides insights into selecting a suitable QM method and a reaction coordinate, valuable for future studies on PET degradation processes.

Mechanism of Nucleic Acid Phosphodiester Bond Cleavage by Human Endonuclease V: MD and QM/MM Calculations Reveal a Versatile Metal Dependence.

Human endonuclease V (EndoV) catalytically removes deaminated nucleobases by cleaving the phosphodiester bond as part of RNA metabolism. Despite being implicated in several diseases (cancers, cardiovascular diseases, and neurological disorders) and potentially being a useful tool in biotechnology, details of the human EndoV catalytic pathway remain unclear due to limited experimental information beyond a crystal structure of the apoenzyme and select mutational data. Since a mechanistic understanding is critical for further deciphering the central roles and expanding applications of human EndoV in medicine and biotechnology, molecular dynamics (MD) simulations and quantum mechanics/molecular mechanics (QM/MM) calculations were used to unveil the atomistic details of the catalytic pathway. Due to controversies surrounding the number of metals required for nuclease activity, enzyme-substrate models with different numbers of active site metals and various metal-substrate binding configurations were built based on structural data for other nucleases. Subsequent MD simulations revealed the structure and stability of the human EndoV-substrate complex for a range of active site metal binding architectures. Four unique pathways were then characterized using QM/MM that vary in metal number (one versus two) and modes of substrate coordination [direct versus indirect (water-mediated)], with several mechanisms being fully consistent with experimental structural, kinetic, and mutational data for related nucleases, including members of the EndoV family. Beyond uncovering key roles for several active site amino acids (D240 and K155), our calculations highlight that while one metal is essential for human EndoV activity, the enzyme can benefit from using two metals due to the presence of two suitable metal binding sites. By directly comparing one- versus two-metal-mediated P-O bond cleavage reactions within the confines of the same active site, our work brings a fresh perspective to the "number of metals" controversy.

Constrained Catalytic Itinerary of a Retaining 3,6‐Anhydro‐D‐Galactosidase, a Key Enzyme in Red Algal Cell Wall Degradation

AbstractThe marine Bacteroidota Zobellia galactanivorans has a polysaccharide utilization locus dedicated to the catabolism of the red algal cell wall galactan carrageenan and its unique and industrially important α‐3,6‐anhydro‐D‐galactose (ADG) monosaccharide. Here we present the first analysis of the specific molecular interactions that the exo‐(α‐1,3)‐3,6‐anhydro‐D‐galactosidase ZgGH129 uses to cope with the strict steric restrictions imposed by its bicyclic ADG substrate — which is ring flipped relative to D‐galactose. Crystallographic snapshots of key catalytic states obtained with the natural substrate and novel chemical tools designed to mimic species along the reaction coordinate, together with quantum mechanics/molecular mechanics (QM/MM) metadynamics methods and kinetic studies, demonstrate a retaining mechanism where the second step is rate limiting. The conformational landscape of the constrained 3,6‐anhydro‐D‐galactopyranose ring proceeds through enzyme glycosylation B1,4→[E4]≠→E4/1C4 and deglycosylation E4/1C4→[E4]≠→B1,4 itineraries limited to the Southern Hemisphere of the Cremer–Pople sphere. These results demonstrate the conformational changes throughout catalysis in a non‐standard, sterically restrained, bicyclic monosaccharide, and provide a molecular framework for mechanism‐based inhibitor design for anhydro‐type carbohydrate‐processing enzymes and for future applications involving carrageenan degradation. In addition, our study provides a rare example of distinct niche‐based conformational itineraries within the same carbohydrate‐active enzyme family.

Constrained Catalytic Itinerary of a Retaining 3,6-Anhydro-D-Galactosidase, a Key Enzyme in Red Algal Cell Wall Degradation.

The marine Bacteroidota Zobellia galactanivorans has a polysaccharide utilization locus dedicated to the catabolism of the red algal cell wall galactan carrageenan and its unique and industrially important α-3,6-anhydro-D-galactose (ADG) monosaccharide. Here we present the first analysis of the specific molecular interactions that the exo-(α-1,3)-3,6-anhydro-D-galactosidase ZgGH129 uses to cope with the strict steric restrictions imposed by its bicyclic ADG substrate - which is ring flipped relative to D-galactose. Crystallographic snapshots of key catalytic states obtained with the natural substrate and novel chemical tools designed to mimic species along the reaction coordinate, together with quantum mechanics/molecular mechanics (QM/MM) metadynamics methods and kinetic studies, demonstrate a retaining mechanism where the second step is rate limiting. The conformational landscape of the constrained 3,6-anhydro-D-galactopyranose ring proceeds through enzyme glycosylation B1,4→[E4]≠→E4/1C4 and deglycosylation E4/1C4→[E4]≠→B1,4 itineraries limited to the Southern Hemisphere of the Cremer-Pople sphere. These results demonstrate the conformational changes throughout catalysis in a non-standard, sterically restrained, bicyclic monosaccharide, and provide a molecular framework for mechanism-based inhibitor design for anhydro-type carbohydrate-processing enzymes and for future applications involving carrageenan degradation. In addition, our study provides a rare example of distinct niche-based conformational itineraries within the same carbohydrate-active enzyme family.

Influence of Selective Deoxyfluorination on the Molecular Structure of Type-2 N-Acetyllactosamine.

N-Acetyllactosamine is a common saccharide motif found in various biologically active glycans. This motif usually works as a backbone for additional modifications and thus significantly influences glycan conformational behavior and biological activity. In this work, we have investigated the type-2 N-acetyllactosamine scaffold using the complete series of its monodeoxyfluorinated analogs. These glycomimetics have been studied by molecular mechanics, quantum mechanics, X-ray crystallography, and various NMR techniques, which have provided a comprehensive and complete insight into the role of individual hydroxyl groups in the conformational behavior and lipophilicity of N-acetyllactosamine.

Molecular modelling of luciferyl adenylate deprotonation in the active site of Photinus pyralis luciferase

Firefly bioluminescence is a product of chemical reactions that involve luciferin chromophore oxidation in the active site of luciferase proteins. We perform a series of classical molecular dynamic simulations and combined quantum and molecular mechanics (QM/MM) calculations to expose the molecular mechanism of C4 carbon atom deprotonation in luciferyl adenylate molecule. QM/MM calculations confirm that ND-protonated His245 residue is a suitable proton acceptor in the wild-type Photinus pyralis luciferase. Classical molecular dynamic simulations reveal oxygen-binding cavities inside the protein including the one located close to the C4 atom of luciferin. In mutant forms that lack direct interactions with the His245 side chain, a proton wire comprising water molecules stimulates either protonation of the phosphate group of the luciferyl adenylate forming an unstable intermediate or keto–enol tautomerisation of luciferin. The comparison of the ionisation potentials of molecular systems with the ‘deprotonated’ C4 carbon atom revealed that the ionisation energy for the enol tautomeric form is close to the system with the His245 proton acceptor. Thus, the existence of the keto–enol tautomerisation channel might explain bioluminescence in case of the absence of the amino acid proton acceptor.

Electron Attachment to the Nucleobase Uracil in Diethylene Glycol: The Signature of a Doorway.

The cellular environment plays a significant role in low energy electron-mediated radiation damage to genetic materials. In this study, we have modeled the effect of the bulk medium on electron attachment to nucleobases in diethylene glycol (DEG) using uracil as a test case, in accordance with recent experimental work on the observation of dissociative quasi-free electron attachment to nucleoside via excited anion radical in solution (in DEG). Our EOM-CCSD-based quantum mechanical/molecular mechanical (QM/MM) simulations indicate that the electron scavenging by uracil in DEG is much slower than that observed in the aqueous medium due to its viscosity. This work also establishes that a doorway mechanism exists in uracil microsolvated and bulk solvated with DEG, with the dipole-bound state and solvent-bound state acting as doorway states, respectively.

Pomotrelvir and Nirmatrelvir Binding and Reactivity with SARS-CoV-2 Main Protease: Implications for Resistance Mechanisms from Computations.

We investigate the inhibition mechanism between pomotrelvir and the SARS-CoV-2 main protease using molecular mechanics and quantum mechanics/molecular mechanics simulations. Alchemical transformations where each Pi group of pomotrelvir was transformed into its counterpart in nirmatrelvir were performed to unravel the individual contribution of each group to the binding and reaction processes. We have shown that while a γ-lactam ring is preferred at position P1, a δ-lactam ring is a good alternative for the design of inhibitors for variants presenting mutations at position 166. For the P2 position, tertiary amines are preferred with respect to secondary amines. Flexible side chains at the P2 position can disrupt the preorganization of the active site, favouring the exploration of non-reactive conformations. The substitution of the P2 group of pomotrelvir by that of nirmatrelvir resulted in a compound, named as C2, that presents a better binding free energy and a higher population of reactive conformations in the Michaelis complex. Analysis of the chemical reaction to form the covalent complex has shown a similar reaction mechanism and activation free energies for pomotrelvir, nirmatrelvir and C2. We hope that these findings could be useful to design better inhibitors to fight present and future variants of the SARS-CoV-2 virus.

Rosmarinic Acid as a Potential Multi-targeted Inhibitor for SAR-CoV-2: An In silico Virtual Screening Approach

Background: Rosmarinic acid, a natural compound found in various plants like rosemary and lemon balm, may have potential as a multi-targeted inhibitor for SARS-CoV-2, a strain of virus responsible for COVID-19. SARS-CoV-2, a fusion protein of S1 and S2 subunits, has multiple precursors angiotensin-converting enzyme2 (ACE2), transmembrane serine protease 2 (TMPRSS2), papain-like protease (PLpro), and 3-chymotrypsin-like protease (3CLpro). The chemical interaction of Rosmarinic acid with SARS-CoV-2 is of major interest reported here. Objective:: The quantitative study of Rosmarinic acid with various precursors of SARS-CoV-2 has been accounted for in detail. Furthermore, the conformational flexibility of Rosmarinic acid has also been investigated during the interaction with four different precursors of SARS-CoV-2. Methods: This investigation delves deeply into the analysis of various aspects, including geometric parameters, atomic charge, the energy gap between the highest occupied and lowest unoccupied molecular orbitals, dipole moments, and the analysis of non-covalent interactions (NCI). Furthermore, the study incorporates molecular docking techniques in conjunction with thorough quantum chemical calculations to provide comprehensive insights. Results: Rosmarinic acid shows promise as a versatile inhibitor of SARS-CoV-2, the virus responsible for COVID-19. It can target multiple key precursors of the virus, including TMPRSS2, angiotensin- converting enzyme2, 3CLpro, and PLpro, found in the fusion protein comprising S1 and S2 subunits. This study delves into the quantitative analysis of Rosmarinic acid's interactions with these precursors. Its adaptable structure allows it to engage with them effectively. Various molecular parameters, including atomic charge, energy gap between molecular orbitals, dipole moment, and noncovalent interactions, are comprehensively explored. Conclusion: Combining molecular docking and quantum mechanics, the findings suggest Rosmarinic acid's potential as a multi-targeted SARS-CoV-2 inhibitor.

Experimental and numerical study of the decomposition, product spectrum, and sooting properties of adamantane fuels

This work combined experimental measurements with two theoretical approaches, reactive Molecular Dynamics (MD) simulations and Quantum Mechanics (QM) calculations, to investigate the combustion and sooting properties of three adamantane fuels, adamantane (AD), 1,3-dimethyl-adamantane (DMAD), and 1-ethyl-adamantane (EAD). These fuels were selected since the adamantane fuel family can potentially be used as sustainable aviation fuels and comparing these fuels would reveal the effects of side chain on their combustion and sooting properties. We determined the bond dissociation energies of the three test fuels using QM calculations and found that the functionalized side chains have the weakest bonds and their presence only slightly affects the bond strength in the AD multi-cyclic core. We performed pyrolysis simulations for all three fuels using ReaxFF-based MD simulations and found that DMAD and EAD have higher decomposition rates than AD and also generate more high-molecular-weight products. For all three fuels, these large products were observed to contribute significantly to hydrocarbon growth processes, which lead to large soot nucleating species and even soot-like structures. A higher yield of such soot nucleating compounds was found during the pyrolysis of DMAD and EAD than AD since the decomposition products of DMAD and EAD are more branched and those of AD have mostly straight chains. These theoretical analyses were supported by experimental measurements of yield-based sooting tendencies, which suggest that DMAD and EAD are sootier than AD and that all three fuels are sootier than standard alkane jet fuel surrogates but less sooty than jet fuel aromatic content.

Theoretical Investigation on the Initial Reaction Mechanism of Hexaethynylbenzene on Au(111) Surface.

Graphyne has attracted considerable interest and attention since its successful synthesis, due to its enormous potential for applications in the fields of electronics, energy, catalysis, information technology, etc. Although various methods for synthesizing graphyne have been explored, single-layer graphynes have not been successfully developed. Hexaethynylbenzene (HEB) is considered an ideal precursor molecule because it can undergo Glaser coupling reactions between molecules to synthesize single layer graphdiyne on single crystal metal surfaces via on-surface reactions. Unfortunately, this method fails to achieve the expected results, and the underlying mechanism is not clear. In this work, we employed a combination of ab initio molecular dynamics (AIMD) and quantum mechanics (QM) methods to investigate the initial reaction mechanism of HEB molecules on a Au(111) surface. We revealed that HEB molecules undergo both intermolecular coupling and intramolecular cyclization on the Au(111) surface. The favorable pathways of these two types of reactions were then distinguished, confirming that the distance between the terminal carbon atoms of the ethynyl groups plays an important role in C-C coupling. The insights revealed from this work could facilitate the rational design of precursor molecules and deepen the understanding of the reaction processes.