Direct calculation of electron density in density-functional theory

Excitation configuration analysis for divide-and-conquer excited-state calculation method using dynamical polarizability

The authors previously developed a divide-and-conquer (DC)-based non-local excited-state calculation method for large systems using dynamical polarizability [Nakai and Yoshikawa, J. Chem. Phys. 146, 124123 (2017)]. This method evaluates the excitation energies and oscillator strengths using information on the dynamical polarizability poles. This article proposes a novel analysis of the previously developed method to obtain further configuration information on excited states, including excitation and de-excitation coefficients of each excitation configuration. Numerical applications to simple molecules, such as ethylene, hydrogen molecule, ammonia, and pyridazine, confirmed that the proposed analysis could accurately reproduce the excitation and de-excitation coefficients. The combination with the DC scheme enables both the local and non-local excited states of large systems with an excited nature to be treated.

The development of the QM/MM interface and its application for the on-the-fly QM/MM nonadiabatic dynamics in JADE package: Theory, implementation, and applications

Understanding the nonadiabatic dynamics of complex systems is a challenging task in computational photochemistry. Herein, we present an efficient and user-friendly quantum mechanics/molecular mechanics (QM/MM) interface to run on-the-fly nonadiabatic dynamics. Currently, this interface consists of an independent set of codes designed for general-purpose use. Herein, we demonstrate the ability and feasibility of the QM/MM interface by integrating it with our long-term developed JADE package. Tailored to handle nonadiabatic processes in various complex systems, especially condensed phases and protein environments, we delve into the theories, implementations, and applications of on-the-fly QM/MM nonadiabatic dynamics. The QM/MM approach is established within the framework of the additive QM/MM scheme, employing electrostatic embedding, link-atom inclusion, and charge-redistribution schemes to treat the QM/MM boundary. Trajectory surface-hopping dynamics are facilitated using the fewest switches algorithm, encompassing classical and quantum treatments for nuclear and electronic motions, respectively. Finally, we report simulations of nonadiabatic dynamics for two typical systems: azomethane in water and the retinal chromophore PSB3 in a protein environment. Our results not only illustrate the power of the QM/MM program but also reveal the important roles of environmental factors in nonadiabatic processes.

Noise reduction of stochastic density functional theory for metals

Density Functional Theory (DFT) has become a cornerstone in the modeling of metals. However, accurately simulating metals, particularly under extreme conditions, presents two significant challenges. First, simulating complex metallic systems at low electron temperatures is difficult due to their highly delocalized density matrix. Second, modeling metallic warm-dense materials at very high electron temperatures is challenging because it requires the computation of a large number of partially occupied orbitals. This study demonstrates that both challenges can be effectively addressed using the latest advances in linear-scaling stochastic DFT methodologies. Despite the inherent introduction of noise into all computed properties by stochastic DFT, this research evaluates the efficacy of various noise reduction techniques under different thermal conditions. Our observations indicate that the effectiveness of noise reduction strategies varies significantly with the electron temperature. Furthermore, we provide evidence that the computational cost of stochastic DFT methods scales linearly with system size for metal systems, regardless of the electron temperature regime.

Dual functionality of the BiN monolayer: unraveling its photocatalytic and piezocatalytic water splitting properties

To achieve scalable and economically viable green hydrogen (H2) production, the photocatalytic and piezocatalytic processes are promising methods. The key to successful overall water splitting (OWS) for H2 production in these processes is using suitable semiconductor catalysts with appropriate band edge potentials, efficient optical absorption, higher mechanical flexibility, and piezoelectric coefficients. Thus, we explore the bismuth nitride (BiN) monolayer using density functional theory simulations, revealing intriguing catalytic properties. The BiN monolayer is a semiconductor with an indirect electronic bandgap (Eg) of 2.08 eV and displays excellent visible light absorption (approximately 105 cm-1). Detailed analyses show that the band edges satisfy the redox potential for photocatalytic OWS via biaxial strain engineering and pH variation. Notably, the solar to hydrogen conversion efficiency (ηSTH) for the BiN monolayer can reach 17.18%, which exceeds the 10% efficiency limit of photocatalysts for economical green H2 production. The obtained in-plane piezoelectric coefficient of e11 = 16.18 Å C m-1 is superior to widely studied 2D materials. Moreover, the generated piezopotential under oscillatory strain stands at 28.34 V, which can initiate the water redox reaction via the piezocatalytic mechanism. This originates from the mechanical flexibility coupled with higher piezoelectric coefficients. The result highlights the BiN monolayer's potential application in photocatalytic, piezocatalytic, and photo-piezo-catalytic OWS.

Surface-Enhanced Raman Spectroscopy: Current Understanding, Challenges, and Opportunities

While surface-enhanced Raman spectroscopy (SERS) has experienced substantial advancements since its discovery in the 1970s, it is an opportunity to celebrate achievements, consider ongoing endeavors, and anticipate the future trajectory of SERS. In this perspective, we encapsulate the latest breakthroughs in comprehending the electromagnetic enhancement mechanisms of SERS, and revisit CT mechanisms of semiconductors. We then summarize the strategies to improve sensitivity, selectivity, and reliability. After addressing experimental advancements, we comprehensively survey the progress on spectrum-structure correlation of SERS showcasing their important role in promoting SERS development. Finally, we anticipate forthcoming directions and opportunities, especially in deepening our insights into chemical or biological processes and establishing a clear spectrum-structure correlation.

Deep Learning-Assisted Spectrum-Structure Correlation: State-of-the-Art and Perspectives

Spectrum-structure correlation is playing an increasingly crucial role in spectral analysis and has undergone significant development in recent decades. With the advancement of spectrometers, the high-throughput detection triggers the explosive growth of spectral data, and the research extension from small molecules to biomolecules accompanies massive chemical space. Facing the evolving landscape of spectrum-structure correlation, conventional chemometrics becomes ill-equipped, and deep learning assisted chemometrics rapidly emerges as a flourishing approach with superior ability of extracting latent features and making precise predictions. In this review, the molecular and spectral representations and fundamental knowledge of deep learning are first introduced. We then summarize the development of how deep learning assist to establish the correlation between spectrum and molecular structure in the recent 5 years, by empowering spectral prediction (i.e., forward structure-spectrum correlation) and further enabling library matching and de novo molecular generation (i.e., inverse spectrum-structure correlation). Finally, we highlight the most important open issues persisted with corresponding potential solutions. With the fast development of deep learning, it is expected to see ultimate solution of establishing spectrum-structure correlation soon, which would trigger substantial development of various disciplines.

