Protein-ligand interactions from a quantum fragmentation perspective: The case of the SARS-CoV-2 main protease interacting with α-ketoamide inhibitors.

We present a hybrid, multi-method, computational scheme for protein/ligand systems well suited to be used on modern and forthcoming massively parallel computing systems. The scheme relies on a multi-scale polarizable molecular modeling, approach to perform molecular dynamics simulations, and on an efficient Density Functional Theory (DFT) linear scaling method to post-process simulation snapshots. We use this scheme to investigate recent α-ketoamide inhibitors targeting the main protease of the SARS-CoV-2 virus. We assessed the reliability and the coherence of the hybrid scheme, in particular, by checking the ability of MM and DFT to reproduce results from high-end ab initio computations regarding such inhibitors. The DFT approach enables an a posteriori fragmentation of the system and an investigation into the strength of interaction among identified fragment pairs. We show the necessity of accounting for a large set of plausible protease/inhibitor conformations to generate reliable interaction data. Finally, we point out ways to further improve α-ketoamide inhibitors to more strongly interact with particular protease domains neighboring the active site.

Complexity reduction in density functional theory: Locality in space and energy.

We present recent developments of the NTChem program for performing large scale hybrid density functional theory calculations on the supercomputer Fugaku. We combine these developments with our recently proposed complexity reduction framework to assess the impact of basis set and functional choice on its measures of fragment quality and interaction. We further exploit the all electron representation to study system fragmentation in various energy envelopes. Building off this analysis, we propose two algorithms for computing the orbital energies of the Kohn-Sham Hamiltonian. We demonstrate that these algorithms can efficiently be applied to systems composed of thousands of atoms and as an analysis tool that reveals the origin of spectral properties.

Probing disorder in 2CzPN using core and valence states.

Molecules which exhibit thermally activated delayed fluorescence (TADF) show great promise for use in efficient, environmentally-friendly OLEDs, and thus the design of new TADF emitters is an active area of research. However, when used in devices, they are typically in the form of disordered thin films, where both the external molecular environment and thermally-induced internal variations in parameters such as the torsion angle can strongly influence their electronic structure. In this work, we use density functional theory and X-ray photoelectron spectroscopy to investigate the impact of disorder on both core and valence states in the TADF emitter 2CzPN (1,2-bis(carbazol-9-yl)-4,5-dicyanobenzene). By simulating gas phase molecules displaying varying levels of disorder, we assess the relative sensitivity of the different states to factors such as varying torsion angle. The theoretical results for both core and valence states show good agreement with experiment, thereby also highlighting the advantages of our approach for interpreting experimental spectra of large aromatic molecules, which are too complex to interpret based solely on experimental data.

Flexibilities of wavelets as a computational basis set for large-scale electronic structure calculations.

The BigDFT project was started in 2005 with the aim of testing the advantages of using a Daubechies wavelet basis set for Kohn-Sham (KS) density functional theory (DFT) with pseudopotentials. This project led to the creation of the BigDFT code, which employs a computational approach with optimal features of flexibility, performance, and precision of the results. In particular, the employed formalism has enabled the implementation of an algorithm able to tackle DFT calculations of large systems, up to many thousands of atoms, with a computational effort that scales linearly with the number of atoms. In this work, we recall some of the features that have been made possible by the peculiar properties of Daubechies wavelets. In particular, we focus our attention on the usage of DFT for large-scale systems. We show how the localized description of the KS problem, emerging from the features of the basis set, is helpful in providing a simplified description of large-scale electronic structure calculations. We provide some examples on how such a simplified description can be employed, and we consider, among the case-studies, the SARS-CoV-2 main protease.

The CECAM electronic structure library and the modular software development paradigm.

