CoChat
Sign Up
PreprintCited by 1Advanced Chemical Physics StudiesAdvanced Physical and Chemical Molecular InteractionsProtein Structure and Dynamics

The DOI matches its Crossref record; 2 of the 2 references it cites were found in a registry or an index.

Bridging Quantum Mechanics and Biology at the Million-Atom Scale

Open access
Garcia, Martin, Wieners, Luc·2025·Communications Chemistry
DOI 10.21203/rs.3.rs-7327472/v1PMC PMC12458572

Preprint — published as 10.1038/s42004-026-02038-y

Open this paper in CoChat — free accountOpen PDF

Abstract

Quantum-mechanical simulations provide the most fundamental description of matter, yet their computational cost has limited applications to systems containing at most thousands of atoms. Here, we present an all-electron quantum-mechanical framework that extends the accessible system size to the multimillion-atom regime. By combining an algorithmically optimized Hartree-Fock formalism with divide-and-conquer, our new approach efficiently handles structures containing millions of atoms with subatomic resolution. We demonstrate this approach on very large biological systems, including a complete bacteriophage in water, totalling over 150 million electrons, representing, to our knowledge, the largest quantum-mechanical calculation performed to date. Our framework allows computing spectral data for DNA and drugs and enables protein structure assessments in strong agreement with structure evaluations by AlphaFold. This approach opens new avenues in quantum physics, structural biology, spectroscopy, bioinformatics, medicine, and materials science.

Overview

The claims and statistics of this paper haven’t been checked yet.

Run the check — free, no account needed. It reads the open-access text and checks the sources it cites.

A few free checks a day · papers already verified open instantly · results are public

Integrity signals

Registry

Crossref · 10.21203/rs.3.rs-7327472/v1

Cites

2 of 2 references verified0 citation errors · 0 not checkable

Statistics

Figures