Long-range proton and hydroxide ion transfer dynamics at the water/CeO2 interface in the nanosecond regime: reactive molecular dynamics simulations and kinetic analysis

The structural properties, dynamical behaviors, and ion transport phenomena at the interface between water and cerium oxide are investigated by reactive molecular dynamics (MD) simulations employing neural network potentials (NNPs). The NNPs are trained to reproduce density functional theory (DFT) results, and DFT-based MD (DFT-MD) simulations with enhanced sampling techniques and refinement schemes are employed to efficiently and systematically acquire training data that include diverse hydrogen-bonding configurations caused by proton hopping events. The water interfaces with two low-index surfaces of (111) and (110) are explored with these NNPs, and the structure and long-range proton and hydroxide ion transfer dynamics are examined with unprecedented system sizes and long simulation times. Various types of proton hopping events at the interface are categorized and analyzed in detail. Furthermore, in order to decipher the proton and hydroxide ion transport phenomena along the surface, a counting analysis based on the semi-Markov process is formulated and applied to the MD trajectories to obtain reaction rates by considering the transport as stochastic jump processes. Through this model, the coupling between hopping events, vibrational motions, and hydrogen bond networks at the interface are quantitatively examined, and the high activity and ion transport phenomena at the water/CeO2 interface are unequivocally revealed in the nanosecond regime.

Fragment-Based Deep Learning for Simultaneous Prediction of Polarizabilities and NMR Shieldings of Macromolecules and Their Aggregates

Simultaneous prediction of the molecular response properties, such as polarizability and the NMR shielding constant, at a low computational cost is an unresolved issue. We propose to combine a linear-scaling generalized energy-based fragmentation (GEBF) method and deep learning (DL) with both molecular and atomic information-theoretic approach (ITA) quantities as effective descriptors. In GEBF, the total molecular polarizability can be assembled as a linear combination of the corresponding quantities calculated from a set of small embedded subsystems in GEBF. In the new GEBF-DL(ITA) protocol, one can predict subsystem polarizabilities based on the corresponding molecular wave function (thus electron density and ITA quantities) and DL model rather than calculate them from the computationally intensive coupled-perturbed Hartree-Fock or Kohn-Sham equations and finally obtain the total molecular polarizability via a linear combination equation. As a proof-of-concept application, we predict the molecular polarizabilities of large proteins and protein aggregates. GEBF-DL(ITA) is shown to be as accurate enough as GEBF, with mean absolute percentage error <1%. For the largest protein aggregate (>4000 atoms), GEBF-DL(ITA) gains a speedup ratio of 3 compared with GEBF. It is anticipated that when more advanced electronic structure methods are used, this advantage will be more appealing. Moreover, one can also predict the NMR chemical shieldings of proteins with reasonably good accuracy. Overall, the cost-efficient GEBF-DL(ITA) protocol should be a robust theoretical tool for simultaneously predicting polarizabilities and NMR shieldings of large systems.

CoeffNet: predicting activation barriers through a chemically-interpretable, equivariant and physically constrained graph neural network

Activation barriers of elementary reactions are essential to predict molecular reaction mechanisms and kinetics. However, computing these energy barriers by identifying transition states with electronic structure methods (e.g., density functional theory) can be time-consuming and computationally expensive. In this work, we introduce CoeffNet, an equivariant graph neural network that predicts activation barriers using coefficients of any frontier molecular orbital (such as the highest occupied molecular orbital) of reactant and product complexes as graph node features. We show that using coefficients as features offer several advantages, such as chemical interpretability and physical constraints on the network's behaviour and numerical range. Model outputs are either activation barriers or coefficients of the chosen molecular orbital of the transition state; the latter quantity allows us to interpret the results of the neural network through chemical intuition. We test CoeffNet on a dataset of SN2 reactions as a proof-of-concept and show that the activation barriers are predicted with a mean absolute error of less than 0.025 eV. The highest occupied molecular orbital of the transition state is visualized and the distribution of the orbital densities of the transition states is described for a few prototype SN2 reactions.

Robustness of Local Predictions in Atomistic Machine Learning Models

Machine learning (ML) models for molecules and materials commonly rely on a decomposition of the global target quantity into local, atom-centered contributions. This approach is convenient from a computational perspective, enabling large-scale ML-driven simulations with a linear-scaling cost and also allows for the identification and posthoc interpretation of contributions from individual chemical environments and motifs to complicated macroscopic properties. However, even though practical justifications exist for the local decomposition, only the global quantity is rigorously defined. Thus, when the atom-centered contributions are used, their sensitivity to the training strategy or the model architecture should be carefully considered. To this end, we introduce a quantitative metric, which we call the local prediction rigidity (LPR), that allows one to assess how robust the locally decomposed predictions of ML models are. We investigate the dependence of the LPR on the aspects of model training, particularly the composition of training data set, for a range of different problems from simple toy models to real chemical systems. We present strategies to systematically enhance the LPR, which can be used to improve the robustness, interpretability, and transferability of atomistic ML models.

Rapidly determining the 3D structure of proteins by surface-enhanced Raman spectroscopy

Despite great advances in protein structure analysis, label-free and ultrasensitive methods to obtain the natural and dynamic three-dimensional (3D) structures are still urgently needed. Surface-enhanced Raman spectroscopy (SERS) can be a good candidate, whereas the complexity originated from the interactions between the protein and the gradient surface electric field makes it extremely challenging to determine the protein structure. Here, we propose a deciphering strategy for accurate determination of 3D protein structure from experimental SERS spectra in seconds by simply summing SERS spectra of isolated amino acids in electric fields of different strength with their orientations in protein. The 3D protein structure can be reconstructed by comparing the experimental spectra obtained in a well-defined gap-mode SERS configuration with the simulated spectra. The gradient electric field endows SERS with a unique advantage to section biomolecules with atomic precision, which makes SERS a competent tool for monitoring biomolecular events under physiological conditions.

Gas phase H+, H3O+ and NH4+ affinities of oxygen-bearing volatile organic compounds; DFT calculations for soft chemical ionisation mass spectrometry