First-principles electronic structure calculations are now accessible to a very large community of users across many disciplines, thanks to many successful software packages, some of which are described in this special issue. The traditional coding paradigm for such packages is monolithic, i.e., regardless of how modular its internal structure may be, the code is built independently from others, essentially from the compiler up, possibly with the exception of linear-algebra and message-passing libraries. This model has endured and been quite successful for decades. The successful evolution of the electronic structure methodology itself, however, has resulted in an increasing complexity and an ever longer list of features expected within all software packages, which implies a growing amount of replication between different packages, not only in the initial coding but, more importantly, every time a code needs to be re-engineered to adapt to the evolution of computer hardware architecture. The Electronic Structure Library (ESL) was initiated by CECAM (the European Centre for Atomic and Molecular Calculations) to catalyze a paradigm shift away from the monolithic model and promote modularization, with the ambition to extract common tasks from electronic structure codes and redesign them as open-source libraries available to everybody. Such libraries include "heavy-duty" ones that have the potential for a high degree of parallelization and adaptation to novel hardware within them, thereby separating the sophisticated computer science aspects of performance optimization and re-engineering from the computational science done by, e.g., physicists and chemists when implementing new ideas. We envisage that this modular paradigm will improve overall coding efficiency and enable specialists (whether they be computer scientists or computational scientists) to use their skills more effectively and will lead to a more dynamic evolution of software in the community as well as lower barriers to entry for new developers. The model comes with new challenges, though. The building and compilation of a code based on many interdependent libraries (and their versions) is a much more complex task than that of a code delivered in a single self-contained package. Here, we describe the state of the ESL, the different libraries it now contains, the short- and mid-term plans for further libraries, and the way the new challenges are faced. The ESL is a community initiative into which several pre-existing codes and their developers have contributed with their software and efforts, from which several codes are already benefiting, and which remains open to the community.

Multipole-preserving quadratures for the discretization of functions in real-space electronic structure calculations.

Discretizing an analytic function on a uniform real-space grid is often done via a straightforward collocation method. This is ubiquitous in all areas of computational physics and quantum chemistry. An example in density functional theory (DFT) is given by the external potential or the pseudo-potential describing the interaction between ions and electrons. The accuracy of the collocation method used is therefore very important for the reliability of subsequent treatments like self-consistent field solutions of the electronic structure problems. By construction, the collocation method introduces numerical artifacts typical of real-space treatments, like the so-called egg-box error, which may spoil the numerical stability of the description when the real-space grid is too coarse. As the external potential is an input of the problem, even a highly precise computational treatment cannot cope this inconvenience. We present in this paper a new quadrature scheme that is able to exactly preserve the moments of a given analytic function even for large grid spacings, while reconciling with the traditional collocation method when the grid spacing is small enough. In the context of real-space electronic structure calculations, we show that this method improves considerably the stability of the results for large grid spacings, opening up the path towards reliable low-accuracy DFT calculations with a reduced number of degrees of freedom.

Accurate and efficient linear scaling DFT calculations with universal applicability.

Density functional theory calculations are computationally extremely expensive for systems containing many atoms due to their intrinsic cubic scaling. This fact has led to the development of so-called linear scaling algorithms during the last few decades. In this way it becomes possible to perform ab initio calculations for several tens of thousands of atoms within reasonable walltimes. However, even though the use of linear scaling algorithms is physically well justified, their implementation often introduces some small errors. Consequently most implementations offering such a linear complexity either yield only a limited accuracy or, if one wants to go beyond this restriction, require a tedious fine tuning of many parameters. In our linear scaling approach within the BigDFT package, we were able to overcome this restriction. Using an ansatz based on localized support functions expressed in an underlying Daubechies wavelet basis - which offers ideal properties for accurate linear scaling calculations - we obtain an amazingly high accuracy and a universal applicability while still keeping the possibility of simulating large system with linear scaling walltimes requiring only a moderate demand of computing resources. We prove the effectiveness of our method on a wide variety of systems with different boundary conditions, for single-point calculations as well as for geometry optimizations and molecular dynamics.

Fragment approach to constrained density functional theory calculations using Daubechies wavelets.

In a recent paper, we presented a linear scaling Kohn-Sham density functional theory (DFT) code based on Daubechies wavelets, where a minimal set of localized support functions are optimized in situ and therefore adapted to the chemical properties of the molecular system. Thanks to the systematically controllable accuracy of the underlying basis set, this approach is able to provide an optimal contracted basis for a given system: accuracies for ground state energies and atomic forces are of the same quality as an uncontracted, cubic scaling approach. This basis set offers, by construction, a natural subset where the density matrix of the system can be projected. In this paper, we demonstrate the flexibility of this minimal basis formalism in providing a basis set that can be reused as-is, i.e., without reoptimization, for charge-constrained DFT calculations within a fragment approach. Support functions, represented in the underlying wavelet grid, of the template fragments are roto-translated with high numerical precision to the required positions and used as projectors for the charge weight function. We demonstrate the interest of this approach to express highly precise and efficient calculations for preparing diabatic states and for the computational setup of systems in complex environments.