Figure 5: Cryo-EM and computed densities. Comparison of electronic densities determined via cryogenic electron microscopy (a,c,e) with quantum-mechanically computed electronic densities (b,d,f) for a flagellar motor (PDB 8WL2)33, a vault protein complex (PDB 7PKR)34 and a bacteriophage (PDB 6Q3G)1. a,c,e were visualised using the Mol* Viewer85 from rcsb.org42.
Figure 5: Cryo-EM and computed densities. Comparison of electronic densities determined via cryogenic electron microscopy (a,c,e) with quantum-mechanically computed electronic densities (b,d,f) for a flagellar motor (PDB 8WL2)33, a vault protein complex (PDB 7PKR)34 and a bacteriophage (PDB 6Q3G)1. a,c,e were visualised using the Mol* Viewer85 from rcsb.org42.CC BY · PubMed Central
Figure 6: Calculated electronic density of the vault protein complex. Detailed depiction of the electronic density of the vault protein complex PDB 7PKR34, computed using the Hartree-Fock method presented in this work. a, Atomistic representation of the complete molecule. b, Projection of the electronic density of one half on the vault structure along the x-axis. c, Close-up view of the region highlighted in b. The protein backbone is displayed as coloured lines, with brown indicating regions of low and blue of high electronic density. For clarity, the electronic density of the solvent is omitted and only contributions from valence basis functions are shown.
Figure 6: Calculated electronic density of the vault protein complex. Detailed depiction of the electronic density of the vault protein complex PDB 7PKR34, computed using the Hartree-Fock method presented in this work. a, Atomistic representation of the complete molecule. b, Projection of the electronic density of one half on the vault structure along the x-axis. c, Close-up view of the region highlighted in b. The protein backbone is displayed as coloured lines, with brown indicating regions of low and blue of high electronic density. For clarity, the electronic density of the solvent is omitted and only contributions from valence basis functions are shown.CC BY · PubMed Central
Figure 7: Comparison of quantum energies and pLDDT scores. Graphical abstract illustrating one of the main findings of this work: the correlation between AlphaFold’s confidence scores and quantum-mechanical atomic energies. A structure prediction from AlphaFold (based on a known amino acid sequence) is analysed using Hartree–Fock calculations to determine atomic energies. These energies are then compared with the pLDDT scores, which are a direct output from AlphaFold. The visualisation of the protein structure PDB AF_AFA0A023FF81F1 was generated using ChimeraX29 with data from rcsb.org42.
Figure 7: Comparison of quantum energies and pLDDT scores. Graphical abstract illustrating one of the main findings of this work: the correlation between AlphaFold’s confidence scores and quantum-mechanical atomic energies. A structure prediction from AlphaFold (based on a known amino acid sequence) is analysed using Hartree–Fock calculations to determine atomic energies. These energies are then compared with the pLDDT scores, which are a direct output from AlphaFold. The visualisation of the protein structure PDB AF_AFA0A023FF81F1 was generated using ChimeraX29 with data from rcsb.org42.CC BY · PubMed Central
Figure 8: Graphical abstract for large-scale quantum-mechanical calculations. Accessing the electronic structure (right) of large-scale biological structures as the bacteriophage P681 (left) is made possible by dividing the structure into clusters and employing an improved algorithm for electronic integral computation.
Figure 8: Graphical abstract for large-scale quantum-mechanical calculations. Accessing the electronic structure (right) of large-scale biological structures as the bacteriophage P681 (left) is made possible by dividing the structure into clusters and employing an improved algorithm for electronic integral computation.CC BY · PubMed Central
Figure 9: Atomic energies versus neighbouring atom count. For the predicted protein structure PDB AF_AFA0A023IKK2F142 (a), the structure evaluation method using atomic energies from Hartree–Fock is compared with an approach that uses the number of neighbouring atoms as a metric. Higher local atom density is typically associated with more accurate predictions, resulting in a loose correlation between the neighbouring atom count (b) and AlphaFold’s pLDDT scores (d). However, the evaluation of protein structures with Hartree–Fock atomic energies (c) provides a more accurate assessment. For example, the alpha helix at the bottom left of the structure is correctly assigned high confidence by the quantum approach but not by the neighbour count method. Similarly, the linker between the two top-left domains is identified as stable by the Hartree–Fock energies but misclassified as a weak region by the atom density method.
Figure 9: Atomic energies versus neighbouring atom count. For the predicted protein structure PDB AF_AFA0A023IKK2F142 (a), the structure evaluation method using atomic energies from Hartree–Fock is compared with an approach that uses the number of neighbouring atoms as a metric. Higher local atom density is typically associated with more accurate predictions, resulting in a loose correlation between the neighbouring atom count (b) and AlphaFold’s pLDDT scores (d). However, the evaluation of protein structures with Hartree–Fock atomic energies (c) provides a more accurate assessment. For example, the alpha helix at the bottom left of the structure is correctly assigned high confidence by the quantum approach but not by the neighbour count method. Similarly, the linker between the two top-left domains is identified as stable by the Hartree–Fock energies but misclassified as a weak region by the atom density method.CC BY · PubMed Central
Figure 10: Computed electronic density of the bacteriophage P68. Detailed visualisation of the electronic density of the bacteriophage P681, calculated using the presented Hartree-Fock method. Colours range from black/blue (low density) to white (high density). For improved clarity, the electronic density of the solvent is omitted, and only contributions from valence basis functions are shown.
Figure 10: Computed electronic density of the bacteriophage P68. Detailed visualisation of the electronic density of the bacteriophage P681, calculated using the presented Hartree-Fock method. Colours range from black/blue (low density) to white (high density). For improved clarity, the electronic density of the solvent is omitted, and only contributions from valence basis functions are shown.CC BY · PubMed Central
Figure 1: Divide-and-conquer clustering for linear-scaling quantum calculations. a, The protein complex (PDB 2M3X)28 is shown as an example system. b, The molecule is spatially partitioned into 12 clusters using the k-means algorithm. c, For larger systems in a solution, such as a flagellar motor in water, the individual clusters are obtained by subdividing the system using a three-dimensional grid. d, Corresponding grid-based clusters are shown. Visualisation of the protein structure in a was done with ChimeraX29.
Figure 1: Divide-and-conquer clustering for linear-scaling quantum calculations. a, The protein complex (PDB 2M3X)28 is shown as an example system. b, The molecule is spatially partitioned into 12 clusters using the k-means algorithm. c, For larger systems in a solution, such as a flagellar motor in water, the individual clusters are obtained by subdividing the system using a three-dimensional grid. d, Corresponding grid-based clusters are shown. Visualisation of the protein structure in a was done with ChimeraX29.CC BY · PubMed Central
Figure 2: Visualisation of the modified all-electron Hartree–Fock—exemplified on the bacteriophage P681 in a solution of water. a,b, Atomic-scale resolved structure of the bacteriophage itself and in solution. c,d, Calculated electronic density as a projection along the x-axis and as a three-dimensional picture (see also Figure 10). e,f, Close-up of the projected electronic density of c after successive enlargements. Brown corresponds to low and blue to high electronic density values. For a better visualisation the electronic density of the solvent is not displayed and only valence basis function contributions are shown. Visualisation of the protein structure in a was done with ChimeraX29.
Figure 2: Visualisation of the modified all-electron Hartree–Fock—exemplified on the bacteriophage P681 in a solution of water. a,b, Atomic-scale resolved structure of the bacteriophage itself and in solution. c,d, Calculated electronic density as a projection along the x-axis and as a three-dimensional picture (see also Figure 10). e,f, Close-up of the projected electronic density of c after successive enlargements. Brown corresponds to low and blue to high electronic density values. For a better visualisation the electronic density of the solvent is not displayed and only valence basis function contributions are shown. Visualisation of the protein structure in a was done with ChimeraX29.CC BY · PubMed Central
Figure 3: Computed spectra of biomolecules in comparison with experiments. Absorption spectra were calculated for Actinomycin5 (a), DNA with 8 and 21 base pairs3,4 (c) and Actinomycin bound to DNA6 (e) and compared to experimental data from Miao et al.37 (b), Ozdemir et al.38 (d) and Toader et al.39 (f). Visualisation of the biomolecular structures in a, c and e was done with ChimeraX29.
Figure 3: Computed spectra of biomolecules in comparison with experiments. Absorption spectra were calculated for Actinomycin5 (a), DNA with 8 and 21 base pairs3,4 (c) and Actinomycin bound to DNA6 (e) and compared to experimental data from Miao et al.37 (b), Ozdemir et al.38 (d) and Toader et al.39 (f). Visualisation of the biomolecular structures in a, c and e was done with ChimeraX29.CC BY · PubMed Central
Figure 4: Comparison of quantum-mechanical atomic energies with AlphaFold confidence scores. Atomic energies from Hartree–Fock calculations are compared with AlphaFold’s pLDDT score for three predicted protein structures: Evasin P1126 (PDB AF_AFA0A023FF81F1) (a), MHC Class II Beta Chain (AF_AFA0A023IKK2F1) (b) and Friend of Echinoid, Isoform H (AF_AFA0A023GPK8F1) (c). For each protein, the left column shows ribbon diagrams and non-hydrogen atoms. The middle column (b,e,h) displays atomic energies computed in a solvated environment (10 Å water shell), smoothened using a Savitzky–Golay filter, scaled by element number, and linearly rescaled to match the pLDDT range. The right column (c,f,i) shows the corresponding AlphaFold pLDDT scores mapped onto atoms. Energies of outlier atoms were clipped to match the [0,100] range of pLDDT values. Protein structures in a, d and g were visualized using ChimeraX19. Data from rcsb.org42.
Figure 4: Comparison of quantum-mechanical atomic energies with AlphaFold confidence scores. Atomic energies from Hartree–Fock calculations are compared with AlphaFold’s pLDDT score for three predicted protein structures: Evasin P1126 (PDB AF_AFA0A023FF81F1) (a), MHC Class II Beta Chain (AF_AFA0A023IKK2F1) (b) and Friend of Echinoid, Isoform H (AF_AFA0A023GPK8F1) (c). For each protein, the left column shows ribbon diagrams and non-hydrogen atoms. The middle column (b,e,h) displays atomic energies computed in a solvated environment (10 Å water shell), smoothened using a Savitzky–Golay filter, scaled by element number, and linearly rescaled to match the pLDDT range. The right column (c,f,i) shows the corresponding AlphaFold pLDDT scores mapped onto atoms. Energies of outlier atoms were clipped to match the [0,100] range of pLDDT values. Protein structures in a, d and g were visualized using ChimeraX19. Data from rcsb.org42.CC BY · PubMed Central

