The SCC-DFTB method and its application to biological systems

dxtb-An efficient and fully differentiable framework for extended tight-binding

Automatic differentiation (AD) emerged as an integral part of machine learning, accelerating model development by enabling gradient-based optimization without explicit analytical derivatives. Recently, the benefits of AD and computing arbitrary-order derivatives with respect to any variable were also recognized in the field of quantum chemistry. In this work, we present dxtb-an open-source, fully differentiable framework for semiempirical extended tight-binding (xTB) methods. Developed entirely in Python and leveraging PyTorch for array operations, dxtb facilitates extensibility and rapid prototyping while maintaining computational efficiency. Through comprehensive code vectorization and optimization, we essentially reach the speed of compiled xTB programs for high-throughput calculations of small molecules. The excellent performance also scales to large systems, and batch operability yields additional benefits for execution on parallel hardware. In particular, energy evaluations are on par with existing programs, whereas the speed of automatically differentiated nuclear derivatives is only 2 to 5 times slower compared to their analytical counterparts. We showcase the utility of AD in dxtb by calculating various molecular and spectroscopic properties, highlighting its capacity to enhance and simplify such evaluations. Furthermore, the framework streamlines optimization tasks and offers seamless integration of semiempirical quantum chemistry in machine learning, paving the way for physics-inspired end-to-end differentiable models. Ultimately, dxtb aims to further advance the capabilities of semiempirical methods, providing an extensible foundation for future developments and hybrid machine learning applications. The framework is accessible at https://github.com/grimme-lab/dxtb.

QM/MM Free Energy Calculations of Long-Range Biological Protonation Dynamics by Adaptive and Focused Sampling

Water-mediated proton transfer reactions are central for catalytic processes in a wide range of biochemical systems, ranging from biological energy conversion to chemical transformations in the metabolism. Yet, the accurate computational treatment of such complex biochemical reactions is highly challenging and requires the application of multiscale methods, in particular hybrid quantum/classical (QM/MM) approaches combined with free energy simulations. Here, we combine the unique exploration power of new advanced sampling methods with density functional theory (DFT)-based QM/MM free energy methods for multiscale simulations of long-range protonation dynamics in biological systems. In this regard, we show that combining multiple walkers/well-tempered metadynamics with an extended system adaptive biasing force method (MWE) provides a powerful approach for exploration of water-mediated proton transfer reactions in complex biochemical systems. We compare and combine the MWE method also with QM/MM umbrella sampling and explore the sampling of the free energy landscape with both geometric (linear combination of proton transfer distances) and physical (center of excess charge) reaction coordinates and show how these affect the convergence of the potential of mean force (PMF) and the activation free energy. We find that the QM/MM-MWE method can efficiently explore both direct and water-mediated proton transfer pathways together with forward and reverse hole transfer mechanisms in the highly complex proton channel of respiratory Complex I, while the QM/MM-US approach shows a systematic convergence of selected long-range proton transfer pathways. In this regard, we show that the PMF along multiple proton transfer pathways is recovered by combining the strengths of both approaches in a QM/MM-MWE/focused US (FUS) scheme and reveals new mechanistic insight into the proton transfer principles of Complex I. Our findings provide a promising basis for the quantitative multiscale simulations of long-range proton transfer reactions in biological systems.

Thermal isomerization rates in retinal analogues using Ab-Initio molecular dynamics

For a detailed understanding of chemical processes in nature and industry, we need accurate models of chemical reactions in complex environments. While Eyring transition state theory is commonly used for modeling chemical reactions, it is most accurate for small molecules in the gas phase. A wide range of alternative rate theories exist that can better capture reactions involving complex molecules and environmental effects. However, they require that the chemical reaction is sampled by molecular dynamics simulations. This is a formidable challenge since the accessible simulation timescales are many orders of magnitude smaller than typical timescales of chemical reactions. To overcome these limitations, rare event methods involving enhanced molecular dynamics sampling are employed. In this work, thermal isomerization of retinal is studied using tight-binding density functional theory. Results from transition state theory are compared to those obtained from enhanced sampling. Rates obtained from dynamical reweighting using infrequent metadynamics simulations were in close agreement with those from transition state theory. Meanwhile, rates obtained from application of Kramers' rate equation to a sampled free energy profile along a torsional dihedral reaction coordinate were found to be up to three orders of magnitude higher. This discrepancy raises concerns about applying rate methods to one-dimensional reaction coordinates in chemical reactions.

Generating High-Precision Force Fields for Molecular Dynamics Simulations to Study Chemical Reaction Mechanisms Using Molecular Configuration Transformer

Theoretical studies on chemical reaction mechanisms have been crucial in organic chemistry. Traditionally, calculating the manually constructed molecular conformations of transition states for chemical reactions using quantum chemical calculations is the most commonly used method. However, this way is heavily dependent on individual experience and chemical intuition. In our previous study, we proposed a research paradigm that used enhanced sampling in molecular dynamics simulations to study chemical reactions. This approach can directly simulate the entire process of a chemical reaction. However, the computational speed limited the use of high-precision potential energy functions for simulations. To address this issue, we presented a scheme for training high-precision force fields for molecular modeling using a previously developed graph-neural-network-based molecular model, molecular configuration transformer. This potential energy function allowed for highly accurate simulations at a low computational cost, leading to more precise calculations of the mechanism of chemical reactions. We applied this approach to study a Claisen rearrangement reaction and a carbonyl insertion reaction catalyzed by manganese.

Decarboxylation Mechanism of iso-Orotate Decarboxylase Revisited