Daubechies wavelets for linear scaling density functional theory.

We demonstrate that Daubechies wavelets can be used to construct a minimal set of optimized localized adaptively contracted basis functions in which the Kohn-Sham orbitals can be represented with an arbitrarily high, controllable precision. Ground state energies and the forces acting on the ions can be calculated in this basis with the same accuracy as if they were calculated directly in a Daubechies wavelets basis, provided that the amplitude of these adaptively contracted basis functions is sufficiently small on the surface of the localization region, which is guaranteed by the optimization procedure described in this work. This approach reduces the computational costs of density functional theory calculations, and can be combined with sparse matrix algebra to obtain linear scaling with respect to the number of electrons in the system. Calculations on systems of 10,000 atoms or more thus become feasible in a systematic basis set with moderate computational resources. Further computational savings can be achieved by exploiting the similarity of the adaptively contracted basis functions for closely related environments, e.g., in geometry optimizations or combined calculations of neutral and charged systems.

Predicting potential SARS-CoV-2 mutations of concern via full quantum mechanical modelling.

Ab initio quantum mechanical models can characterize and predict intermolecular binding, but only recently have models including more than a few hundred atoms gained traction. Here, we simulate the electronic structure for approximately 13 000 atoms to predict and characterize binding of SARS-CoV-2 spike variants to the human ACE2 (hACE2) receptor using the quantum mechanics complexity reduction (QM-CR) approach. We compare four spike variants in our analysis: Wuhan, Omicron, and two Omicron-based variants. To assess binding, we mechanistically characterize the energetic contribution of each amino acid involved, and predict the effect of select single amino acid mutations. We validate our computational predictions experimentally by comparing the efficacy of spike variants binding to cells expressing hACE2. At the time we performed our simulations (December 2021), the mutation A484K which our model predicted to be highly beneficial to ACE2 binding had not been identified in epidemiological surveys; only recently (August 2023) has it appeared in variant BA.2.86. We argue that our computational model, QM-CR, can identify mutations critical for intermolecular interactions and inform the engineering of high-specificity interactors.

Selecting boron fullerenes by cage-doping mechanisms.

So far, no boron fullerenes were synthesized: more compact sp(3)-bonded clusters are energetically preferred. To circumvent this, metallic clusters have been suggested by Pochet et al. [Phys. Rev. B 83, 081403(R) (2011)] as "seeds" for a possible synthesis which would topologically protect the sp(2) sector of the configuration space. In this paper, we identify a basic pentagonal unit which allows a balance between the release of strain and the self-doping rule. We formulate a guiding principle for the stability of boron fullerenes, which takes the form of an isolated filled pentagon rule (IFPR). The role of metallic clusters is then reexamined. It is shown that the interplay of the IFPR and the seed-induced doping breaks polymorphism and its related problems: it can effectively select between different isomers and reduce the reactivity of the boron shells. The balance between self and exterior doping represents the best strategy for boron buckyball synthesis.

Accurate complex scaling of three dimensional numerical potentials.

The complex scaling method, which consists in continuing spatial coordinates into the complex plane, is a well-established method that allows to compute resonant eigenfunctions of the time-independent Schrödinger operator. Whenever it is desirable to apply the complex scaling to investigate resonances in physical systems defined on numerical discrete grids, the most direct approach relies on the application of a similarity transformation to the original, unscaled Hamiltonian. We show that such an approach can be conveniently implemented in the Daubechies wavelet basis set, featuring a very promising level of generality, high accuracy, and no need for artificial convergence parameters. Complex scaling of three dimensional numerical potentials can be efficiently and accurately performed. By carrying out an illustrative resonant state computation in the case of a one-dimensional model potential, we then show that our wavelet-based approach may disclose new exciting opportunities in the field of computational non-Hermitian quantum mechanics.

Norm-conserving pseudopotentials with chemical accuracy compared to all-electron calculations.