Quantum chemistry calculations were performed using the density functional theory, DFT, to understand the structures and energetics of organic ions relevant to gas phase ion chemistry in soft chemical ionisation mass spectrometry analytical methods. Geometries of a range of neutral volatile organic compound molecules and ions resulting from protonation, the addition of H3O+ and the addition of NH4+ were optimised using the B3LYP hybrid DFT method. Then, the total energies and the normal mode vibrational frequencies were determined, and the total enthalpies of the neutral molecules and ions were calculated for the standard temperature and pressure. The calculations were performed for several feasible structures of each of the ions. The proton affinities of several benchmark molecules agree with the accepted values within ±4 kJ mol-1, indicating that B3LYP/6-311++G(d,p) provides chemical accuracy for oxygen-containing volatile organic compounds. It was also found that the binding energies of H3O+ and NH4+ to molecules correlate with their proton affinities. The results contribute to the understanding of ligand switching ion-molecule reactions important for secondary electrospray ionisation, SESI, and selected ion flow tube, SIFT, mass spectrometries.

Spectral-partitioned Kohn-Sham density functional theory

We introduce a general, variational scheme for systematic approximation of a given Kohn-Sham free-energy functional by partitioning the density matrix into distinct spectral domains, each of which may be spanned by an independent diagonal representation without requirement of mutual orthogonality. It is shown that by generalizing the entropic contribution to the free energy to allow for independent representations in each spectral domain, the free energy becomes an upper bound to the exact (unpartitioned) Kohn-Sham free energy, attaining this limit as the representations approach Kohn-Sham eigenfunctions. A numerical procedure is devised for calculation of the generalized entropy associated with spectral partitioning of the density matrix. The result is a powerful framework for Kohn-Sham calculations of systems whose occupied subspaces span multiple energy regimes. As a case in point, we apply the proposed framework to warm- and hot-dense matter described by finite-temperature density functional theory, where at high energies the density matrix is represented by that of the free-electron gas, while at low energies it is variationally optimized. We derive expressions for the spectral-partitioned Kohn-Sham Hamiltonian, atomic forces, and macroscopic stresses within the projector-augmented wave (PAW) and the norm-conserving pseudopotential methods. It is demonstrated that at high temperatures, spectral partitioning facilitates accurate calculations at dramatically reduced computational cost. Moreover, as temperature is increased, fewer exact Kohn-Sham states are required for a given accuracy, leading to further reductions in computational cost. Finally, it is shown that standard multiprojector expansions of electronic orbitals within atomic spheres in the PAW method lack sufficient completeness at high temperatures. Spectral partitioning provides a systematic solution for this fundamental problem.

Toward a general neural network force field for protein simulations: Refining the intramolecular interaction in protein

Molecular dynamics (MD) is an extremely powerful, highly effective, and widely used approach to understanding the nature of chemical processes in atomic details for proteins. The accuracy of results from MD simulations is highly dependent on force fields. Currently, molecular mechanical (MM) force fields are mainly utilized in MD simulations because of their low computational cost. Quantum mechanical (QM) calculation has high accuracy, but it is exceedingly time consuming for protein simulations. Machine learning (ML) provides the capability for generating accurate potential at the QM level without increasing much computational effort for specific systems that can be studied at the QM level. However, the construction of general machine learned force fields, needed for broad applications and large and complex systems, is still challenging. Here, general and transferable neural network (NN) force fields based on CHARMM force fields, named CHARMM-NN, are constructed for proteins by training NN models on 27 fragments partitioned from the residue-based systematic molecular fragmentation (rSMF) method. The NN for each fragment is based on atom types and uses new input features that are similar to MM inputs, including bonds, angles, dihedrals, and non-bonded terms, which enhance the compatibility of CHARMM-NN to MM MD and enable the implementation of CHARMM-NN force fields in different MD programs. While the main part of the energy of the protein is based on rSMF and NN, the nonbonded interactions between the fragments and with water are taken from the CHARMM force field through mechanical embedding. The validations of the method for dipeptides on geometric data, relative potential energies, and structural reorganization energies demonstrate that the CHARMM-NN local minima on the potential energy surface are very accurate approximations to QM, showing the success of CHARMM-NN for bonded interactions. However, the MD simulations on peptides and proteins indicate that more accurate methods to represent protein-water interactions in fragments and non-bonded interactions between fragments should be considered in the future improvement of CHARMM-NN, which can increase the accuracy of approximation beyond the current mechanical embedding QM/MM level.

Theoretical and Numerical Comparison of Quantum- and Classical Embedding Models for Optical Spectra

Quantum-mechanical (QM) and classical embedding models approximate a supermolecular quantum-chemical calculation. This is particularly useful when the supermolecular calculation has a size that is out of reach for present QM models. Although QM and classical embedding methods share the same goal, they approach this goal from different starting points. In this study, we compare the polarizable embedding (PE) and frozen-density embedding (FDE) models. The former is a classical embedding model, whereas the latter is a density-based QM embedding model. Our comparison focuses on solvent effects on optical spectra of solutes. This is a typical scenario where super-system calculations including the solvent environment become prohibitively large. We formulate a common theoretical framework for PE and FDE models and systematically investigate how PE and FDE approximate solvent effects. Generally, differences are found to be small, except in cases where electron spill-out becomes problematic in the classical frameworks. In these cases, however, atomic pseudopotentials can reduce the electron-spill-out issue.

Toward Accurate Prediction of Ion Mobility in Organic Semiconductors by Atomistic Simulation

A multiscale scheme (MLMS: Multi-Level Multi-Scale) to predict the ion mobility (μ) of amorphous organic semiconductors is proposed, which was successfully applied to the hole mobility predictions of 14 organic systems. An inverse relationship between μ and reorganization energy is observed due to local polaronic distortions. Another moderate inverse correlation between μ and distribution of site energy change exists, representing the effects of geometric flexibility. The former and the latter represent the intramolecular and intermolecular geometric effects, respectively. In addition, a linear correlation between transfer coupling and μ is observed, showing the importance of orbital overlaps between monomers. Especially, the highest hole mobility of C6-2TTN is due to its large transfer coupling. On the other hand, another high hole mobility of CBP turned out to come from the high first neighbor density (ρ_FND) of its first self-solvation, emphasizing the proper description of amorphous structural configurations with a sufficiently large number of monomers. In general, systems with either unusually high transfer coupling or high first neighbor density can potentially have high μ regardless of geometric effects. Especially, the newly suggested design parameter, ρ_FND, is pointing to a new direction as opposed to the traditional π-conjugation strategy.

Nanoscale and Real-Time Nuclear-Electronic Dynamics Simulation Study of Charge Transfer at the Donor-Acceptor Interface in Organic Photovoltaics