iso-Orotate decarboxylase (IDCase), which is involved in the thymidine salvage pathway, has attracted considerable interest owing to its chemical similarity to a hypothetical DNA decarboxylase in mammals. Although valuable insights into the active DNA demethylation of 5-methyl-cytosine can be obtained from the decarboxylation mechanism of 5-carboxyl-uracil (5caU) catalyzed by IDCase, this mechanism remains under debate. In this study, the catalytic mechanism of 5caU decarboxylation by IDCase was studied using hybrid quantum mechanics/molecular mechanics (QM/MM) methodologies and density functional theory (DFT) calculations with a truncated model. The calculations supported a mechanism involving three sequential stages: activation of the 5caU substrate via proton transfer from an arginine (R262') to the carboxyl group of 5caU, formation of a tetrahedral intermediate, and decarboxylation of the tetrahedral intermediate to generate uracil as the product. The reaction pathways and structures obtained using the QM/MM and DFT methods coincided with each other. These simulations provided detailed insights into the unique mechanism of IDCase, clarifying various unresolved issues, such as the critical role of R262'. In addition, aspartate D323 was found to act as a general base in the tetrahedral intermediate formation step and a general acid in the later C-C bond cleavage step.

Exploring Molecular Energy Landscapes by Coupling the DFTB Potential with a Tree-Based Stochastic Algorithm: Investigation of the Conformational Diversity of Phthalates

Exploring the global energy landscape of relatively large molecules at the quantum level is a challenging problem. In this work, we report the coupling of a nonredundant conformational space exploration method, namely, the robotics-inspired iterative global exploration and local optimization (IGLOO) algorithm, with the quantum-chemical density functional tight binding (DFTB) potential. The application of this fast and efficient computational approach to three close-sized molecules of the phthalate family (DBP, BBP, and DEHP) showed that they present different conformational landscapes. These differences have been rationalized by making use of descriptors based on distances and dihedral angles. Coulomb interactions, steric hindrance, and dispersive interactions have been found to drive the geometric properties. A strong correlation has been evidenced between the two dihedral angles describing the side-chain orientation of the phthalate molecules. Our approach identifies low-energy minima without prior knowledge of the potential energy surface, paving the way for future investigations into transition paths and states.

Assessment of Hydrazone Derivatives for Enhanced Steel Corrosion Resistance in 15 wt.% HCl Environments: A Dual Experimental and Theoretical Perspective

This study evaluates the corrosion inhibition capabilities of two novel hydrazone derivatives, (E)-2-(5-methoxy-2-methyl-1H-indol-3-yl)-N'-(4-methylbenzylidene)acetohydrazide (MeHDZ) and (E)-N'-benzylidene-2-(5-methoxy-2-methyl-1H-indol-3-yl)acetohydrazide (HHDZ), on carbon steel in a 15 wt.% HCl solution. A comprehensive suite of analytical techniques, including gravimetric analysis, potentiodynamic polarization (PDP), electrochemical impedance spectroscopy (EIS), and scanning electron microscopy (SEM), demonstrates their significant inhibition efficiency. At an optimal concentration of 5 × 10-3 mol/L, MeHDZ and HHDZ achieve remarkable inhibition efficiencies of 98% and 94%, respectively. EIS measurements reveal a dramatic reduction in effective double-layer capacitance (from 236.2 to 52.8 and 75.3 µF/cm2), strongly suggesting inhibitor adsorption on the steel surface. This effect is further corroborated by an increase in polarization resistance and a significant decrease in corrosion current density at optimal concentrations. Moreover, these inhibitors demonstrate sustained corrosion mitigation over extended exposure durations and maintain effectiveness even under elevated temperatures, highlighting their potential for diverse operational conditions. The adsorption process of these inhibitors aligns well with the Langmuir adsorption isotherm, implying physicochemical interactions at the carbon steel surface. Density functional tight-binding (DFTB) calculations and molecular dynamics simulations provide insights into the inhibitor-surface interaction mechanism, further elucidating the potential of these hydrazone derivatives as highly effective corrosion inhibitors in acidic environments.

Quantum Mechanical/Molecular Mechanical Elucidation of the Catalytic Mechanism of Leukotriene A4 Hydrolase as an Epoxidase

Leukotriene A4 hydrolase (LTA4H) functions as a mono-zinc bifunctional enzyme with aminopeptidase and epoxidase activities. While the aminopeptidase mechanism is well understood, the epoxidase mechanism remains less clear. In continuation of our prior research, we undertook an in-depth exploration of the LTA4H catalytic role as an epoxidase, employing a combined SCC-DFTB/CHARMM method. In the current work, we found that the conversion of LTA4 to leukotriene B4 (LTB4) involves three successive steps: epoxy ring opening (RO), nucleophilic attack (NA), and proton transfer (PT) reactions at the epoxy oxygen atom. Among these steps, the RO and NA stages constitute the potential rate-limiting step within the entire epoxidase mechanism. Notably, the NA step implicates D375 as the general base catalyst, while the PT step engages protonated E271 as the general acid catalyst. Additionally, we delved into the mechanism behind the formation of the isomer product, Δ6-trans-Δ8-cis-LTB4. Our findings debunked the feasibility of a direct LTB4 to iso-LTB4 conversion. Instead, we highlight the possibility of isomerization from LTA4 to its isomeric conjugate (iso-LTA4), showing comparable energy barriers of 5.1 and 5.5 kcal/mol in aqueous and enzymatic environments, respectively. The ensuing dynamics of iso-LTA4 hydrolysis subsequently yield iso-LTB4 via a mechanism akin to LTA4 hydrolysis, albeit with a heightened barrier. Our computations firmly support the notion that substrate isomerization exclusively takes place prior to or during the initial substrate-binding phase, while LTA4 remains the dominant conformer. Notably, our simulations suggest that irrespective of the active site's constrained L-shape, isomerization from LTA4 to its isomeric conjugate remains plausible. The mechanistic insights garnered from our simulations furnish a valuable understanding of LTA4H's role as an epoxidase, thereby facilitating potential advancements in inhibitor design.

Quantum mechanical assessment on the optical properties of capsanthin conformers

As optical properties, the ultraviolet-visible (UV-Vis) absorption spectra of capsanthin-based red natural dye are a decisive parameter for their usage in various applications. Thus, accurately predicting the maximum UV-Vis wavelength (λ max ) values is critical in designing dye-conjugated material. Extensive metadynamics simulations were carried out to generate capsanthin conformers at various levels of the extended tight-binding method. Benchmarking the time-dependent density-functional theory (TD-DFT) methods help understand the results of a particular functional and allows a comparison between results obtained with different functional. The long-range correction (LC) scheme in LC-TD-DFT-D4/ωB97X/def2-SVP has been found to reproduce the experimentalλ max , and exhibited the effect of conformational changes to the calculated wavelengths. On the other hand, an inexpensive yet efficient LC-TD-DFTB method reproduced the experimentalλ max insensitive to conformational changes.

Accurate SCC-DFTB Parametrization of Liquid Water with Improved Atomic Charges and Iterative Boltzmann Inversion

This work presents improvements of the description of liquid water within the self-consistent-charge density-functional based tight-binding scheme combining the use of Weighted Mulliken (WMull) charges and optimized O-H repulsive potential through the iterative Boltzmann inversion (IBI) process. The quality of the newly developed models is validated considering pair radial distribution functions (RDFs), as well as other structural, energetic, thermodynamic, and dynamic properties. The use of WMull charges certainly improves the agreement with experimental data, however leading to over-structured RDFs at short distance, that can be further improved by considering an optimized O-H repulsive potential obtained by the IBI process. Three different schemes were used to optimize this potential: (i) optimization including short O-H distances. This led to accurate RDFs as well as improved self-diffusion coefficient and heat of vaporization, while the proton transfer energy barrier is severely deteriorated; (ii) optimization starting at long distance. The proton transfer energy barrier is recovered while the heat of vaporization is deteriorated and the O-H RDF is less accurate at short distance; (iii) optimization within the path-integral molecular dynamics scheme which allows us to exclude nuclear quantum effects from the repulsive potential. The latter potential, in conjunction with the WMull improved atomic charges, provides similar results as (i) for structural, dynamic, and thermodynamic properties while recovering a large part of the proton transfer energy barrier. It therefore offers a good compromise to study both dynamic properties and chemistry within liquid water at a quantum chemical level.

Two excited-state datasets for quantum chemical UV-vis spectra of organic molecules

We present two open-source datasets that provide time-dependent density-functional tight-binding (TD-DFTB) electronic excitation spectra of organic molecules. These datasets represent predictions of UV-vis absorption spectra performed on optimized geometries of the molecules in their electronic ground state. The GDB-9-Ex dataset contains a subset of 96,766 organic molecules from the original open-source GDB-9 dataset. The ORNL_AISD-Ex dataset consists of 10,502,904 organic molecules that contain between 5 and 71 non-hydrogen atoms. The data reveals the close correlation between the magnitude of the gaps between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO), and the excitation energy of the lowest singlet excited state energies quantitatively. The chemical variability of the large number of molecules was examined with a topological fingerprint estimation based on extended-connectivity fingerprints (ECFPs) followed by uniform manifold approximation and projection (UMAP) for dimension reduction. Both datasets were generated using the DFTB+ software on the "Andes" cluster of the Oak Ridge Leadership Computing Facility (OLCF).