By adding a nonlinear core correction to the well established dual space Gaussian type pseudopotentials for the chemical elements up to the third period, we construct improved pseudopotentials for the Perdew-Burke-Ernzerhof [J. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996)] functional and demonstrate that they exhibit excellent accuracy. Our benchmarks for the G2-1 test set show average atomization energy errors of only half a kcal/mol. The pseudopotentials also remain highly reliable for high pressure phases of crystalline solids. When supplemented by empirical dispersion corrections [S. Grimme, J. Comput. Chem. 27, 1787 (2006); S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010)] the average error in the interaction energy between molecules is also about half a kcal/mol. The accuracy that can be obtained by these pseudopotentials in combination with a systematic basis set is well superior to the accuracy that can be obtained by commonly used medium size Gaussian basis sets in all-electron calculations.

Biochemical and Phylogenetic Analysis of Italian Phaseolus vulgaris Cultivars as Sources of α-Amylase and α-Glucosidase Inhibitors.

Phaseolus vulgaris α-amylase inhibitor (α-AI) is a protein that has recently gained commercial interest, as it inhibits mammalian α-amylase activity, reducing the absorption of dietary carbohydrates. Numerous studies have reported the efficacy of preparations based on this protein on the control of glycaemic peaks in type-2 diabetes patients and in overweight subjects. A positive influence on microbiota regulation has also been described. In this work, ten insufficiently studied Italian P. vulgaris cultivars were screened for α-amylase- and α-glucosidase-inhibiting activity, as well as for the absence of antinutritional compounds, such as phytohemagglutinin (PHA). All the cultivars presented α-glucosidase-inhibitor activity, while α-AI was missing in two of them. Only the Nieddone cultivar (ACC177) had no haemagglutination activity. In addition, the partial nucleotide sequence of the α-AI gene was identified with the degenerate hybrid oligonucleotide primer (CODEHOP) strategy to identify genetic variability, possibly linked to functional α-AI differences, expression of the α-AI gene, and phylogenetic relationships. Molecular studies showed that α-AI was expressed in all the cultivars, and a close similarity between the Pisu Grogu and Fasolu cultivars' α-AI and α-AI-4 isoform emerged from the comparison of the partially reconstructed primary structures. Moreover, mechanistic models revealed the interaction network that connects α-AI with the α-amylase enzyme characterized by two interaction hotspots (Asp38 and Tyr186), providing some insights for the analysis of the α-AI primary structure from the different cultivars, particularly regarding the structure-activity relationship. This study can broaden the knowledge about this class of proteins, fuelling the valorisation of Italian agronomic biodiversity through the development of commercial preparations from legume cultivars.

Experimental-theoretical study of laccase as a detoxifier of aflatoxins.

We investigate laccase-mediated detoxification of aflatoxins, fungal carcinogenic food contaminants. Our experimental comparison between two aflatoxins with similar structures (AFB1 and AFG2) shows significant differences in laccase-mediated detoxification. A multi-scale modeling approach (Docking, Molecular Dynamics, and Density Functional Theory) identifies the highly substrate-specific changes required to improve laccase detoxifying performance. We employ a large-scale density functional theory-based approach, involving more than 7000 atoms, to identify the amino acid residues that determine the affinity of laccase for aflatoxins. From this study we conclude: (1) AFB1 is more challenging to degrade, to the point of complete degradation stalling; (2) AFG2 is easier to degrade by laccase due to its lack of side products and favorable binding dynamics; and (3) ample opportunities to optimize laccase for aflatoxin degradation exist, especially via mutations leading to π-π stacking. This study identifies a way to optimize laccase for aflatoxin bioremediation and, more generally, contributes to the research efforts aimed at rational enzyme optimization.

Efficient and accurate solver of the three-dimensional screened and unscreened Poisson's equation with generic boundary conditions.

We present an explicit solver of the three-dimensional screened and unscreened Poisson's equation, which combines accuracy, computational efficiency, and versatility. The solver, based on a mixed plane-wave/interpolating scaling function representation, can deal with any kind of periodicity (along one, two, or three spatial axes) as well as with fully isolated boundary conditions. It can seamlessly accommodate a finite screening length, non-orthorhombic lattices, and charged systems. This approach is particularly advantageous because convergence is attained by simply refining the real space grid, namely without any adjustable parameter. At the same time, the numerical method features O(NlogN) scaling of the computational cost (N being the number of grid points) very much like plane-wave methods. The methodology, validated on model systems, is tailored for leading-edge computer simulations of materials (including ab initio electronic structure computations), but it might as well be beneficial for other research domains.