Charge-transfer (CT) processes in donor-acceptor interfaces of organic photovoltaics have been challenging targets for computational chemistry owing to their nanoscale and ultrafast nature. Herein, we report real-time nuclear-electronic dynamics simulations of CT processes in a nanometer-scale donor-acceptor interface model composed of a donor poly(3-hexylthiophene-2,5-diyl) crystal and an acceptor [6,6]-phenyl-C₆₁-butyric acid methyl ester aggregate. The simulations were realized using our original reduced-scaling computational technique, namely, patchwork-approximation-based Ehrenfest dynamics. The results illustrated the CT pathway with atomic resolution, thereby rationalizing the observed excitation-energy dependence of the quantity of CT. Further, nuclear motion, which is affected by the electronic dynamics, was observed to play a significant role in the CT process by modulating molecular orbital energies. The present study suggests that microscopic CT processes strongly depend on local structures of disordered donor-acceptor interfaces as well as coupling between nuclear and electronic dynamics.

Graph-based quantum response theory and shadow Born-Oppenheimer molecular dynamics

Graph-based linear scaling electronic structure theory for quantum-mechanical molecular dynamics simulations [A. M. N. Niklasson et al., J. Chem. Phys. 144, 234101 (2016)] is adapted to the most recent shadow potential formulations of extended Lagrangian Born-Oppenheimer molecular dynamics, including fractional molecular-orbital occupation numbers [A. M. N. Niklasson, J. Chem. Phys. 152, 104103 (2020) and A. M. N. Niklasson, Eur. Phys. J. B 94, 164 (2021)], which enables stable simulations of sensitive complex chemical systems with unsteady charge solutions. The proposed formulation includes a preconditioned Krylov subspace approximation for the integration of the extended electronic degrees of freedom, which requires quantum response calculations for electronic states with fractional occupation numbers. For the response calculations, we introduce a graph-based canonical quantum perturbation theory that can be performed with the same natural parallelism and linear scaling complexity as the graph-based electronic structure calculations for the unperturbed ground state. The proposed techniques are particularly well-suited for semi-empirical electronic structure theory, and the methods are demonstrated using self-consistent charge density-functional tight-binding theory both for the acceleration of self-consistent field calculations and for quantum-mechanical molecular dynamics simulations. Graph-based techniques combined with the semi-empirical theory enable stable simulations of large, complex chemical systems, including tens-of-thousands of atoms.

Extending multi-layer energy-based fragment method for excited-state calculations of large covalently bonded fragment systems

Recently, we developed a low-scaling Multi-Layer Energy-Based Fragment (MLEBF) method for accurate excited-state calculations and nonadiabatic dynamics simulations of nonbonded fragment systems. In this work, we extend the MLEBF method to treat covalently bonded fragment ones. The main idea is cutting a target system into many fragments according to chemical properties. Fragments with dangling bonds are first saturated by chemical groups; then, saturated fragments, together with the original fragments without dangling bonds, are grouped into different layers. The accurate total energy expression is formulated with the many-body energy expansion theory, in combination with the inclusion-exclusion principle that is used to delete the contribution of chemical groups introduced to saturate dangling bonds. Specifically, in a two-layer MLEBF model, the photochemically active and inert layers are calculated with high-level and efficient electronic structure methods, respectively. Intralayer and interlayer energies can be truncated at the two- or three-body interaction level. Subsequently, through several systems, including neutral and charged covalently bonded fragment systems, we demonstrate that MLEBF can provide accurate ground- and excited-state energies and gradients. Finally, we realize the structure, conical intersection, and path optimizations by combining our MLEBF program with commercial and free packages, e.g., ASE and SciPy. These developments make MLEBF a practical and reliable tool for studying complex photochemical and photophysical processes of large nonbonded and bonded fragment systems.

Divide-and-Conquer Linear-Scaling Quantum Chemical Computations

Fragmentation and embedding schemes are of great importance when applying quantum-chemical calculations to more complex and attractive targets. The divide-and-conquer (DC)-based quantum-chemical model is a fragmentation scheme that can be connected to embedding schemes. This feature article explains several DC-based schemes developed by the authors over the last two decades, which was inspired by the pioneering study of DC self-consistent field (SCF) method by Yang and Lee (J. Chem. Phys. 1995, 103, 5674-5678). First, the theoretical aspects of the DC-based SCF, electron correlation, excited-state, and nuclear orbital methods are described, followed by the two-component relativistic theory, quantum-mechanical molecular dynamics simulation, and the introduction of three programs, including DC-based schemes. Illustrative applications confirmed the accuracy and feasibility of the DC-based schemes.

Structure optimization with stochastic density functional theory

Linear-scaling techniques for Kohn-Sham density functional theory are essential to describe the ground state properties of extended systems. Still, these techniques often rely on the localization of the density matrix or accurate embedding approaches, limiting their applicability. In contrast, stochastic density functional theory (sDFT) achieves linear- and sub-linear scaling by statistically sampling the ground state density without relying on embedding or imposing localization. In return, ground state observables, such as the forces on the nuclei, fluctuate in sDFT, making optimizing the nuclear structure a highly non-trivial problem. In this work, we combine the most recent noise-reduction schemes for sDFT with stochastic optimization algorithms to perform structure optimization within sDFT. We compare the performance of the stochastic gradient descent approach and its variations (stochastic gradient descent with momentum) with stochastic optimization techniques that rely on the Hessian, such as the stochastic Broyden-Fletcher-Goldfarb-Shanno algorithm. We further provide a detailed assessment of the computational efficiency and its dependence on the optimization parameters of each method for determining the ground state structure of bulk silicon with varying supercell dimensions.

Adaptive Subsystem Density Functional Theory

Subsystem density functional theory (DFT) is emerging as a powerful electronic structure method for large-scale simulations of molecular condensed phases and interfaces. Key to its computational efficiency is the use of approximate nonadditive noninteracting kinetic energy functionals. Unfortunately, currently available nonadditive functionals lead to inaccurate results when the subsystems interact strongly such as when they engage in chemical reactions. This work disrupts the status quo by devising a workflow that extends subsystem DFT's applicability also to strongly interacting subsystems. This is achieved by implementing a fully automated adaptive definition of subsystems which is realized during geometry optimizations or ab initio molecular dynamics simulations. The new method prescribes subsystem merging and splitting events redistributing the resources (both for work and data) in an efficient way making use of modern parallelization strategies and object-oriented programming. We showcase the method with examples probing from moderate-to-strong inter-subsystem interactions, opening the door to using subsystem DFT for modeling chemical reactions in molecular condensed phases with a black box computational tool.

Density Embedding Method for Nanoscale Molecule-Metal Interfaces