Exploring Routes to Enhance the Calculation of Free Energy Differences via Non-Equilibrium Work SQM/MM Switching Simulations Using Hybrid Charge Intermediates between MM and SQM Levels of Theory or Non-Linear Switching Schemes

Non-equilibrium work switching simulations and Jarzynski's equation are a reliable method for computing free energy differences, ΔAlow→high, between two levels of theory, such as a pure molecular mechanical (MM) and a quantum mechanical/molecular mechanical (QM/MM) description of a system of interest. Despite the inherent parallelism, the computational cost of this approach can quickly become very high. This is particularly true for systems where the core region, the part of the system to be described at different levels of theory, is embedded in an environment such as explicit solvent water. We find that even for relatively simple solute-water systems, switching lengths of at least 5 ps are necessary to compute ΔAlow→high reliably. In this study, we investigate two approaches towards an affordable protocol, with an emphasis on keeping the switching length well below 5 ps. Inserting a hybrid charge intermediate state with modified partial charges, which resembles the charge distribution of the desired high level, makes it possible to obtain reliable calculations with 2 ps switches. Attempts using step-wise linear switching paths, on the other hand, did not lead to improvement, i.e., a faster convergence for all systems. To understand these findings, we analyzed the solutes' properties as a function of the partial charges used and the number of water molecules in direct contact with the solute, and studied the time needed for water molecules to reorient themselves upon a change in the solute's charge distribution.

Tautomer-Specific Deacylation and Ω-Loop Flexibility Explain the Carbapenem-Hydrolyzing Broad-Spectrum Activity of the KPC-2 β-Lactamase

KPC-2 (Klebsiella pneumoniae carbapenemase-2) is a globally disseminated serine-β-lactamase (SBL) responsible for extensive β-lactam antibiotic resistance in Gram-negative pathogens. SBLs inactivate β-lactams via a mechanism involving a hydrolytically labile covalent acyl-enzyme intermediate. Carbapenems, the most potent β-lactams, evade the activity of many SBLs by forming long-lived inhibitory acyl-enzymes; however, carbapenemases such as KPC-2 efficiently deacylate carbapenem acyl-enzymes. We present high-resolution (1.25-1.4 Å) crystal structures of KPC-2 acyl-enzymes with representative penicillins (ampicillin), cephalosporins (cefalothin), and carbapenems (imipenem, meropenem, and ertapenem) obtained utilizing an isosteric deacylation-deficient mutant (E166Q). The mobility of the Ω-loop (residues 165-170) negatively correlates with antibiotic turnover rates (kcat), highlighting the role of this region in positioning catalytic residues for efficient hydrolysis of different β-lactams. Carbapenem-derived acyl-enzyme structures reveal the predominance of the Δ1-(2R) imine rather than the Δ2 enamine tautomer. Quantum mechanics/molecular mechanics molecular dynamics simulations of KPC-2:meropenem acyl-enzyme deacylation used an adaptive string method to differentiate the reactivity of the two isomers. These identify the Δ1-(2R) isomer as having a significantly (7 kcal/mol) higher barrier than the Δ2 tautomer for the (rate-determining) formation of the tetrahedral deacylation intermediate. Deacylation is therefore likely to proceed predominantly from the Δ2, rather than the Δ1-(2R) acyl-enzyme, facilitated by tautomer-specific differences in hydrogen-bonding networks involving the carbapenem C-3 carboxylate and the deacylating water and stabilization by protonated N-4, accumulating a negative charge on the Δ2 enamine-derived oxyanion. Taken together, our data show how the flexible Ω-loop helps confer broad-spectrum activity upon KPC-2, while carbapenemase activity stems from efficient deacylation of the Δ2-enamine acyl-enzyme tautomer.

Structures and stabilities of PAH clusters solvated by water aggregates: The case of the pyrene dimer