Transition-Based Constrained DFT for the Robust and Reliable Treatment of Excitations in Supramolecular Systems.

Despite the variety of available computational approaches, state-of-the-art methods for calculating excitation energies, such as time-dependent density functional theory (TDDFT), are computationally demanding and thus limited to moderate system sizes. Here, we introduce a new variation of constrained DFT (CDFT), wherein the constraint corresponds to a particular transition (T), or a combination of transitions, between occupied and virtual orbitals, rather than a region of the simulation space as in traditional CDFT. We compare T-CDFT with TDDFT and ΔSCF results for the low-lying excited states (S1 and T1) of a set of gas-phase acene molecules and OLED emitters and with reference results from the literature. At the PBE level of theory, T-CDFT outperforms ΔSCF for both classes of molecules, while also proving to be more robust. For the local excitations seen in the acenes, T-CDFT and TDDFT perform equally well. For the charge transfer (CT)-like excitations seen in the OLED molecules, T-CDFT also performs well, in contrast to the severe energy underestimation seen with TDDFT. In other words, T-CDFT is equally applicable to both local excitations and CT states, providing more reliable excitation energies at a much lower computational cost than TDDFT cost. T-CDFT is designed for large systems and has been implemented in the linear-scaling BigDFT code. It is therefore ideally suited for exploring the effects of explicit environments on excitation energies, paving the way for future simulations of excited states in complex realistic morphologies, such as those which occur in OLED materials.

Probing the mutational landscape of the SARS-CoV-2 spike protein via quantum mechanical modeling of crystallographic structures.

We employ a recently developed complexity-reduction quantum mechanical (QM-CR) approach, based on complexity reduction of density functional theory calculations, to characterize the interactions of the SARS-CoV-2 spike receptor binding domain (RBD) with ACE2 host receptors and antibodies. QM-CR operates via ab initio identification of individual amino acid residue's contributions to chemical binding and leads to the identification of the impact of point mutations. Here, we especially focus on the E484K mutation of the viral spike protein. We find that spike residue 484 hinders the spike's binding to the human ACE2 receptor (hACE2). In contrast, the same residue is beneficial in binding to the bat receptor Rhinolophus macrotis ACE2 (macACE2). In agreement with empirical evidence, QM-CR shows that the E484K mutation allows the spike to evade categories of neutralizing antibodies like C121 and C144. The simulation also shows how the Delta variant spike binds more strongly to hACE2 compared to the original Wuhan strain, and predicts that a E484K mutation can further improve its binding. Broad agreement between the QM-CR predictions and experimental evidence supports the notion that ab initio modeling has now reached the maturity required to handle large intermolecular interactions central to biological processes.

Energy landscape of fullerene materials: a comparison of boron to boron nitride and carbon.

Using the minima hopping global geometry optimization method on the density functional potential energy surface we show that the energy landscape of boron clusters is glasslike. Larger boron clusters have many structures which are lower in energy than the cages. This is in contrast to carbon and boron nitride systems which can be clearly identified as structure seekers. The differences in the potential energy landscape explain why carbon and boron nitride systems are found in nature whereas pure boron fullerenes have not been found. We thus present a methodology which can make predictions on the feasibility of the synthesis of new nanostructures.

Use of the dual potential to rationalize the occurrence of some DNA lesions (pyrimidic dimers).

Exploiting the locality of the chemical potential of an excited state when it is evaluated using the ground state Density Functional Theory (DFT), a new local descriptor for excited states has been proposed (J. Chem. Theory Comput.2009, 5, 2274). This index is based on the assumption that the relaxation of the electronic density toward that of the ground state drives the chemical reactivity of excited states. The sign of the descriptor characterizes the electrophilic or nucleophilic behavior of atomic regions. Through an exact excited state DFT formalism provided by Gross, Oliveira, and Kohn, a mathematical argument is given for this descriptor only for the first excited state. It is afterward used to rationalize the occurrence and the regioselectivity of some DNA lesions based on the [2 + 2] cycloaddition between two adjacent bases.