In this work, we extend the applicability of standard Kohn-Sham DFT (KS-DFT) to model realistically sized molecule-metal interfaces where the metal slabs venture into the tens of nanometers in size. Employing state-of-the-art noninteracting kinetic energy functionals, we describe metallic subsystems with orbital-free DFT and combine their electronic structure with molecular subsystems computed at the KS-DFT level resulting in a multiscale subsystem DFT method. The method reproduces within a few millielectronvolts the binding energy difference of water and carbon dioxide molecules adsorbed on the top and hollow sites of an Al(111) surface compared to KS-DFT of the combined supersystem. It is also robust for Born-Oppenheimer molecular dynamics simulations. Very large system sizes are approached with standard computing resources thanks to a parallelization scheme that avoids accumulation of memory at the gather-scatter stage. The results as presented are encouraging and open the door to ab initio simulations of realistically sized, mesoscopic molecule-metal interfaces.

Body-Ordered Approximations of Atomic Properties

We show that the local density of states (LDOS) of a wide class of tight-binding models has a weak body-order expansion. Specifically, we prove that the resulting body-order expansion for analytic observables such as the electron density or the energy has an exponential rate of convergence both at finite Fermi-temperature as well as for insulators at zero Fermi-temperature. We discuss potential consequences of this observation for modelling the potential energy landscape, as well as for solving the electronic structure problem.

Quantum Algorithm of the Divide-and-Conquer Unitary Coupled Cluster Method with a Variational Quantum Eigensolver

The variational quantum eigensolver (VQE) with shallow or constant-depth quantum circuits is one of the most pursued approaches in the noisy intermediate-scale quantum (NISQ) devices with incoherent errors. In this study, the divide-and-conquer (DC) linear scaling technique, which divides the entire system into several fragments, is applied to the VQE algorithm based on the unitary coupled cluster (UCC) method, denoted as DC-qUCC/VQE, to reduce the number of required qubits. The unitarity of the UCC ansatz that enables the evaluation of the total energy as well as various molecular properties as expectation values can be easily implemented on quantum devices because the quantum gates are unitary operators themselves. Based on this feature, the present DC-qUCC/VQE algorithm is designed to conserve the total number of electrons in the entire system using the density matrix evaluated on a quantum computer. Numerical assessments clarified that the energy errors of the DC-qUCC/VQE calculations decrease by using the constraint of the total number of electrons. Furthermore, the DC-qUCC/VQE algorithm could reduce the number of quantum gates and shows the possibility of decreasing incoherent errors.

Combined fragment-based machine learning force field with classical force field and its application in the NMR calculations of macromolecules in solutions

We have developed a combined fragment-based machine learning (ML) force field and molecular mechanics (MM) force field for simulating the structures of macromolecules in solutions, and then compute its NMR chemical shifts with the generalized energy-based fragmentation (GEBF) approach at the level of density functional theory (DFT). In this work, we first construct Gaussian approximation potential based on GEBF subsystems of macromolecules for MD simulations and then a GEBF-based neural network (GEBF-NN) with deep potential model for the studied macromolecule. Then, we develop a GEBF-NN/MM force field for macromolecules in solutions by combining the GEBF-NN force field for the solute molecule and ff14SB force field for solvent molecules. Using the GEBF-NN/MM MD simulation to generate snapshot structures of solute/solvent clusters, we then perform the NMR calculations with the GEBF approach at the DFT level to calculate NMR chemical shifts of the solute molecule. Taking a heptamer of oligopyridine-dicarboxamides in chloroform solution as an example, our results show that the GEBF-NN force field is quite accurate for this heptamer by comparing with the reference DFT results. For this heptamer in chloroform solution, both the GEBF-NN/MM and classical MD simulations could lead to helical structures from the same initial extended structure. The GEBF-DFT NMR results indicate that the GEBF-NN/MM force field could lead to more accurate NMR chemical shifts on hydrogen atoms by comparing with the experimental NMR results. Therefore, the GEBF-NN/MM force field could be employed for predicting more accurate dynamical behaviors than the classical force field for complex systems in solutions.

Constructing high-accuracy theoretical Raman spectra of SARS-CoV-2 spike proteins based on a large fragment method

In order to control COVID-19, rapid and accurate detection of the pathogenic, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), is an urgent task. The target spike proteins of SARS-CoV-2 have been detected experimentally via Raman spectroscopy. However, there lacks high-accuracy theoretical Raman spectra of the spike proteins to as a standard reference for the clinic diagnostic purpose. In this paper, we propose a large fragment method to construct the high-precision Raman spectra for the spike proteins. The large fragment method not only reduces the calculation error but also improves the accuracy of the protein Raman spectra by completely calculating the interactions within the large fragment. The Pearson correlation coefficient of theoretical Raman spectra is greater than 0.929 or more. Compared with the experimental spectra, the characteristic patterns are easily visible. This work provides a detection standard for the spike proteins which shall bring a step closer to the fast recognition of SARS-CoV-2 via Raman spectroscopy method.

Effective electronic forces and potentials from ab initio path integral Monte Carlo simulations

The rigorous description of correlated quantum many-body systems constitutes one of the most challenging tasks in contemporary physics and related disciplines. In this context, a particularly useful tool is the concept of effective pair potentials that take into account the effects of the complex many-body medium consistently. In this work, we present extensive, highly accurate ab initio path integral Monte Carlo (PIMC) results for the effective interaction and the effective force between two electrons in the presence of the uniform electron gas. This gives us a direct insight into finite-size effects, thereby, opening up the possibility for novel domain decompositions and methodological advances. In addition, we present unassailable numerical proof for an effective attraction between two electrons under moderate coupling conditions, without the mediation of an underlying ionic structure. Finally, we compare our exact PIMC results to effective potentials from linear-response theory, and we demonstrate their usefulness for the description of the dynamic structure factor. All PIMC results are made freely available online and can be used as a thorough benchmark for new developments and approximations.

A double exponential coupled cluster theory in the fragment molecular orbital framework

Fragmentation-based methods enable electronic structure calculations for large chemical systems through partitioning them into smaller fragments. Here, we have developed and benchmarked a dual exponential operator-based coupled cluster theory to account for high-rank electronic correlation of large chemical systems within the fragment molecular orbital (FMO) framework. Upon partitioning the molecular system into several fragments, the zeroth order reference determinants for each fragment and fragment pair are constructed in a self-consistent manner with two-body FMO expansion. The dynamical correlation is induced through a dual exponential ansatz with a set of fragment-specific rank-one and rank-two operators that act on the individual reference determinants. While the single and double excitations for each fragment are included through the conventional rank-one and rank-two cluster operators, the triple excitation space is spanned via the contraction between the cluster operators and a set of rank-two scattering operators over a few optimized fragment-specific occupied and virtual orbitals. Thus, the high-rank dynamical correlation effects within the FMO framework are computed with rank-one and rank-two parametrization of the wave operator, leading to significant reduction in the number of variables and associated computational scaling over the conventional methods. Through a series of pilot numerical applications on various covalent and non-covalently bonded systems, we have shown the quantitative accuracy of the proposed methodology compared to canonical, as well as FMO-based coupled-cluster single double triple. The accuracy of the proposed method is shown to be systematically improvable upon increasing the number of contractible occupied and virtual molecular orbitals employed to simulate triple excitations.