Although clusters made of polycyclic aromatic hydrocarbon and water monomers are relevant objects in both atmospheric and astrophysical science, little is known about their energetic and structural properties. In this work, we perform global explorations of the potential energy landscapes of neutral clusters made of two pyrene units and one to ten water molecules using a density-functional-based tight-binding (DFTB) potential followed by local optimizations at the density-functional theory level. We discuss the binding energies with respect to various dissociation channels. It shows that cohesion energies of the water clusters interacting with a pyrene dimer are larger than those of the pure water clusters, reaching for the largest clusters an asymptotic limit similar to that of pure water clusters and that, although the hexamer and octamer can be considered magic numbers for isolated water clusters, it is not the case anymore when they are interacting with a pyrene dimer. Ionization potentials are also computed by making use of the configuration interaction extension of DFTB, and we show that in cations, the charge is mostly carried by the pyrene molecules.

Accelerating the density-functional tight-binding method using graphical processing units

Acceleration of the density-functional tight-binding (DFTB) method on single and multiple graphical processing units (GPUs) was accomplished using the MAGMA linear algebra library. Two major computational bottlenecks of DFTB ground-state calculations were addressed in our implementation: the Hamiltonian matrix diagonalization and the density matrix construction. The code was implemented and benchmarked on two different computer systems: (1) the SUMMIT IBM Power9 supercomputer at the Oak Ridge National Laboratory Leadership Computing Facility with 1-6 NVIDIA Volta V100 GPUs per computer node and (2) an in-house Intel Xeon computer with 1-2 NVIDIA Tesla P100 GPUs. The performance and parallel scalability were measured for three molecular models of 1-, 2-, and 3-dimensional chemical systems, represented by carbon nanotubes, covalent organic frameworks, and water clusters.

Graph-convolutional neural networks for (QM)ML/MM molecular dynamics simulations

To accurately study the chemical reactions in the condensed phase or within enzymes, both quantum-mechanical description and sufficient configurational sampling are required to reach converged estimates. Here, quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations play an important role, providing QM accuracy for the region of interest at a decreased computational cost. However, QM/MM simulations are still too expensive to study large systems on longer time scales. Recently, machine learning (ML) models have been proposed to replace the QM description. The main limitation of these models lies in the accurate description of long-range interactions present in condensed-phase systems. To overcome this issue, a recent workflow has been introduced combining a semi-empirical method (i.e. density functional tight binding (DFTB)) and a high-dimensional neural network potential (HDNNP) in a Δ-learning scheme. This approach has been shown to be capable of correctly incorporating long-range interactions within a cutoff of 1.4 nm. One of the promising alternative approaches to efficiently take long-range effects into account is the development of graph-convolutional neural networks (GCNNs) for the prediction of the potential-energy surface. In this work, we investigate the use of GCNN models - with and without a Δ-learning scheme - for (QM)ML/MM MD simulations. We show that the Δ-learning approach using a GCNN and DFTB as a baseline achieves competitive performance on our benchmarking set of solutes and chemical reactions in water. This method is additionally validated by performing prospective (QM)ML/MM MD simulations of retinoic acid in water and S-adenoslymethionine interacting with cytosine in water. The results indicate that the Δ-learning GCNN model is a valuable alternative for the (QM)ML/MM MD simulations of condensed-phase systems.

Thermal decomposition mechanisms of energetic CL-20-based co-crystals: quantum molecular dynamics simulations

The decomposition mechanisms of energetic CL-20:2,4-dinitro-2,4-diazapentane (DNP) and CL-20:2,4-dinitro-2,4-diazaheptane (DNG) co-crystals at high temperatures (1000, 2000, and 3000 K) were studied by density functional tight-binding molecular dynamics (DFTB-MD) simulations. At different temperatures, their decomposition mechanisms are very different. At 1000 K, conformational changes are observed only for the CL-20:DNG co-crystal, in which the CL-20 changes from β-CL-20 to γ-CL-20. When the temperature is increased to 2000 K, CL-20, DNP, and DNG begin to decompose, and there are five paths for the main initial mechanisms. Further increasing the temperature to 3000 K promotes a more complete decomposition. The initial reactions of CL-20 in the two co-crystals have two channels. There are two initial decomposition channels in the DNP molecule and only one channel in the DNG molecule. As the temperature increases, the decomposition products of the two co-crystals are different. Our work may provide the in-depth understanding of the decomposition mechanisms of high-energy CL-20-based co-crystals at high temperatures.

Efficient simulation of resonance Raman spectra with tight-binding approximations to Density Functional Theory

Resonance Raman spectroscopy has long been established as one of the most sensitive techniques for detection, structure characterization and probing the excited-state dynamics of biochemical systems. However, the analysis of resonance Raman spectra is much facilitated when measurements are accompanied by Density Functional Theory (DFT) calculations which are expensive for large biomolecules. In this work, resonance Raman spectra are therefore computed with the Density Functional Tight-Binding (DFTB) method in the time-dependent excited-state gradient approximation. To test the accuracy of the tight-binding approximations, this method is first applied to typical resonance Raman benchmark molecules like β-carotene and compared to results obtained with pure and range-separated exchange-correlation (xc) functionals. We then demonstrate the efficiency of the approach by considering a computationally challenging heme variation. Overall, we find that the vibrational frequencies and excited-state properties (energies and gradients) which are needed to simulate the spectra are reasonably accurate and suitable for interpretation of experiments. We can therefore recommend DFTB as a fast computational method to interpret resonance Raman spectra.

Comprehensive Study of the Chemistry behind the Stability of Carboxylic SWCNT Dispersions in the Development of a Transparent Electrode

Single-walled carbon nanotubes (SWCNTs) are well-known for their excellent electrical conductivity. One promising application for SWCNT-based thin films is as transparent electrodes for uncooled mid-IR detectors (MIR). In this paper, a combination of computational and experimental studies were performed to understand the chemistry behind the stability of carboxylic SWCNTs (SWCNTs-COOH) dispersions in different solvents. A computational study based on the density functional tight-binding (DFTB) method was applied to understand the interactions of COOH-functionalized carbon nanotubes with selected solvents. Attention was focused on understanding how the protonation of COOH groups influences the binding energies between SWCNTs and different solvents. Thin film electrodes were prepared by alternately depositing PEI and SWCNT-COOH on soda lime glass substrates. To prepare a stable SWCNT dispersion, different solvents were tested, such as deionized (DI) water, ethanol and acetone. The SWCNT-COOH dispersion stability was tested in different solvents. Samples were prepared to study the relationship between the number of depositions, transparency in the MIR range (2.5–5 µm) and conductivity, looking for the optimal thickness that would satisfy the application. The MIR transparency of the electrode was reduced by 20% for the thickest SWCNT layers, whereas sheet resistance values were reduced to 150–200 kΩ/sq.