Verification

Claims the paper makes

The claims and statistics of this paper haven’t been checked yet.

Run the check — free, no account needed. It reads the open-access text and checks the sources it cites.

A few free checks a day · papers already verified open instantly · results are public

Self-check

Percentages, sums and reported differences re-computed from the paper’s own numbers.

The self-check needs the paper’s full text.

Citations

Cites

Show all 2 references

Research works printed but not linked: 83

  1. 1. Hrebík, D., et al. (2019). Structure and genome ejection mechanism of Staphylococcus aureus phage P68. Science advances, 5(10). · DOI
  2. 3. Thota, N. A. R. A. S. A. I. A. H., Li, X. H., Bingman, C., & Sundaralingam, M. (1993). High-resolution refinement of the hexagonal A-DNA octamer d (GTGTACAC) at 1.4 Å. Biological Crystallography, 49(2), 282–291. · DOI
  3. 4. Masliah, G., René, B., Zargarian, L., Fermandjian, S., & Mauffret, O. (2008). Identification of intrinsic dynamics in a DNA sequence preferentially cleaved by topoisomerase II enzyme. Journal of molecular biology, 381(3), 692–706. · DOI
  4. 5. Schäfer, M., Sheldrick, G. M., Bahner, I., & Lackner, H. (1998). Crystal structures of actinomycin D and actinomycin Z3. Angewandte Chemie International Edition, 37(17), 2381–2384. · DOI
  5. 6. Hou, M. H., Robinson, H., Gao, Y. G., & Wang, A. H. J. (2002). Crystal structure of actinomycin D bound to the CTG triplet repeat sequences linked to neurological diseases. Nucleic acids research, 30(22), 4910–4917. · DOI
  6. 7. Van Der Spoel, D., et al. (2005). GROMACS: fast, flexible, and free. Journal of computational chemistry, 26(16), 1701–1718. · DOI
  7. 8. Case, D. A., et al. (2005). The Amber biomolecular simulation programs. Journal of computational chemistry, 26(16), 1668–1688. · DOI · Paper page
  8. 9. Vanommeslaeghe, K., et al. (2010). CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all‐atom additive biological force fields. Journal of computational chemistry, 31(4), 671–690. · DOI
  9. 10. Senn, H. M., & Thiel, W. (2009). QM/MM methods for biomolecular systems. Angewandte Chemie International Edition, 48(7), 1198–1229. · DOI
  10. 11. Dauber-Osguthorpe, P., & Hagler, A. T. (2019). Biomolecular force fields: where have we been, where are we now, where do we need to go and how do we get there? Journal of computer-aided molecular design, 33(2), 133–203. · DOI
  11. 12. Soler, J. M., et al. (2002). The SIESTA method for ab initio order-N materials simulation. Journal of Physics: Condensed Matter, 14(11), 2745. · DOI
  12. 13. Neese, F. (2012). The ORCA program system. Wiley Interdisciplinary Reviews: Computational Molecular Science, 2(1), 73–78. · DOI
  13. 14. Stocks, R., et al. (2024, November). Breaking the million-electron and 1 EFLOP/s barriers: biomolecular-Scale ab initio molecular dynamics using MP2 potentials. In SC24: International Conference for High Performance Computing, Networking, Storage and Analysis (pp. 1–12). IEEE. · DOI
  14. 15. Andermatt, S., Cha, J., Schiffmann, F., & VandeVondele, J. (2016). Combining linear-scaling DFT with subsystem DFT in Born–Oppenheimer and Ehrenfest molecular dynamics simulations: From molecules to a virus in solution. Journal of chemical theory and computation, 12(7), 3214–3227. · DOI
  15. 16. Mohr, S., et al. (2015). Accurate and efficient linear scaling DFT calculations with universal applicability. Physical Chemistry Chemical Physics, 17(47), 31360–31370. · DOI
  16. 17. de Luzuriaga, I. O., et al. (2022). Semi-empirical and linear-scaling DFT methods to characterize duplex DNA and G-quadruplexes in the presence of interacting small molecules. Physical Chemistry Chemical Physics, 24(19), 11510–11519. · DOI
  17. 18. Gundelach, L., Fox, T., Tautermann, C. S., & Skylaris, C. K. (2021). Protein-ligand free energies of binding from full-protein DFT calculations: convergence and choice of exchange–correlation functional. Physical Chemistry Chemical Physics, 23(15), 9381–9393. · DOI
  18. 19. He, X., & Merz Jr, K. M. (2010). Divide and conquer Hartree–Fock calculations on proteins. Journal of chemical theory and computation, 6(2), 405–411. · DOI
  19. 20. Guidon, M., Hutter, J., & VandeVondele, J. (2009). Robust periodic Hartree–Fock exchange for large-scale simulations using Gaussian basis sets. Journal of chemical theory and computation, 5(11), 3010–3021. · DOI
  20. 21. Cole, D. J., & Hine, N. D. (2016). Applications of large-scale density functional theory in biology. Journal of Physics: Condensed Matter, 28(39), 393001. · DOI
  21. 22. Shimojo, F., Kalia, R. K., Nakano, A., & Vashishta, P. (2008). Divide-and-conquer density functional theory on hierarchical real-space grids: Parallel implementation and applications. Physical Review B, 77(8), 085103. · DOI
  22. 23. Ordejón, P., Artacho, E., & Soler, J. M. (1996). Self-consistent order-N density-functional calculations for very large systems. Physical review B, 53(16), R10441. · DOI
  23. 24. Fonseca Guerra, C., Snijders, J. G., Te Velde, G. T., & Baerends, E. J. (1998). Towards an order-N DFT method. Theoretical Chemistry Accounts, 99, 391–403. · DOI
  24. 25. Goletto, L., Kjønstad, E. F., Folkestad, S. D., Høyvik, I. M., & Koch, H. (2021). Linear-scaling implementation of multilevel hartree–fock theory. Journal of Chemical Theory and Computation, 17(12), 7416–7427. · DOI
  25. 26. Nakata, A., et al. (2020). Large scale and linear scaling DFT with the CONQUEST code. The Journal of chemical physics, 152(16). · DOI
  26. 27. Bowler, D. R., & Miyazaki, T. (2010). Calculations for millions of atoms with density functional theory: linear scaling shows its potential. Journal of Physics: Condensed Matter, 22(7), 074207. · DOI
  27. 28. Varnay, I., Truffault, V., Kessler, H., Coles, M., Solution structure of Ph1500: a homohexameric protein centered on a 12-bladed beta-propeller (2013) https://doi.org/10.2210/pdb2M3X/pdb. · DOI
  28. 29. Meng, E. C., et al. (2023). UCSF ChimeraX: Tools for structure building and analysis. Protein Science, 32(11), e4792. · DOI
  29. 30. Strout, D. L., & Scuseria, G. E. (1995). A quantitative study of the scaling properties of the Hartree–Fock method. The Journal of chemical physics, 102(21), 8448–8452. · DOI
  30. 31. Royant, A., et al. (2001). X-ray structure of sensory rhodopsin II at 2.1-Å resolution. Proceedings of the National Academy of Sciences, 98(18), 10131–10136. · DOI
  31. 32. Whitten, J. L. (1966). Gaussian Lobe Function Expansions of Hartree—Fock Solutions for the First-Row Atoms and Ethylene. The Journal of Chemical Physics, 44(1), 359–364. · DOI
  32. 33. Tan, J., et al. (2024). Structural basis of the bacterial flagellar motor rotational switching. Cell Research, 34(11), 788–801. · DOI
  33. 34. Guerra, P., et al. (2022). Symmetry disruption commits vault particles to disassembly. Science Advances, 8(6). · DOI
  34. 35. Müller, C., Sharma, M., & Sierka, M. (2020). Real-time time‐dependent density functional theory using density fitting and the continuous fast multipole method. Journal of Computational Chemistry, 41(30), 2573–2582. · DOI
  35. 36. Liu, Y., et al. (2021). Ultra-Fast Computation of Excited-States Spectra for Large Systems: Ultraviolet and Fluorescence Spectra of Proteins. Interdisciplinary Sciences: Computational Life Sciences, 13, 140–146. · DOI
  36. 37. Miao, L., et al. (2017). The anti-quorum sensing activity and bioactive substance of a marine derived Streptomyces. Biotechnology & Biotechnological Equipment, 31(5), 1007–1015. · DOI
  37. 38. Ozdemir, A., Sarioglu, O. F., & Tekinay, T. (2015). Spectroscopic evaluation of DNA–borate interactions. Biological trace element research, 168, 508–515. · DOI
  38. 39. Toader, A. M., Dascalu, I., & Enache, M. (2022). Effect of SDS Micelles on Actinomycin D-DNA Complexes. Acta Chimica Slovenica, 69(2). · DOI
  39. 40. David, A., Islam, S., Tankhilevich, E., & Sternberg, M. J. (2022). The AlphaFold database of protein structures: a biologist’s guide. Journal of molecular biology, 434(2), 167336. · DOI · Paper page
  40. 41. Dill, K. A., Ozkan, S. B., Shell, M. S., & Weikl, T. R. (2008). The protein folding problem. Annual Review of Biophysics, 37, 289–316 · DOI · Paper page
  41. 42. Berman, H. M., et al. (2000). The protein data bank. Nucleic acids research, 28(1), 235–242. · DOI · Paper page
  42. 44. Koch, D., et al. (2024). The analysis of electron densities: From basics to emergent applications. Chemical Reviews, 124(22), 12661–12737. · DOI
  43. 45. Felder, C. E., Prilusky, J., Silman, I., & Sussman, J. L. (2007). A server and database for dipole moments of proteins. Nucleic acids research, 35, 512–521. · DOI · Paper page
  44. 46. Kulik, M., & Dominiak, P. M. (2022). Electron density is not spherical: the many applications of the transferable aspherical atom model. Computational and Structural Biotechnology Journal, 20, 6237–6243. · DOI
  45. 47. Rackers, J. A., Tecot, L., Geiger, M., & Smidt, T. E. (2023). A recipe for cracking the quantum scaling limit with machine learned electron densities. Machine Learning: Science and Technology, 4(1), 015027. Methods Modifications to the Hartree–Fock algorithm Unlike density functional theory (DFT), Hartree–Fock is a wavefunction-based method and does not directly operate on the electronic density. Each electronic interaction integral contains products of four wave functions. As the system grows, the amount of involved wave functions increases, and the number of electronic repulsion integrals (ERIs) scales with the fourth power of the system size, making their evaluation one of the main computational bottlenecks of Hartree–Fock calculations. To address this, we introduce systematic approximations that significantly reduce the computational effort. Both wavefunction-wavefunction30 and density-density interactions are truncated using physically motivated thresholds. Since wave functions decay exponentially with distance, only nearby orbitals contribute significantly to overlaps (see subsection “Density cut-offs”). We assess the relevance of a density by calculating the overlap of its two constituent wave functions. For orbital types with negative signs, we use the absolute of the product to capture significant interactions that may otherwise cancel numerically but still produce meaningful contributions, especially to energy gradients, which are important for force calculations. Once densities below the threshold are excluded, the scaling of ERI evaluations drops from quartic to quadratic. To further reduce complexity to linear scaling (order N), we introduce a cut-off for the Coulomb interaction based on inter-density distance (see subsection “Coulomb cut-offs”). Since electrostatic interactions decay with the inverse of the distance between two charges, we use a Coulomb cut-off between 8 and 12 Å, similar to that employed in classical forces fields for biophysical simulations. The distance of two densities is defined as the distance between the centres of the two wavefunction pairs that generate them. Only ERIs with Coulomb distances below the Coulomb cut-off are calculated. For ERIs with distances near the cut-off, we apply a smoothening function to ensure continuity. An ERI is evaluated if the overlap threshold of both densities is exceeded and their centres lie within the Coulomb cut-off. Additionally, we compute the product of the two density relevance values (see the next subsections), divide it by the distance between their centres, and compare it to a final relevance threshold. This avoids evaluating interactions between weakly overlapping densities and/or those separated by large distances. To further accelerate the process, we adopt an efficient and parallelizable ERI algorithm48,49,50 based on the so-called Gaussian lobe functions32,51 (see subsection “Basis function fitting”). Unlike conventional approaches52,53,54,55, which use angular prefactors (x, y, z) to construct directional orbitals, our method builds positive and negative lobes from distinct Gaussian components. This removes angular dependencies and simplifies the integral expressions, making them uniform across orbital types. As a result, all ERI calculations can be performed by using the same expression, which enables a simple parallelisation for both ERIs and their derivatives (necessary to calculate atomic forces). We also precompute components of the ERIs which only depend on pairs of wave functions. This reduces the full calculation of the electronic repulsion between two densities consisting of in total of four Gaussian functions to only 6 additions, 8 multiplications, 10 memory accesses, 1 square root and the evaluation of the error function (see below). For most interactions with large separations, the error function can be approximated as unity, further reducing computational cost. The following subsections detail the generation of the lobe-based basis functions, the calculation of the density relevance, the calculation of the Coulomb interaction relevance and the mathematical expressions for the modified ERI algorithm. Basis function fitting In electronic structure calculations, wave functions are commonly modelled as Slater-type orbitals (STOs), which for an s-orbital have the form. To enable efficient integration, STOs are typically approximated by using Gaussian functions of the form. Each STO is approximated by a linear combination of Gaussians as This approximation enables an analytic evaluation of key integrals52,53,54,55, relevant for the overlap matrix, the kinetic and nuclei energy as well as the Coulomb and exchange interaction of the electrons in the system. Orbitals with higher angular momenta are typically expressed via Gaussians containing cartesian prefactors (). Although these integrals can also be evaluated analytically, their complexity increases due to the recursive nature of the associated expressions52,53,54,55. Furthermore, different angular momentum combinations require distinct equations, reducing parallelisation efficiency. To avoid this bottleneck, we here adopt Gaussian lobe function expansions32,51, which allow the construction of basis sets without cartesian prefactors. We show in this work that this approach, though currently not widely used, is particularly well-suited to modern computing architectures optimized for high parallelisation. It also enables the evaluation of ERIs using a unified matrix formalism, which is potentially compatible with GPU acceleration56. The algorithms for ERI evaluation are described in the section “Electron repulsion integral algorithm”. For fitting, we start with a Slater-type orbital of the STO-6G type taken from the Basis Set Exchange57. In the case of p-orbitals, we approximate them with 6 Gaussian functions, three paced symmetrically slightly to the left and three to the right of the orbital centre. This corresponds in accuracy to a conventional STO-3G orbital. The fitting is performed on a three-dimensional real space grid. Fit parameters include 3 contraction exponents and 3 coefficients, with the latter multiplied by on one side to model the nodal character of p-orbitals. For the optimisation of the parameters during the fitting process, the trust region algorithm is used. This is done via the function curve_fit from the mathematic Python library scipy58. Density cut-offs A density is defined as the product of two wave functions: To determine whether a given density is relevant for further calculations, we compute a scalar value given by the integral over the absolute value of all pairwise products of Gaussian functions: This equation is valid for s-type orbitals, where the integrand remains positive. Orbitals with nonzero angular momentum (e.g. p-orbitals) include prefactors which would require using the absolute values to ensure that the integrand is always positive. However, using these absolute values makes an analytic evaluation of the overlap integrals impossible. To resolve this, we developed an approach to represent each lobe of orbitals with a separate set of Gaussians. The previously described lobe-based basis functions (which achieve a nodal plane by summing Gaussians of different signs) are not suitable here, as their absolute-value integrals do not correspond to the absolute value of the full orbital. Instead, we independently fit each lobe of a p-type orbital using Gaussian functions, and then reconstruct the full orbital by combining lobes with different signs and different centres. The representation of a p-type orbital used for the density relevance calculations is given by: where is the sign of lobe and each cartesian direction has an individual contraction factor (e.g. ) and offset (e.g. ) for higher flexibility upon fitting. The fitting process is identical to that described in the section “Basis function fitting”. Coulomb cut-offs In the Hartree–Fock formalism, the Coulomb interaction is modelled as the repulsion between two electronic densities. Following the method described above, Coulomb interactions are neglected if the separation between densities exceeds a defined cut-off. However, since densities are not point-like, a suitable distance metric must be established to estimate their effective separation. For this purpose, each density is assigned an origin, calculated as the arithmetic mean of the atomic positions associated with the two basis functions forming that density. The interaction distance between two densities and is defined as: where are the positions of the atoms involved in the four basis functions. A lower and upper cut-off threshold and (typically 8 and 10 Å), are introduced to control the inclusion of interactions. For density pairs with distances between these two thresholds, a smooth transition is applied using a cut-off function to avoid discontinuities. The scaling function (cut-off function) used is which smoothly interpolates between 1 and 0 as increases in the interval The distance interval is linearly mapped to the domain of. Electron repulsion integral algorithm An electron repulsion integral (ERI) describes the interaction between two electronic densities: where each orbital is expanded in terms of Gaussian functions. Each ERI of the individual Gaussian functions is given by: with the confluent hypergeometric function, two overlaps and, a factor and a distance which are defined as: We define and evaluate the confluent hypergeometric function using the following equation: where the error function can be efficiently approximated using: (with coefficients, ,, , and )59. For large the error function rapidly approaches 1, since and, and can be approximated accordingly. We recall that we evaluate the ERIs of the form It will only be necessary to compute a small part of all possible combinations, since we use both a density cut-off and a Coulomb cut-off. Therefore, and to reduce redundant computation, the quadruple summation over all Gaussians is reformulated as a double summation over precomputed Gaussian pairs and, for which the quantities and as well as and, respectively, are stored in lookup tables. This is practical because number of density pairs is considerably smaller than the total number of ERIs. The complete ERI evaluation algorithm can be visualized in terms of a three-dimensional tensor: the first and second axis are the sums described above and the third axis is the list of all relevant combinations. Summing up this tensor along the first two axes yields the full list of ERIs. It must be noted that each integral in this tensor is computed by using the same equation. However, the lengths of the first two axes and are not always the same as the number of used Gaussian functions is not the same for all wave functions. To account for this, we modify the equation for by setting the limits of the sums to and which leads to: where is if and —otherwise it is. Therefore, describes if the index combination is defined for the set of Gaussian functions within. This method allows the evaluation of all with the same equation and is therefore highly parallelizable. Additionally, it saves memory, since no intermediate values are stored (all Gaussian ERIs become contracted immediately in each step of the loop over. Precomputing and storing intermediate values dependent on only one density (,, and ) reduce the computational effort as well. Validation of cut-offs and basis To ensure that the proposed cut-off methods and basis function manipulations do not significantly compromise the accuracy of the calculations, results were benchmarked against the quantum-chemistry package ORCA13,60,61,62,63,64. Specifically, we assessed the agreement between our results and a standard calculation using the STO-3G basis set, focusing on both total and orbital energies. This allowed us to quantify the impact of lobe function expansion and interaction truncations on the computed energies. As a test system, we selected the molecule beta-carotene (C40H56), which involves 541,089,856 unique ERIs, of which 30,572,290 were considered to have relevant density contributions. Among these, 3,946,243 integrals met both density and Coulomb interaction relevance criteria. The total electronic energy calculated using our implementation was − 1528.547 Hartree (Ha), compared with − 1528.663 Ha obtained via ORCA. This corresponds to an absolute deviation of 0.115 Ha, or 0.0075%. For the orbital energies, we compared all valence orbitals and the ten lowest virtual orbitals, yielding a root mean square deviation of 7.98 mHa. Using only the Coulomb cut-off leads to a value of 7.57 mHa, whereas only density cut-offs result in an error of 6.07 mHa. No cut-offs at all result in a deviation of 3.23 mHa which is attributed to the modified basis functions in the lobe function expansion which do not match the STO-3G basis exactly and therefore lead to a small error. In terms of computational performance, ORCA required 192 seconds for a beta-carotene calculation on an Intel i5 processor. In comparison, our implementation with all cut-offs was completed in 21 seconds. When disabling the Coulomb cut-off and using a minimal density threshold of 10− 10, typical of many Hartree–Fock codes, the runtime increased to 131 seconds. The accuracy of subsystem splitting was tested using the protein Evasin P1126 (PDB ID: AF_AFA0A023FF81F1)42. For atomic energy calculations, only valence electrons were considered. Subsystems were constructed by segmenting the amino acid chain of the whole protein into smaller peptide fragments. Comparing a full Hartree-Fock calculation with the subsystem-based approach, we obtained a root mean square deviation of 0.385 Ha for the atomic energies, corresponding to a relative error of 0.444%. Parallel Hartree–Fock on large systems The density cut-off was set to 10− 4 and the Coulomb cut-off to 10 Å with a smooth transition to zero starting at 8 Å. The self-consistent cycle of the Hartree–Fock calculation used a convergence threshold of 10− 6 for the RMSD of the density matrices. The STO-3G basis set was split into individual Gaussian functions to allow an ERI evaluation according to the above-mentioned Gaussian lobe function expansion algorithm. System preparation involved solvation in water using the solvation tool of the programme Visual Molecular Dynamics (VMD)65. The solvated system was then partitioned into subsystems on a 3D mesh. Hartree-Fock calculations were carried out independently for each subsystem. To reconstruct the total electronic density, the density matrix elements were chosen from the subsystem containing both atoms corresponding to the two basis functions. If the basis functions were associated with different subsystems, the mean of the relevant density matrix elements was used. The total electronic density was then calculated as where is the total number of basis functions, is a relevance determining whether the pair basis functions and contributes significantly, and is the corresponding element of the density matrix. Spectral calculation parameters For the quantum spectral calculations, we employed the real-time time-dependent Hartree–Fock (RT-TDHF) formalism, as it scales more favourably for larger systems than the linear-response time-dependent Hartree–Fock (LR-TDHF) approach, which theoretically scales with the sixth power of the system size. In RT-TDHF the system is excited by an electric field pulse and propagated in time using quantum-mechanical time evolution. During the propagation, the dipole moments are recorded at each time step, and the frequency-dependent polarizability and absorption spectra are obtained via discrete Fourier transform of the dipole time series. The systems were propagated over 2,000 steps with a time step of 0.25 atomic units (a.u.), corresponding to a total propagation time of 500 a.u. (12.1 fs). The applied electric field pulse had a strength of 10− 5 a.u., a standard deviation of 0.2 a.u., and a starting time of 2 a.u. During the Fourier transform, an attenuation factor of 0.01 a.u. was used to supress numerical noise. Following established practice, the resulting spectra were rescaled66,67 to account for the systematic shifts in excitation energies inherent to time-dependent Hartree–Fock and density functional theory calculations. A rescaling factor of 1.335 was applied, consistent with the benchmark study by Jacquemin et al., who computed TDHF spectra for a wide range of organic dyes67. Atomic energy processing The atomic energies used for the comparison with AlphaFold’s pLDDT scores were smoothed using a Savitzky-Golay filter68 with a window length of 150. To match the scale of pLDDT values, the computed energies were multiplied by the factor of 15 and shifted by a structure-specific offset (–430, − 440 and − 460 for PDB entries AF_AFA0A023FF81F142, PDB AF_AFA0A023IKK2F142 and AF_AFA0A023GPK8F142, respectively). The rescaled energies were truncated at the corresponding structure’s maximum and minimum pLDDT values to eliminate extreme atomic energy outliers. A modified expression for atomic energies was employed, defined as where is the atomic energy of atom, is the atomic charge (used for normalisation), is the number valence basis functions for atom (excluding core contributions), and is the total number of basis functions., and denote the density matrix, core Hamiltonian, and two-electron interaction matrix elements, respectively. System preparation was based on the aforementioned PDB entries42, with solvation performed using a 10 Å water shell generated via the solvation tool of Chimera69. Programme implementation The described programme is written in the Python programming language70. To achieve high computational performance, all computational demanding functions are written either in matrix form or compiled to machine code with just-in-time (JIT) compilation. For JIT compilation, the library Numba71 was used, and matrix operations were performed with Numpy72 or PyTorch73. Visualisation of results and plotting were performed using the library Matplotlib74. Biomolecular structures and visualisation All calculations in this article were based on biomolecular structures obtained from the RCSB Protein Data Bank (PDB)42. The structures and their full references (including original research articles and PDB entries) are as follows: PDB 2M3X28 (12-bladed propeller protein), PDB 1H6831,75 (rhodopsin), PDB 6Q3G1,76 (bacteriophage P68), PDB 1A7Y5,77 (Actinomycin D), PDB 1D793,78 (DNA with 8 base pairs), PDB 2JYK4,79 (DNA with 21 base pairs), PDB 1MNV6,80 (Actinomycin D bound to DNA), PDB 8WL233,81 (flagellar motor), PDB 7PKR34,82 (vault protein complex). Structures used for comparing AlphaFold’s pLDDT scores to atomic energies include: PDB AF_AFA0A023FF81F1 (Evasin P1126), PDB AF_AFA0A023IKK2F1 (MHC Class II Beta Chain) and PDB AF_AFA0A023GPK8F1 (Friend of Echinoid, Isoform H). These were obtained from the Computed Structure Models (CSM) section of the RCSB PDB42. Visualisations of biomolecular structures were created using ChimeraX29,69,83,84. All other figure elements (including biomolecular structures at atomic resolution, electronic densities, absorption spectra and atomic energies) were generated with Matplotlib74 for python. References
  46. 48. Liu, X., Patel, A., & Chow, E. (2014, May). A new scalable parallel algorithm for Fock matrix construction. In: 2014 IEEE 28th international parallel and distributed processing symposium (pp. 902–914). IEEE. · DOI
  47. 49. Umeda, H., et al. (2010). Parallel Fock matrix construction with distributed shared memory model for the FMO-MO method. Journal of computational chemistry, 31(13), 2381–2388 · DOI
  48. 50. Chow, E., Liu, X., Smelyanskiy, M., & Hammond, J. R. (2015). Parallel scalability of Hartree–Fock calculations. The Journal of chemical physics, 142(10). · DOI
  49. 51. Petke, J. D., Whitten, J. L., & Douglas, A. W. (1969). Gaussian Lobe Function Expansions of Hartree–Fock Solutions for the Second-Row Atoms. The Journal of Chemical Physics, 51(1), 256–262. · DOI
  50. 52. Helgaker, T., & Taylor, P. R. (1995). Gaussian basis sets and molecular integrals. Modern Electronic Structure Theory: Part II, 725–856. · DOI
  51. 53. Gill, P. M. (1994). Molecular integrals over gaussian basis functions. In: Advances in quantum chemistry (Vol. 25, pp. 141–205). Academic Press. · DOI
  52. 54. Dupuis, M., Rys, J., & King, H. F. (1976). Evaluation of molecular integrals over Gaussian basis functions. The Journal of Chemical Physics, 65(1), 111–116. · DOI
  53. 55. Pople, J. A., & Hehre, W. J. (1978). Computation of electron repulsion integrals involving contracted Gaussian basis functions. Journal of Computational Physics, 27(2), 161–168. · DOI
  54. 56. Yasuda, K. (2008). Two-electron integral evaluation on the graphics processor unit. Journal of Computational Chemistry, 29(3), 334–342. · DOI
  55. 57. Pritchard, B. P., Altarawy, D., Didier, B., Gibson, T. D., & Windus, T. L. (2019). New basis set exchange: An open, up-to-date resource for the molecular sciences community. Journal of chemical information and modeling, 59(11), 4814–4820. · DOI
  56. 58. Virtanen, P., et al. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature methods, 17(3), 261–272.
  57. 59. Abramowitz, M., & Stegun, I. A. (Eds.). (1968). Handbook of mathematical functions with formulas, graphs, and mathematical tables (Vol. 55). US Government printing office.
  58. 60. Neese, F., Wennmohs, F., Becker, U., & Riplinger, C. (2020). The ORCA quantum chemistry program package. The Journal of chemical physics, 152(22). · DOI
  59. 61. Neese, F. (2018). Software update: the ORCA program system, version 4.0. Wiley Interdisciplinary Reviews: Computational Molecular Science, 8(1), e1327. · DOI
  60. 62. Neese, F. (2022). Software update: The ORCA program system—Version 5.0. Wiley Interdisciplinary Reviews: Computational Molecular Science, 12(5), e1606. · DOI
  61. 63. Neese, F. (2023). The SHARK integral generation and digestion system. Journal of Computational Chemistry, 44(3), 381–396. · DOI
  62. 64. Neese, F. (2000). Approximate second-order SCF convergence for spin unrestricted wavefunctions. Chemical Physics Letters, 325(1–3), 93–98. · DOI
  63. 65. Humphrey, W., Dalke, A., & Schulten, K. (1996). VMD: visual molecular dynamics. Journal of molecular graphics, 14(1), 33–38. · DOI · Paper page
  64. 66. Wu, X., Xie, X., & Troisi, A. (2024). Calibration of several first excited state properties for organic molecules through systematic comparison of TDDFT with experimental spectra. Journal of Materials Chemistry C, 12(46), 18886–18892. · DOI
  65. 67. Jacquemin, D., Perpete, E. A., Scuseria, G. E., Ciofini, I., & Adamo, C. (2008). TD-DFT performance for the visible absorption spectra of organic dyes: conventional versus long-range hybrids. Journal of chemical theory and computation, 4(1), 123–135. · DOI
  66. 68. Savitzky, A., & Golay, M. J. (1964). Smoothing and differentiation of data by simplified least squares procedures. Analytical chemistry, 36(8), 1627–1639. · DOI · Paper page
  67. 69. Pettersen, E. F., et al. (2004). UCSF Chimera—a visualization system for exploratory research and analysis. Journal of computational chemistry, 25(13), 1605–1612. · DOI
  68. 70. Van Rossum, G., & Drake Jr, F. L. (1995). Python tutorial (Vol. 620). Amsterdam, The Netherlands: Centrum voor Wiskunde en Informatica.
  69. 71. Lam, S. K., Pitrou, A., & Seibert, S. (2015, November). Numba: A llvm-based python jit compiler. In: Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC (pp. 1–6). · DOI
  70. 72. Harris, C. R., et al. (2020). Array programming with NumPy. Nature, 585(7825), 357–362. · DOI · Paper page
  71. 73. Paszke, A. (2019). Pytorch: An imperative style, high-performance deep learning library. arXiv preprint arXiv:1912.01703. · arXiv
  72. 74. Hunter, J. D. (2007). Matplotlib: A 2D graphics environment. Computing in science & engineering, 9(03), 90–95. · DOI · Paper page
  73. 75. Royant, A., et al. Sensory rhodopsin II (2001) https://doi.org/10.2210/pdb1H68/pdb · DOI
  74. 76. Dominik, H., Karel, S., Fuzik, T., Plevka, P., Structure of native bacteriophage P68 (2018) https://doi.org/10.2210/pdb6Q3G/pdb · DOI
  75. 77. Schafer, M., Sheldrick, G.M., Bahner, I., Lackner, H., Crystal structure of Actinomycin D (1998) https://doi.org/10.2210/pdb1A7Y/pdb · DOI
  76. 78. Thota, N., Li, X.H., Bingman, C.A., Sundaralingam, M., High resolution refinement of the hexagonal A-DNA octamer D(GTGTACAC) at 1.4 Angstroms resolution (1992) https://doi.org/10.2210/pdb1D79/pdb · DOI
  77. 79. Masliah, G., Mauffret, O., NMR Structure of a 21 bp DNA duplex preferentially cleaved by Human Topoisomerase II (2007) https://doi.org/10.2210/pdb2JYK/pdb · DOI
  78. 80. Hou, M.-H., Robinson, H., Gao, Y.-G., Wang, A.H.-J., Actinomycin D binding to ATGCTGCAT (2002) https://doi.org/10.2210/pdb1MNV/pdb · DOI
  79. 81. Tan, J.X., Zhang, L., Zhou, Y., Zhu, Y.Q., Cryo-EM structure of the membrane-anchored part of the flagellar motor-hook complex in the CW state. (2023) https://doi.org/10.2210/pdb8WL2/pdb · DOI
  80. 82. Guerra, P. et al., Vault structure in primmed conformation (2021) https://doi.org/10.2210/pdb7PKR/pdb · DOI
  81. 83. Pettersen, E. F., et al. UCSF ChimeraX: Structure visualization for researchers, educators, and developers. Protein science, 30(1), 70–82. · DOI
  82. 84. Goddard, T. D., et al. (2018). UCSF ChimeraX: Meeting modern challenges in visualization and analysis. Protein science, 27(1), 14–25. · DOI · Paper page
  83. 85. Sehnal, D., et al. (2021). Mol* Viewer: modern web app for 3D visualization and analysis of large biomolecular structures. Nucleic acids research, 49(1), 431–437. · DOI · Paper page

Metadata

Authors
Garcia, M., & Wieners, L.
Affiliations
  • Universität Kassel
  • University of Kassel (DE)
ORCID
Venue
Communications Chemistry
Publisher
Springer Science and Business Media LLC
Published
2025-09-16
Volume
9
Issue
1
ISSN
2399-3669 (ISSN-L)
DOI
10.21203/rs.3.rs-7327472/v1
PMC
PMC12458572
OpenAlex
W7159083704
Cited by
1
References
2
Citations by year
2026: 1
FWCI
2.56
Type
journal article
Language
en
Open access
Licence
Locations
  • Communications Chemistry (publishedVersion, cc-by)
  • Research Square (acceptedVersion, cc-by)
  • PubMed (publishedVersion)
  • PubMed Central (publishedVersion, cc-by)
  • Europe PMC (PubMed Central) (submittedVersion, cc-by)
  • DOAJ (DOAJ: Directory of Open Access Journals) (submittedVersion, cc-by-sa)
  • Kobra (Universitätsbibliothek Kassel) (cc-by)
Funders
  • Deutsche Forschungsgemeinschaft (RTG 2749/1, RTG 2749, GRK2749-1)
  • Technische Universität Darmstadt
  • Universität Kassel
  • National Institutes of Health
  • University of California, San Francisco (R01-GM129325)
  • National Institute of General Medical Sciences (R01GM129325)
  • National Institute of Allergy and Infectious Diseases
Topics
  • Physical Sciences › Physics and Astronomy › Atomic and Molecular Physics, and Optics › Advanced Chemical Physics Studies
  • Physical Sciences › Chemistry › Physical and Theoretical Chemistry › Advanced Physical and Chemical Molecular Interactions
  • Life Sciences › Biochemistry, Genetics and Molecular Biology › Molecular Biology › Protein Structure and Dynamics
Keywords
divide-and-conquer, minimal basis sets, computational efficiency

About these checks

Registry, retraction, Cites and statistics results are computed, not judged, with one exception: a printed reference that no index lists is read once by a model before it counts as a citation error. Claims are extracted by a model and kept only when their quote appears verbatim in the text; each cited source is then checked, with a model’s help, for whether it says what the claim relies on. Those verdicts can be wrong — read the claim’s quote and the cited source before relying on one.

Statistics and arithmetic are recomputed, not judged.

Confidence bands are shown only for fields calibrated on our annotated test set. None is calibrated yet, so no band is shown.

What counts as a claim ›How we measure confidence ›How we check this ›

Are you an author, or do you see a problem with this page? Request a review

How to cite

Cite CoChat for the checks on this page, and the paper for its own findings.

This page
CoChat. (2026). Citation check: Bridging Quantum Mechanics and Biology at the Million-Atom Scale. https://app.cochat.ai/papers/10.21203/rs.3.rs-7327472/v1
The paper
Garcia, M., & Wieners, L. (2025). Bridging Quantum Mechanics and Biology at the Million-Atom Scale. Communications Chemistry. https://doi.org/10.21203/rs.3.rs-7327472/v1