Local electronic structure analysis by ab initio elongation method: A benchmark using DNA block polymers

The ab initio elongation (ELG) method based on a polymerization concept is a feasible way to perform linear-scaling electronic structure calculations for huge aperiodic molecules while maintaining computational accuracy. In the method, the electronic structures are sequentially elongated by repeating (1) the conversion of canonical molecular orbitals (CMOs) to region-localized MOs (RLMOs), that is, active RLMOs localized onto a region close to an attacking monomer or frozen RLMOs localized onto the remaining region, and the subsequent (2) partial self-consistent-field calculations for an interaction space composed of the active RLMOs and the attacking monomer. For each ELG process, one can obtain local CMOs for the interaction space and the corresponding local orbital energies. Local site information, such as the local highest-occupied/lowest-unoccupied MOs, can be acquired with linear-scaling efficiency by correctly including electronic effects from the frozen region. In this study, we performed a local electronic structure analysis using the ELG method for various DNA block polymers with different sequential patterns. This benchmark aimed to confirm the effectiveness of the method toward the efficient detection of a singular local electronic structure in unknown systems as a future practical application. We discussed the high-throughput efficiency of our method and proposed a strategy to detect singular electronic structures by combining with a machine learning technique.

Regularized Localized Molecular Orbitals in a Divide-and-Conquer Approach for Linear Scaling Calculations

Non-orthogonal localized molecular orbitals (NOLMOs) have been employed as building blocks for the divide-and-conquer (DC) linear scaling method. The NOLMOs are calculated from subsystems and used for constructing the density matrix (DM) of the entire system, instead of the subsystem DM in the original DC approach. Also, unlike the original DC method, the inverse electronic temperature parameter β is not needed anymore. Furthermore, a new regularized localization approach for NOLMOs has been developed, in which the localization cost function is a sum of the spatial spread function, as in the Boys method, and the kinetic energy, as a regularization measure to limit the oscillation of the NOLMOs. The optimal weight of the kinetic energy can be determined by optimization with analytical gradients. The resulting regularized NOLMOs have enhanced smoothness and better transferability because of reduced kinetic energies. Compared with the original DC, while NOLMO-DC has a similar computational linear scaling cost, the accuracy of NOLMO-DC is better by several orders of magnitude for large conjugated systems and by about 1 order of magnitude for other systems. The NOLMO-DC method is thus a promising development of the DC approach for linear scaling calculations.

Fragmentation Method for Computing Quantum Mechanics and Molecular Mechanics Gradients for Force Matching: Validation with Hydration Free Energy Predictions Using Adaptive Force Matching

A fragmentation approach referred to as a simple overlapping region method for force matching (SORForM) is presented. SORForM is designed to enable efficient computation of quantum mechanical (QM) forces for large molecules and is validated in the framework of adaptive force matching (AFM) to develop solute models in water. The SORForM method divides a molecule into overlapping QM regions with each region containing a gradient zone and a buffer zone. The buffer zone ensures that the atoms in the gradient zone have their surroundings unchanged with fragmentation. The performance of the method is validated with mefenamic acid and linalyl acetate by comparing the hydration free energies of AFM models developed with and without SORForM. The AFM hydration free energies are also compared with that of the experiments. The models developed with B3LYP-D3(BJ) and def2-TZVP are in excellent agreement with experiments. Our work shows that PBE-D3(BJ) provides less satisfactory results when compared to B3LYP-D3(BJ). The def2-TZVP basis set is found to greatly improve the agreement with experiments when compared to a double-zeta quality basis set.

Study of crystal properties based on attention mechanism and crystal graph convolutional neural network

The prediction of crystal properties has always been limited by huge computational costs. In recent years, the rise of machine learning methods has gradually made it possible to study crystal properties on a large scale. We propose an attention mechanism-based crystal graph convolutional neural network, which builds a machine learning model by inputting crystallographic information files and target properties. In our research, the attention mechanism is introduced in the crystal graph convolutional neural network (CGCNN) to learn the local chemical environment, and node normalization is added to reduce the risk of overfitting. We collect structural information and calculation data of about 36 000 crystals and examine the prediction performance of the models for the formation energy, total energy, bandgap, and Fermi energy of crystals in our research. Compared with the CGCNN, it is found that the accuracy (ACCU) of the predicted properties can be further improved to varying degrees by the introduction of the attention mechanism. Moreover, the total magnetization and bandgap can be classified under the same neural network framework. The classification ACCU of wide bandgap semiconductor crystals with a bandgap threshold of 2.3 eV reaches 93.2%, and the classification ACCU of crystals with a total magnetization threshold of 0.5 μ_Breaches 88.8%. The work is helpful to realize large-scale prediction and classification of crystal properties, accelerating the discovery of new functional crystal materials.

Forces from Stochastic Density Functional Theory under Nonorthogonal Atom-Centered Basis Sets

We develop a formalism for calculating forces on the nuclei within the linear-scaling stochastic density functional theory (sDFT) in a nonorthogonal atom-centered basis set representation (Fabian et al. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2019, 9, e1412, 10.1002/wcms.1412) and apply it to the Tryptophan Zipper 2 (Trp-zip2) peptide solvated in water. We use an embedded-fragment approach to reduce the statistical errors (fluctuation and systematic bias), where the entire peptide is the main fragment and the remaining 425 water molecules are grouped into small fragments. We analyze the magnitude of the statistical errors in the forces and find that the systematic bias is of the order of 0.065 eV/Å (∼1.2 × 10-3Eh/a0) when 120 stochastic orbitals are used, independently of system size. This magnitude of bias is sufficiently small to ensure that the bond lengths estimated by stochastic DFT (within a Langevin molecular dynamics simulation) will deviate by less than 1% from those predicted by a deterministic calculation.

All-Atom Nonadiabatic Dynamics Simulation of Hybrid Graphene Nanoribbons Based on Wannier Analysis and Machine Learning