QM/MM investigation of the catalytic mechanism of processive endoglucanase Cel9G from Clostridium cellulovorans

Carbohydrate degradation catalyzed by glucoside hydrolases (GHs) is a major mechanism in biomass conversion. GH family 9 endoglucanase (Cel9G) from Clostridium cellulovorans, a typical multimodular enzyme, contains a catalytic domain closely linked to a family 3c carbohydrate-binding module (CBM3c). Unlike the conventional behavior proposed for other carbohydrate-binding modules, CBM3c has a direct impact on catalytic activity. In this work, extensive molecular dynamics (MD) simulations were employed to clarify the functional role of CBM3c. Furthermore, the detailed catalytic mechanism of Cel9G was investigated at the atomistic level using the combined quantum mechanical and molecular mechanical (QM/MM) method. Based on these simulations, owing to the rigidity of the peptide linker, CBM3c may affect the enzymatic activity via direct interactions with alpha helix 4 of GH9, especially with the K123 and H125 residues. In addition, using cellohexaose as a substrate, the QM/MM MD simulations confirmed that this enzyme can cleave the β-1,4-glycosidic linkage via an inverting mechanism. An oxocarbenium ion-like transition state was located with a barrier height of 19.6 kcal mol^-1. Furthermore, the G(-1) pyranose unit preferentially adopted a distorted ¹S₅/⁴H₅ conformer in the enzyme-substrate complex. For the cleavage of the glycosidic bond, we were able to identify a plausible route (¹S₅/⁴H₅ → [⁴H₅/⁴E]^# → ⁴C₁) from the reactant to the product at the G(-1) site.

Mechanistic Modeling of Monoglyceride Lipase Covalent Modification Elucidates the Role of Leaving Group Expulsion and Discriminates Inhibitors with High and Low Potency

Inhibition of monoglyceride lipase (MGL), also known as monoacylglycerol lipase (MAGL), has emerged as a promising approach for treating neurological diseases. To gain useful insights in the design of agents with balanced potency and reactivity, we investigated the mechanism of MGL carbamoylation by the reference triazole urea SAR629 (IC₅₀ = 0.2 nM) and two recently described inhibitors featuring a pyrazole (IC₅₀ = 1800 nM) or a 4-cyanopyrazole (IC₅₀ = 8 nM) leaving group (LG), using a hybrid quantum mechanics/molecular mechanics (QM/MM) approach. Opposite to what was found for substrate 2-arachidonoyl-sn-glycerol (2-AG), covalent modification of MGL by azole ureas is controlled by LG expulsion. Simulations indicated that changes in the electronic structure of the LG greatly affect reaction energetics with triazole and 4-cyanopyrazole inhibitors following a more accessible carbamoylation path compared to the unsubstituted pyrazole derivative. The computational protocol provided reaction barriers able to discriminate between MGL inhibitors with different potencies. These results highlight how QM/MM simulations can contribute to elucidating structure–activity relationships and provide insights for the design of covalent inhibitors.

Theoretical Characterization of the Step-by-Step Mechanism of Conversion of Leukotriene A4 to Leukotriene B4 Catalysed by the Enzyme Leukotriene A4 Hydrolase

LTA₄H is a bifunctional zinc metalloenzyme that converts leukotriene A₄ (LTA₄) into leukotriene B₄ (LTB₄), one of the most potent chemotactic agents involved in acute and chronic inflammatory diseases. In this reaction, LTA₄H acts as an epoxide hydrolase with a unique and fascinating mechanism, which includes the stereoselective attachment of one water molecule to the carbon backbone of LTA₄ several methylene units away from the epoxide moiety. By combining Molecular Dynamics simulations and Quantum Mechanics/Molecular Mechanics calculations, we obtained a very detailed molecular picture of the different consecutive steps of that mechanism. By means of a rather unusual 1,7-nucleophilic substitution through a clear S_N1 mechanism, the epoxide opens and the triene moiety of the substrate twists in such a way that the bond C₆-C₇ adopts its cis (Z) configuration, thus exposing the R face of C₁₂ to the addition of a water molecule hydrogen-bonded to ASP375. Thus, the two stereochemical features that are required for the bioactivity of LTB₄ appear to be closely related. The noncovalent π-π stacking interactions between the triene moiety and two tyrosines (TYR267 and, especially, TYR378) that wrap the triene system along the whole reaction explain the preference for the cis configuration inside LTA₄H.

Absorption Properties of Large Complex Molecular Systems: The DFTB/Fluctuating Charge Approach

We report on the first formulation of a novel polarizable QM/MM approach, where the density functional tight binding (DFTB) is coupled to the fluctuating charge (FQ) force field. The resulting method (DFTB/FQ) is then extended to the linear response within the TD-DFTB framework and challenged to study absorption spectra of large condensed-phase systems.

Development and Evaluation of Geometry Optimization Algorithms in Conjunction with ANI Potentials

An efficient yet accurate method for producing a large amount of energy data for molecular mechanical force field (MMFF) parameterization is on demand, especially for torsional angle parameters which are typically derived to reproduce ab initio rotational profiles or torsional potential energy surfaces (PESs). Recently, an active learning potential (ANI-1x) for organic molecules which can produce smooth and physically meaningful PESs has been developed. The high efficiency and accuracy make ANI-1x especially attractive for geometry optimization at low cost. To apply the ANI-1x potential in MMFF parameterization, one needs to perform constrained geometry optimization. In this work, we first developed a computational protocol to constrain rotatable torsional angles and other geometric parameters for a molecule whose geometry is described by Cartesian coordinates. The constraint is successfully achieved by force projection for the two conjugated gradient (CG) algorithms. We then conducted large-scale assessments on ANI-1x along with four different optimization algorithms in reproducing DFT energies and geometries for two CG algorithms, CG backtracking line search (CG-BS) and CG Wolfe line search (CG-WS), and two quasi-Newton algorithms, Broyden-Fletcher-Goldfarb-Shanno (BFGS) and low-memory BFGS (L-BFGS). Note that CG-BS is a new algorithm we developed in this work. All four algorithms take the ANI energies and forces to optimize a molecule geometry. Last, we conducted a large-scale assessment of applying ANI-1x in MMFF development in three aspects. First, we performed full optimizations for 100 drug molecules, each consisting of five distinct conformations. The average root-mean-square error (RMSE) between ANI-1x and DFT is about 1.3 kcal/mol, and the root-mean-square displacement (RMSD) of heavy atoms is about 0.35 Å. Second, we generated torsional PESs for 160 organic molecules, and constrained optimizations were performed for up to 18 conformations for each PES. We found that the RMSE of all the conformers is 1.23 kcal/mol. Last, we carried out constrained optimizations for alanine dipeptide with both ϕ and φ angles being frozen. The Ramachandran plots indicate that the two CG algorithms in conjunction with the ANI-1x potential could well reproduce the DFT-optimized geometries and torsional PESs. We concluded that CG-BS and CG-WS are good choices for generating PESs, while CG-WS or BFGS is ideal for performing full geometry optimization. With the continuously increased quality of ANI, it is expected that the computational algorithms and protocols presented in this work will have great applications in improving the quality of an existing small-molecule MMFF.

Structural stability and initial decomposition mechanisms of BTF crystal induced by vacancy defects: a computational study

Density functional tight binding (DFTB) and DFTB-based molecular dynamics (DFTB-MD) were used to study the effects of vacancy defects on the structure, stability, and initial decomposition mechanisms of condensed phase...

Generation of Quantum Configurational Ensembles Using Approximate Potentials

Conformational analysis is of paramount importance in drug design: it is crucial to determine pharmacological properties, understand molecular recognition processes, and characterize the conformations of ligands when unbound. Molecular Mechanics (MM) simulation methods, such as Monte Carlo (MC) and molecular dynamics (MD), are usually employed to generate ensembles of structures due to their ability to extensively sample the conformational space of molecules. The accuracy of these MM-based schemes strongly depends on the functional form of the force field (FF) and its parametrization, components that often hinder their performance. High-level methods, such as ab initio MD, provide reliable structural information but are still too computationally expensive to allow for extensive sampling. Therefore, to overcome these limitations, we present a multilevel MC method that is capable of generating quantum configurational ensembles while keeping the computational cost at a minimum. We show that FF reparametrization is an efficient route to generate FFs that reproduce QM results more closely, which, in turn, can be used as low-cost models to achieve the gold standard QM accuracy. We demonstrate that the MC acceptance rate is strongly correlated with various phase space overlap measurements and that it constitutes a robust metric to evaluate the similarity between the MM and QM levels of theory. As a more advanced application, we present a self-parametrizing version of the algorithm, which combines sampling and FF parametrization in one scheme, and apply the methodology to generate the QM/MM distribution of a ligand in aqueous solution.

Surface tension of liquids and binary mixtures from molecular dynamics simulations

In this work we assess and extend strategies for calculating surface tension of complex liquids from molecular dynamics simulations: the mechanical route and the instantaneous liquid interface (ILI) approach. The former employs the connection between stress tensor and surface tension, whereas the latter involves computation of instantaneous density field. Whereas the mechanical route is general, the ILI method involves system-dependent parameters restricting its original application to liquid water only. Here we generalize the approach to complex molecular liquids using atomic van der Waals radii. The performance of the approaches is evaluated on two liquid systems: acetonitrile and water-methanol mixture. In addition, we compare the effect of the computational models for interaction potentials based on semi-empirical electronic structure theory and classical force fields on the estimate of the surface tension within both stress tensor and ILI approaches.

Insights into Coke Formation and Removal under Operating Conditions with a Quantum Nanoreactor Approach

The in situ formation and removal of coke is a critical problem in heterogeneous catalysis, but its mechanism is not well understood. This work investigates the mechanism of carbon deposition and hydrogenation on an Fe cluster under high-temperature conditions with the density functional tight-binding (DFTB) based nanoreactor molecular dynamics (NMD) method. Our study shows that successive formation of carbon chains, rings, and fused rings occurred during the carbon deposition on Fe clusters. Hydrogenation of activated carbon happens through direct C-H coupling, while the hydrogenation of graphitic carbon involves hydrogenation of the edge carbon, ring-opening reaction, and dealkylation reaction. The main function of the Fe catalyst is to provide the active sites for H₂ dissociation and dissociated H spillover, while its activity toward C-C bond breaking is limited. These results highlight the role of the DFTB-NMD method as an effective tool to investigate reaction mechanisms under operating conditions in heterogeneous catalysis.

Structural Factors Determining the Absorption Spectrum of Channelrhodopsins: A Case Study of the Chimera C1C2

Channelrhodopsins are photosensitive proteins that trigger flagella motion in single-cell algae and have been successfully utilized in optogenetic applications. In optogenetics, light is used to activate neural cells in living organisms, which can be achieved by exploiting the ion channel signaling of channelrhodopsins. Tailoring channelrhodopsins for such applications includes the tuning of the absorption maximum. In order to establish rational design and to obtain a desired spectral shift, a basic understanding of the absorption spectrum is required. We have studied the chimera C1C2 as a representative of this protein family and the first member with an available crystal structure. For this purpose, we sampled the conformations of C1C2 using quantum mechanical/molecular mechanical molecular dynamics and subjected the resulting snapshots of the trajectory to excitation energy calculations using ADC(2) and simplified time-dependent density functional theory. In contrast to previous reports, we found that different hydrogen-bonding networks-involving the retinal protonated Schiff base, the putative counterions E162 and D292, and water molecules-had only a small impact on the absorption spectrum. However, in the case of deprotonated E162, increasing the distance to the Schiff base hydrogen-bonding partner led to a systematic blue shift. The β-ionone ring rotation was identified as another important contributor. Yet the most important factors were found to be the bond length alternation and bond order alternation that were linearly correlated to the absorption maximum by up to 62 and 82%, respectively. We ascribe this novel insight into the structural basis of the absorption spectrum to our enhanced protein setup that includes membrane embedding as well as long and extensive sampling.

Modelling biocompatible ionic liquids based on organic acids and amino acids: challenges for computational models and future perspectives

In this short review I shall highlight the basic principle and the difficulties that arise in attempting the computational modeling of seemingly simple systems which hide an unexpected complexity. Biocompatible ionic liquids which are based on the coupling of organic or amino acid anions with metabolic cations such as cholinium are the target of this review. These substances have been the subject of intense research activities in the last few years and have attracted the attention of computational chemists. I shall show that the computational description of these substances is far from trivial and requires the use of sophisticated techniques in order to account for a surprisingly rich chemistry that is due to several phenomena such as polarization, charge transfer, proton transfer equilibria and tautomerization reactions.

Docking Paradigm in Drug Design

Docking is in demand for the rational computer aided structure based drug design. A review of docking methods and programs is presented. Different types of docking programs are described. They include docking of non-covalent small ligands, protein-protein docking, supercomputer docking, quantum docking, the new generation of docking programs and the application of docking for covalent inhibitors discovery. Taking into account the threat of COVID-19, we present here a short review of docking applications to the discovery of inhibitors of SARS-CoV and SARS-CoV-2 target proteins, including our own result of the search for inhibitors of SARS-CoV-2 main protease using docking and quantum chemical post-processing. The conclusion is made that docking is extremely important in the fight against COVID-19 during the process of development of antivirus drugs having a direct action on SARS-CoV-2 target proteins.

Cis-Trans Isomerization and Thermal Decomposition Mechanisms of a Series of Nx (x = 4, 8, 10, 11) Chain-Catenated Energetic Crystals

Nitrogen-rich compounds based on heteroaromatic rings with different lengths of nitrogen chains are at the forefront of the energetic materials field. We studied the decomposition processes and reaction kinetics of a series of N_x (x = 4, 8, 10, 11) chain-catenated energetic crystals at various temperatures (2400-3000 K) based on a combinational strategy based on density functional tight binding molecular dynamics (DFTB-MD) simulations and density functional theory (DFT). The results show that the thermal decomposition and reaction kinetics are dependent on both the temperature and nitrogen chain's length. There are two sequential stages in the initial decomposition process for the crystals N₈ and N₁₀: (i) competition between cis-trans isomerization and initial unimolecular decomposition and (ii) subsequent complicated global decomposition reactions. Increasing either the temperature or nitrogen chain's length will accelerate the competition and make initial decomposition dominate. However, cis-trans isomerization does not occur in the crystals N₄ and N₁₁. The dominant initiation paths for N₄, N₈, and N₁₀ occur in the heterocycle and in the bond between the heterocycle and azo group, while that for N₁₁ is ring elimination. The decomposition reactions exhibit a clear first-order kinetics character. The energy paths based on DFT calculations are determined as an addition to the DFTB-MD results. Our findings provide insights into the comprehensive understanding of thermal decomposition behaviors of nitrogen chain-catenated and even all-nitrogen energetic materials.

Conformational Shifts of Stacked Heteroaromatics: Vacuum vs. Water Studied by Machine Learning

Stacking interactions play a crucial role in drug design, as we can find aromatic cores or scaffolds in almost any available small molecule drug. To predict optimal binding geometries and enhance stacking interactions, usually high-level quantum mechanical calculations are performed. These calculations have two major drawbacks: they are very time consuming, and solvation can only be considered using implicit solvation. Therefore, most calculations are performed in vacuum. However, recent studies have revealed a direct correlation between the desolvation penalty, vacuum stacking interactions and binding affinity, making predictions even more difficult. To overcome the drawbacks of quantum mechanical calculations, in this study we use neural networks to perform fast geometry optimizations and molecular dynamics simulations of heteroaromatics stacked with toluene in vacuum and in explicit solvation. We show that the resulting energies in vacuum are in good agreement with high-level quantum mechanical calculations. Furthermore, we show that using explicit solvation substantially influences the favored orientations of heteroaromatic rings thereby emphasizing the necessity to include solvation properties starting from the earliest phases of drug design.

Machine learning-accelerated quantum mechanics-based atomistic simulations for industrial applications

Atomistic simulations have become an invaluable tool for industrial applications ranging from the optimization of protein-ligand interactions for drug discovery to the design of new materials for energy applications. Here we review recent advances in the use of machine learning (ML) methods for accelerated simulations based on a quantum mechanical (QM) description of the system. We show how recent progress in ML methods has dramatically extended the applicability range of conventional QM-based simulations, allowing to calculate industrially relevant properties with enhanced accuracy, at reduced computational cost, and for length and time scales that would have otherwise not been accessible. We illustrate the benefits of ML-accelerated atomistic simulations for industrial R&D processes by showcasing relevant applications from two very different areas, drug discovery (pharmaceuticals) and energy materials. Writing from the perspective of both a molecular and a materials modeling scientist, this review aims to provide a unified picture of the impact of ML-accelerated atomistic simulations on the pharmaceutical, chemical, and materials industries and gives an outlook on the exciting opportunities that could emerge in the future.

ParaMol: A Package for Automatic Parameterization of Molecular Mechanics Force Fields

The ensemble of structures generated by molecular mechanics (MM) simulations is determined by the functional form of the force field employed and its parameterization. For a given functional form, the quality of the parameterization is crucial and will determine how accurately we can compute observable properties from simulations. While accurate force field parameterizations are available for biomolecules, such as proteins or DNA, the parameterization of new molecules, such as drug candidates, is particularly challenging as these may involve functional groups and interactions for which accurate parameters may not be available. Here, in an effort to address this problem, we present ParaMol, a Python package that has a special focus on the parameterization of bonded and nonbonded terms of druglike molecules by fitting to ab initio data. We demonstrate the software by deriving bonded terms' parameters of three widely known drug molecules, viz. aspirin, caffeine, and a norfloxacin analogue, for which we show that, within the constraints of the functional form, the methodologies implemented in ParaMol are able to derive near-ideal parameters. Additionally, we illustrate the best practices to follow when employing specific parameterization routes. We also determine the sensitivity of different fitting data sets, such as relaxed dihedral scans and configurational ensembles, to the parameterization procedure, and discuss the features of the various weighting methods available to weight configurations. Owing to ParaMol's capabilities, we propose that this software can be introduced as a routine step in the protocol normally employed to parameterize druglike molecules for MM simulations.

Analytical Time-Dependent Long-Range Corrected Density Functional Tight Binding (TD-LC-DFTB) Gradients in DFTB+: Implementation and Benchmark for Excited-State Geometries and Transition Energies

The absorption and emission of light is a ubiquitous process in chemical and biological processes, making a theoretical description inevitable for understanding and predicting such properties. Although ab initio and DFT methods are capable of describing excited states with good accuracy in many cases, the investigation of dynamical processes and the need to sample the phase space in complex systems often requires methods with reduced computational costs but still sufficient accuracy. In the present work, we report the derivation and implementation of analytical nuclear gradients for time-dependent long-range corrected density functional tight binding (TD-LC-DFTB) in the DFTB+ program. The accuracy of the TD-LC-DFTB potential-energy surfaces is benchmarked for excited-state geometries and adiabatic as well as vertical transition energies. The benchmark set consists of more than 100 organic molecules taken as subsets from available benchmark sets. The reported method yields a mean deviation of 0.31 eV for adiabatic excitation energies with respect to CC2. In order to study more subtle effects, seminumerical second derivatives based on the analytical gradients are employed to simulate vibrationally resolved UV/vis spectra. This extensive test exhibits few problematic cases, which can be traced back to the parametrization of the repulsive potential.

Exploring of paritaprevir and glecaprevir resistance due to A156T mutation of HCV NS3/4A protease: molecular dynamics simulation study

Hepatitis C virus (HCV) NS3/4A serine protease is a promising drug target for the discovery of anti-HCV drugs. However, its amino acid mutations, particularly A156T, commonly lead to rapid emergence of drug resistance. Paritaprevir and glecaprevir, the newly FDA-approved HCV drugs, exhibit distinct resistance profiles against the A156T mutation of HCV NS3/4A serine protease. To illustrate their different molecular resistance mechanisms, molecular dynamics simulations and binding free energy calculations were carried out on the two compounds complexed with both wild-type (WT) and A156T variants of HCV NS3/4A protease. QM/MM-GBSA-based binding free energy calculations revealed that the binding affinities of paritaprevir and glecaprevir towards A156T NS3/4A were significantly reduced by ∼4 kcal/mol with respect to their WT complexes, which were in line with the experimental resistance folds. Moreover, the relatively weak intermolecular interactions with amino acids such as H57, R155, and T156 of NS3 protein, the steric effect and the destabilized protein binding surface, which is caused by the loss of salt bridge between R123 and D168, are the main contributions for the higher fold-loss in potency of glecaprevir due to A156T mutation. An insight into the difference of molecular mechanism of drug resistance against the A156T substitution among the two classes of serine protease inhibitors could be useful for further optimization of new generation HCV NS3/4A inhibitors with enhanced inhibitory potency.Communicated by Ramaswamy H. Sarma.

Structure, intermolecular interactions, and dynamic properties of NTO crystals with impurity defects: a computational study

Frontier orbitals distribute in the position of impurity molecules, whose adjacent NTO molecules begin to decompose first.

Effect of pore-size disorder on the electronic properties of semiconducting graphene nanomeshes

Graphene nanomeshes (GNMs) are novel materials that recently raised a lot of interest. They are fabricated by forming a lattice of pores in graphene. Depending on the pore size and pore lattice constant, GNMs can be either semimetallic or semiconducting with a gap large enough (∼ 0.5 eV) to be considered for transistor applications. The fabrication process is bound to produce some structural disorder due to variations in pore size. Recent electronic transport measurements in GNM devices (ACS Appl. Mater. Interfaces 10, 10 362, 2018) show a degradation of their bandgap in devices having pore-size disorder. It is therefore important to understand the effect of such variability on the electronic properties of semiconducting GNMs. In this work we use the density functional-based tight binding formalism to calculate the electronic properties of GNM structures with different pore sizes, pore densities, and with hydrogen and oxygen pore edge passivations. We find that structural disorder reduces the electronic gap and the carrier group velocity, which may interpret recent transport measurements in GNM devices. Furthermore, the trend of the bandgap with structural disorder is not significantly affected by the change in pore edge passivation. Our results show that even with structural disorder, GNMs are still attractive from a transistor device perspective.

Adventures in DFTB: Toward an Automatic Parameterization Scheme

As we push forward on understanding the fate of chemicals in the environment, we need a method that will allow for the simulation of the inherent heterogeneity. Density functional tight binding (DFTB) is a methodology that allows for a detailed electronic description and would be ideal for this problem. While many parameters can be derived directly from DFT, empirical parameters still exist in the confinement and repulsion potentials. In this manuscript, we examine these potentials and present solutions that will minimize the degree of empiricism. Our results show that it is possible to construct confinement potentials from examining the atomic radial wavefunctions. Moreover, we found that the heterogeneous repulsion potentials can be derived from using only homogeneous repulsion curves.

Mechanistic Insights into the Hydrolysis of O-GlcNAcylation Catalyzed by Human O-GlcNAcase

O-GlcNAc glycosylation occurs on specific serine/threonine residues of intracellular proteins, which is widely related to various diseases, including type II diabetes, cancer, and Alzheimer's disease. Human O-GlcNAcase (hOGA) is responsible for the removal of O-GlcNAc modification and thus serves as the main target for inhibitor design. In this work, we systematically investigated the mechanism catalyzed by hOGA using the quantum mechanical/molecular mechanical method. Based on calculated free energy profiles, two essential steps named cyclization (Cyc) step and ring opening step are required to generate the final hemiacetal product. The Cyc of the 2-acetamido group, the rate-limiting step, leads to the generation of the intermediate of a bicyclic oxazolinium ion (EI1). Next, the oxazoline ring could be broken via the nucleophilic attack of a water molecule at the C1 position, which generates the final product. Along with this, our simulations clearly suggest the existence of an oxazoline intermediate (EI2), which is produced via proton transfer (PT) from the 2-acetamido group (EI1) to D174. This PT step features a reversible process with a low energy barrier, which could be attributed to a low barrier hydrogen bond between the donor and acceptor. The stabilizing effect of the low barrier hydrogen bond on EI1 is proposed to be very important for accelerating the overall reaction. In fact, the site-directed mutagenesis simulations of D174A and D175A strongly indicate that the catalytic residues mainly affect the observed reaction rate by affecting the stability of the intermediate.