Trajectory surface hopping combined with ab initio electronic structure calculations is a popular and powerful approach for on-the-fly nonadiabatic dynamics simulations. For large systems, however, this remains a significant challenge because of the unaffordable computational cost of large-scale electronic structure calculations. Here, we present an efficient divide-and-conquer approach to construct the system Hamiltonian based on Wannier analysis and machine learning. In detail, the large system under investigation is first decomposed into small building blocks, and then all possible segments formed by building blocks within a cutoff distance are found out. Ab initio molecular dynamics is carried out to generate a sequence of geometries for each equivalent segment with periodicity. The Hamiltonian matrices in the maximum localized Wannier function (MLWF) basis are obtained for all geometries and utilized to train artificial neural networks (ANNs) for the structure-dependent Hamiltonian elements. Taking advantage of the orthogonality and spatial locality of MLWFs, the one-electron Hamiltonian of a large system at arbitrary geometry can be directly constructed by the trained ANNs. As demonstrations, we study charge transport in a zigzag graphene nanoribbon (GNR), a coved GNR, and a series of hybrid GNRs with a state-of-the-art surface hopping method. The interplay between delocalized and localized states is found to determine the electron dynamics in hybrid GNRs. Our approach has successfully studied GNRs with >10 000 atoms, paving the way for efficient and reliable all-atom nonadiabatic dynamics simulation of general systems.

Scalable Ehrenfest Molecular Dynamics Exploiting the Locality of Density-Functional Tight-Binding Hamiltonian

To explore the science behind excited-state dynamics in high-complexity chemical systems, a scalable nonadiabatic molecular dynamics (MD) technique is indispensable. In this study, by treating the electronic degrees of freedom at the density-functional tight-binding level, we developed and implemented a reduced scaling and multinode-parallelizable Ehrenfest MD method. To achieve this goal, we introduced a concept called patchwork approximation (PA), where the effective Hamiltonian for real-time propagation of the electronic density matrix is partitioned into a set of local parts. Numerical results for giant icosahedral fullerenes, which comprise up to 6000 atoms, suggest that the scaling of the present PA-based method is less than quadratic, which yields a significant advantage over the conventional cubic scaling method in terms of computational time. The acceleration by the parallelization on multiple nodes was also assessed. Furthermore, the electronic and structural dynamics resulting from the perturbation by the external electric field were accurately reproduced with the PA, even when the electronic excitation was spatially delocalized.

Computational and data driven molecular material design assisted by low scaling quantum mechanics calculations and machine learning

Electronic structure methods based on quantum mechanics (QM) are widely employed in the computational predictions of the molecular properties and optoelectronic properties of molecular materials. The computational costs of these QM methods, ranging from density functional theory (DFT) or time-dependent DFT (TDDFT) to wave-function theory (WFT), usually increase sharply with the system size, causing the curse of dimensionality and hindering the QM calculations for large sized systems such as long polymer oligomers and complex molecular aggregates. In such cases, in recent years low scaling QM methods and machine learning (ML) techniques have been adopted to reduce the computational costs and thus assist computational and data driven molecular material design. In this review, we illustrated low scaling ground-state and excited-state QM approaches and their applications to long oligomers, self-assembled supramolecular complexes, stimuli-responsive materials, mechanically interlocked molecules, and excited state processes in molecular aggregates. Variable electrostatic parameters were also introduced in the modified force fields with the polarization model. On the basis of QM computational or experimental datasets, several ML algorithms, including explainable models, deep learning, and on-line learning methods, have been employed to predict the molecular energies, forces, electronic structure properties, and optical or electrical properties of materials. It can be conceived that low scaling algorithms with periodic boundary conditions are expected to be further applicable to functional materials, perhaps in combination with machine learning to fast predict the lattice energy, crystal structures, and spectroscopic properties of periodic functional materials.

Graph-Theory-Based Molecular Fragmentation for Efficient and Accurate Potential Surface Calculations in Multiple Dimensions

We present a multitopology molecular fragmentation approach, based on graph theory, to calculate multidimensional potential energy surfaces in agreement with post-Hartree-Fock levels of theory but at the density functional theory cost. A molecular assembly is coarse-grained into a set of graph-theoretic nodes that are then connected with edges to represent a collection of locally interacting subsystems up to an arbitrary order. Each of the subsystems is treated at two levels of electronic structure theory, the result being used to construct many-body expansions that are embedded within an ONIOM scheme. These expansions converge rapidly with the many-body order (or graphical rank) of subsystems and capture many-body interactions accurately and efficiently. However, multiple graphs, and hence multiple fragmentation topologies, may be defined in molecular configuration space that may arise during conformational sampling or from reactive, bond breaking and bond formation, events. Obtaining the resultant potential surfaces is an exponential scaling proposition, given the number of electronic structure computations needed. We utilize a family of graph-theoretic representations within a variational scheme to obtain multidimensional potential surfaces at a reduced cost. The fast convergence of the graph-theoretic expansion with increasing order of many-body interactions alleviates the exponential scaling cost for computing potential surfaces, with the need to only use molecular fragments that contain a fewer number of quantum nuclear degrees of freedom compared to the full system. This is because the dimensionality of the conformational space sampled by the fragment subsystems is much smaller than the full molecular configurational space. Additionally, we also introduce a multidimensional clustering algorithm, based on physically defined criteria, to reduce the number of energy calculations by orders of magnitude. The molecular systems benchmarked include coupled proton motion in protonated water wires. The potential energy surfaces and multidimensional nuclear eigenstates obtained are shown to be in very good agreement with those from explicit post-Hartree-Fock calculations that become prohibitive as the number of quantum nuclear dimensions grows. The developments here provide a rigorous and efficient alternative to this important chemical physics problem.

Protein Response Effects on Cofactor Excitation Energies from First Principles: Augmenting Subsystem Time-Dependent Density-Functional Theory with Many-Body Expansion Techniques

We investigate the possibility of describing protein response effects on a chromophore excitation by means of subsystem time-dependent density-functional theory (sTDDFT) in combination with a many-body expansion (MBE) approach. While sTDDFT is in principle intrinsically able to include such contributions, addressing cofactor excitations in protein models or entire proteins with full environment-response treatments is currently out of reach. Taking different model structures of the green fluorescent protein (GFP) and bovine rhodopsin as examples, we demonstrate that an embedded-MBE approach based on sTDDFT in its simplest version leads to a good agreement of the predicted protein response effect already at second order. To reproduce reference response effects from nonsubsystem TDDFT calculations quantitatively (error ≤ 5%), however, a third- or even fourth-order MBE may be required. For the latter case, we explore a selective inclusion of fourth-order terms that drastically reduces the computational burden. In addition, we demonstrate how this sTDDFT-MBE treatment can be utilized as an analysis tool to identify residues with dominant response contributions. This, in turn, can be employed to arrive at smaller structural models for light-absorbing proteins, which still feature all of the main characteristics in terms of photoresponse properties.

Electronic Couplings and Electrostatic Interactions Behind the Light Absorption of Retinal Proteins

The photo-functional chromophore retinal exhibits a wide variety of optical absorption properties depending on its intermolecular interactions with surrounding proteins and other chromophores. By utilizing these properties, microbial and animal rhodopsins express biological functions such as ion-transport and signal transduction. In this review, we present the molecular mechanisms underlying light absorption in rhodopsins, as revealed by quantum chemical calculations. Here, symmetry-adapted cluster-configuration interaction (SAC-CI), combined quantum mechanical and molecular mechanical (QM/MM), and transition-density-fragment interaction (TDFI) methods are used to describe the electronic structure of the retinal, the surrounding protein environment, and the electronic coupling between chromophores, respectively. These computational approaches provide successful reproductions of experimentally observed absorption and circular dichroism (CD) spectra, as well as insights into the mechanisms of unique optical properties in terms of chromophore-protein electrostatic interactions and chromophore-chromophore electronic couplings. On the basis of the molecular mechanisms revealed in these studies, we also discuss strategies for artificial design of the optical absorption properties of rhodopsins.

Generalized energy-based fragmentation approach for calculations of solvation energies of large systems

A generalized energy-based fragmentation (GEBF) approach has been combined with a universal solvation model based on solute electron density (SMD) to compute the solvation energies of general large systems (such as protein molecules) in solutions. In the GEBF-SMD method, the solvation energy of a target system could be combined by the corresponding solvation energies of various subsystems, each of which is embedded in the background point charges and surface charges on the surface of solute cavity at the positions of its atoms and neighbouring atoms outside of the subsystem. Our results show that the GEBF-SMD model could reproduce the conventional SMD solvation energies quite well for various proteins in solutions, and could significantly reduce the computational costs for the SMD calculations of large proteins. In addition, the GEBF-SMD approach is almost independent of the basis sets and the types of solvents (including protic, polar, and nonpolar ones). Also, the GEBF-SMD approach could reproduce the relative energies of various conformers of large systems in solutions. Therefore, the GEBF-SMD method is expected to be applicable for computing the solvation energies of a broad range of large systems.

Ab initio multi-level layered elongation method and its application to local interaction analysis between DNA bulge and ligand molecules

A multi-level layered elongation method was developed for efficiently analyzing the electronic states of local structures in large bio/nano-systems at the full ab initio level of theory. The original elongation method developed during the last three decades in our group has focused on the system in one direction from one terminal to the other terminal to sequentially construct the electronic states of a polymer, called a theoretical synthesis of polymers. In this study, an important region termed the central (C) part is targeted in a large polymer and the remainder are terminal (T) parts. The electronic structures along with polymer elongation are calculated repeatedly from both end T parts to the C central part at the same time. The important C part is treated with large basis sets (high level) and the other regions are treated with small basis sets (low level) in the ab initio theoretical framework. The electronic structures besides the C part can be reused for other systems with different structures at the C part, which renders the method computationally efficient. This multi-level layered elongation method was applied to the investigation on DNA single bulge recognition of small molecules (ligands). The reliability and validity of our approach were examined in comparison with the results obtained by direct calculations using a conventional quantum chemical method for the entire system. Furthermore, stabilization energies by the formation of the complex of bulge DNA and a ligand were estimated with basis set superposition error corrections incorporated into the elongation method.

The Middle Science: Traversing Scale In Complex Many-Body Systems

A roadmap is developed that integrates simulation methodology and data science methods to target new theories that traverse the multiple length- and time-scale features of many-body phenomena.

iOI: An Iterative Orbital Interaction Approach for Solving the Self-Consistent Field Problem

An iterative orbital interaction (iOI) approach is proposed to solve, in a bottom-up fashion, the self-consistent field problem in quantum chemistry. While it belongs grossly to the family of fragment-based quantum chemical methods, iOI is distinctive in that (1) it divides and conquers not only the energy but also the wave function and that (2) the subsystem sizes are automatically determined by successively merging neighboring small subsystems until they are just enough for converging the wave function to a given accuracy. Orthonormal occupied and virtual localized molecular orbitals are obtained in a natural manner, which can be used for all post-SCF purposes.

Physics-Inspired Structural Representations for Molecules and Materials

The first step in the construction of a regression model or a data-driven analysis, aiming to predict or elucidate the relationship between the atomic-scale structure of matter and its properties, involves transforming the Cartesian coordinates of the atoms into a suitable representation. The development of atomic-scale representations has played, and continues to play, a central role in the success of machine-learning methods for chemistry and materials science. This review summarizes the current understanding of the nature and characteristics of the most commonly used structural and chemical descriptions of atomistic structures, highlighting the deep underlying connections between different frameworks and the ideas that lead to computationally efficient and universally applicable models. It emphasizes the link between properties, structures, their physical chemistry, and their mathematical description, provides examples of recent applications to a diverse set of chemical and materials science problems, and outlines the open questions and the most promising research directions in the field.

Three-Layer Multiscale Approach Based on Extremely Localized Molecular Orbitals to Investigate Enzyme Reactions

Quantum mechanics/molecular mechanics (QM/MM) calculations are widely used embedding techniques to computationally investigate enzyme reactions. In most QM/MM computations, the quantum mechanical region is treated through density functional theory (DFT), which offers the best compromise between chemical accuracy and computational cost. Nevertheless, to obtain more accurate results, one should resort to wave function-based methods, which however lead to a much larger computational cost already for relatively small QM subsystems. To overcome this drawback, we propose the coupling of our QM/ELMO (quantum mechanics/extremely localized molecular orbital) approach with molecular mechanics, thus introducing the three-layer QM/ELMO/MM technique. The QM/ELMO strategy is an embedding method in which the chemically relevant part of the system is treated at the quantum mechanical level, while the rest is described through frozen ELMOs. Since the QM/ELMO method reproduces results of fully QM computations within chemical accuracy and with a much lower computational effort, it can be considered a suitable strategy to extend the range of applicability and accuracy of the QM/MM scheme. In this paper, other than briefly presenting the theoretical bases of the QM/ELMO/MM technique, we will also discuss its validation on the well-tested deprotonation of acetyl coenzyme A by aspartate in citrate synthase.