Modular implementation of royalsocietypublishing.org/journal/rsos the linear- and cubic-scaling orbital minimization methods Research in electronic structure codes Cite this article: Lebedeva IV, García A, Artacho using atomic orbitals Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 E, Ordejón P. 2023 Modular implementation of the linear- and cubic-scaling orbital minimization methods in electronic structure codes using atomic Irina V. Lebedeva1,2,3, Alberto García4, orbitals. R. Soc. Open Sci. 10: 230063. Emilio Artacho1,5,6,7 and Pablo Ordejón2 https://doi.org/10.1098/rsos.230063 1 CIC nanoGUNE BRTA, Donostia-San Sebastián 20018, Spain 2 Catalan Institute of Nanoscience and Nanotechnology—ICN2 (CSIC and BIST), Campus UAB, Bellaterra 08193, Spain 3 Simune Atomistics, Avenida de Tolosa 76, Donostia-San Sebastián 20018, Spain Received: 19 January 2023 4 Institut de Ciència de Materials de Barcelona (ICMAB-CSIC), Bellaterra 08193, Spain Accepted: 29 March 2023 5 Donostia International Physics Center DIPC, Donostia-San Sebastián 20018, Spain 6 Theory of Condensed Matter, Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, UK 7 Ikerbasque, Basque Foundation for Science, Bilbao 48011, Spain Subject Category: IVL, 0000-0002-2880-0275; AG, 0000-0001-5138-9579; EA, 0000-0001-9357-1547; PO, 0000-0002-2353-2793 Physics and biophysics We present a code modularization approach to design efficient and Subject Areas: massively parallel cubic- and linear-scaling solvers for electronic computational physics structure calculations using atomic orbitals. The modular implementation of the orbital minimization method, in which Keywords: linear algebra and parallelization issues are handled via external density functional theory, linear-scaling methods, libraries, is demonstrated in the SIESTA code. The distributed block compressed sparse row (DBCSR) and scalable linear algebra code modularization package (ScaLAPACK) libraries are used for algebraic operations with sparse and dense matrices, respectively. The MatrixSwitch and libOMM libraries, recently developed within the Electronic Author for correspondence: Structure Library, facilitate switching between different matrix Irina V. Lebedeva formats and implement the energy minimization. We show results comparing the performance of several cubic-scaling algorithms, e-mail: liv_ira@hotmail.com and also demonstrate the parallel performance of the linear-scaling solvers, and their supremacy over the cubic-scaling solvers for insulating systems with sizes of several hundreds of atoms. 1. Introduction The success of electronic structure theory [1] in modelling new Electronic supplementary material is available materials and devices [2,3] has stimulated the development of online at https://doi.org/10.6084/m9.figshare.c. 6601948. © 2023 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited. hundreds of electronic structure codes [4]. Historically almost all of these individual software packages 2 are written in distinct ways, although many tasks performed are similar. Except for numerical and royalsocietypublishing.org/journal/rsos performance-related libraries such as basic linear algebra subroutines (BLAS) [5], higher-level linear algebra utilities (serial linear algebra package—LAPACK [6] and its parallel counterpart, scalable LAPACK—ScaLAPACK [7]), message passing interface (MPI) level [8] etc., significant parts of the codes are replicated with some variations. Electronic structure packages are growing rapidly, incorporating more and more new features. Also, the codes have to adapt to the constant hardware evolution, which in the case of monolithic code architecture implies significant efforts on re- engineering. In this situation, it seems more efficient to change the traditional monolithic paradigm of software development to the modular one in which common tasks arise [9]. In addition, such an approach allows to separate the tasks related to high-level routines focused on the calculation of physical properties from the implementation of the underlying routines for parallelization and algebra. On one hand, this means that the implementation of new models and Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 algorithms becomes much simpler and does not require the knowledge of technical details related to R. Soc. Open Sci. 10: 230063 parallelization. On the other hand, much better performance is achieved using specialized external libraries, and thus much larger systems can be modelled. Significant efforts (e.g. the European project MAX [10]) are underway to stimulate the paradigm change in software design and to facilitate exascale computing. Here, we show the benefits of modularization by the example of SIESTA [11–15]. SIESTA [11–15] was specifically designed for linear-scaling calculations [16,17] in which the computational time grows linearly with the number of atoms [16–21]. Such methods make possible calculations of large systems at a considerably less computational cost compared with common cubic- scaling approaches. SIESTA uses strictly localized atomic-like functions for basis sets in which the Hamiltonian and overlap matrices, H and S, are sparse. If, additionally, the confinement of the wave functions is imposed, the coefficient matrix C expanding wave functions in the basis is also sparse. Reducing the problem of solving the Kohn–Sham equations to the minimization of a properly constructed energy functional within the Ordejón–Mauri [16–18,22] and Kim [23] approaches, the inversion of the overlap matrix is avoided and only expressions involving products and sums of sparse H, S and C matrices need to be computed, all in linear-scaling effort. The linear-scaling solvers in SIESTA, although available from the start [12], are not widely used in practice. One of the reasons is that the implementation of these physical methods involved also coding of the algebra and parallelization of sparse matrices, which inevitably increased the code complexity and hindered progress. Recent efforts on linear-scaling methods have produced the distributed block compressed sparse row (DBCSR) library [24] that efficiently handles algebraic operations for sparse matrices and is massively parallelized [25,26]. Using this library, we have implemented an improved and more reliable version of linear-scaling solvers in SIESTA (figure 1). Another recent initiative that has helped re-designing SIESTA is the Electronic Structure Library (ESL) [9,31], a collaboration platform for shared software development. We use ESL’s libOMM library [9,27] distributed within the omm-bundle [32]. It encodes the Ordejón–Mauri [16–18,22] and Kim [23] functionals, originally without the additional approximation of wave function confinement, rendering dense C matrices and cubic scaling. Such an approach provides an alternative to conventional cubic- scaling methods, which can be faster in long simulations by avoiding computationally expensive orthonormalization and using history on previous steps [28]. We refer to unconstrained minimization methods of suitable energy functionals, with either linear or cubic scaling, as the orbital minimization method (OMM) [21,28,33,34]. In libOMM [9,27,28], the minimization is customarily performed via conjugate gradients (CG). The parameters of the quartic function describing the energy dependence along the search direction are computed analytically [17,28]. Although the original libOMM library provides cubic scaling [9,27,28], it has been straightforward to extend it to linear scaling, the equations being almost the same, the key difference being the use of sparse matrices instead of dense. Normally two separate pieces of the code dealing with sparse and dense matrices would be used for the same equations. This code duplication can be avoided in libOMM thanks to the MatrixSwitch (MS) library [29], an interface between high-level physical routines and low-level routines for matrix algebra. MS, which is also distributed within the omm-bundle [32] of ESL [9,31], simplifies the coding of matrix operations and allows a single code independent of matrix format, by means of format-independent high-level commands. Depending on matrix format, MS calls the appropriate linear algebra library. An example of calculations using the MS library is shown in listing 1 (see electronic supplementary material for MS overview). Recently, the MS library was extended to support sparse matrices [9,30] via the DBCSR [24–26] library. Here, we consider dense and sparse matrices in the parallel-distributed dense block cyclic 3 SIESTA royalsocietypublishing.org/journal/rsos ESL MatrixSwitch libOMM dense sparse ScaLAPACK DBCSR Figure 1. The use of libraries within the revised orbital minimization method (OMM) solver in the electronic structure code SIESTA [11–15]. The red rectangular box corresponds to SIESTA. Blue ellipses indicate the libraries used [7,9,24,25,27–30]. The libraries in Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 R. Soc. Open Sci. 10: 230063 the dashed frame belong to the Electronic Structure Library (ESL) [9,31]. The arrows demonstrate calls to the libraries. Listing 1. An example of the calculation of the total charge in the OMM approach using the MatrixSwitch library. ( pddbc) and parallel-distributed compressed sparse row ( pdcsr) MS formats for which algebraic operations are handled with the help of the ScaLAPACK [7] and DBCSR [24] libraries, respectively. Although basic functionality for sparse matrices was already provided in this recent MS version [9,30], a revision of the library was needed towards treating sparse and dense matrices on the same footing, and to enable linear-scaling calculations. The incorporation of the solver library into an electronic structure code also implies additional matrix manipulations such as conversions between the matrix formats supported by the code and the solver library as well as reading and writing of restart files. The corresponding subroutines have been here implemented in MS and are discussed below. After a brief overview of the OMM approaches, the new implementation of linear- and cubic-scaling OMM 4 in SIESTA is presented, including the necessary changes in the MS and libOMM libraries forming part of ESL. royalsocietypublishing.org/journal/rsos The results of the first tests are discussed, and recommendations on the efficient use of OMM are given. 2. Overview of orbital minimization method approaches In density functional theory (DFT) [35,36], the problem of finding the ground state of a many-electron system is reduced to an energy minimization for the system of 2n non-interacting electrons moving in an effective potential and described by one-particle states {|ψi〉} (i = 1, …, n) each of which is occupied by two electrons of opposite spin (assuming no spin polarization, for simplicity). The set of states {|ψi〉} is one of the many possible bases in the occupied subspace of the Hilbert space of the system and can be chosen orthonormal or not. In the latter case [37], the overlap matrix SW with the elements Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 (SW)ij = 〈ψi|ψj〉 is not the identity matrix (SW ≠ I, (IW)ij = δij ) and the density matrix operator that R. Soc. Open Sci. 10: 230063 determines the projection onto the occupied subspace is then given by X n ^r ¼ 2 jci i(SW 1 )ij hcj j, ð2:1Þ i,j¼1 involving the inverse of SW [17,18,22]. Here and below, we limit our consideration to insulating systems. The linear-scaling methods applicable to metals are discussed e.g. in [38–41]. The corresponding band structure energy becomes [16–18,22] ^ r] ¼ 2Tr[SW 1 HW ], E ¼ Tr[H^ ð2:2Þ where H ^ is the Hamiltonian operator, and HW is the matrix with the elements ðHW Þ ¼ hci jHj ^ cj i. Note that ij the traces in equation (2.2) are taken on spaces of different dimensions: the size of the basis set for the first, and of the occupied states in the second. Also, the second equality holds for zero temperature. In the basis of m functions {|ϕi〉} (strictly localized atomic orbitals in SIESTA) X m jci i ¼ Cmi jfm i, ð2:3Þ m¼1 where we refer to C as the coefficient matrix. Then HW ¼ Cy HC and SW ¼ Cy SC, where Hij ¼ hfi jHj ^ fj i, y Sij = 〈ϕi|ϕj〉 and C is the Hermitian conjugate of C. The energy functional in equation (2.2) is minimized to find the ground state energy. The most common approach is direct diagonalization of the Hamiltonian matrix H (an m × m matrix for the basis set of size m). Energy and charge density are then obtained using the wave functions and energies of the n lowest eigenstates. By contrast, in the iterative approaches [42], the energy is minimized with respect to variations in the states {|ψi〉}. Here one needs to calculate the inverse of the overlap matrix SW −1 or impose the orthonormality condition (SW)ij = δij. In any case, the computational time increases as O(n 3) with the system size, while the memory required to store the wave functions grows as O(n 2). In OMM approaches [16–18,22], the expensive orthonormalization step is avoided via the modification of the energy functional in such a way that it automatically induces the orthonormalization of the wave functions during minimization ~ ¼ 2Tr[ðIW þ ðIW  SW ÞÞHW ] E ¼ 2Tr[ð2IW  Cy SCÞCy HC]: ð2:4Þ This expression can be derived from consideration of Lagrange multipliers [16,17] or expansion of the inverse overlap matrix to first order in the deviation from the identity [18,22]: SW −1 ≈ IW + (IW − SW). The solution obtained from equation (2.4) is the same as from equation (2.2). Within the same approximation, the density matrix of equation (2.1) is computed as [17] r ¼ C(IW þ ðIW  SW Þ)Cy ¼ 2Cð2IW  Cy SCÞCy ð2:5Þ and the forces on atom I as [17]     @H @S FI ¼ Tr r þ Tr rE , ð2:6Þ @ RI @ RI where we refer to rE ¼ 2CHW Cy as the ‘energy density’. If the basis functions and wave functions are chosen to be strictly localized, the Hamiltonian, overlap 5 and coefficient matrices, H, S and C, are sparse and O(n) scaling with system size is achieved [16–18,22]. Note that this is not the case for equations (2.1) and (2.2) as the inverse of S is not sparse (although sub- royalsocietypublishing.org/journal/rsos cubic scaling can be achieved using selected inversion to compute just the needed elements of the inverse [43]). In the case of periodic systems, localized wave functions are close to the Wannier functions that decay exponentially with the distance from the centre of localization in insulators and in metals at a finite temperature. Imposing localization constraints on the wave functions, however, leads to a deviation from the exact solution of equations (2.2) and (2.4). Also the localized wave functions obtained are not strictly orthonormal and do not comply with the system symmetries [23]. However, the degree of approximation can be controlled with the cut-off radius RC for the wave functions. Here, we limit our consideration to insulators with a substantial band gap, where RC of several Å is normally enough [17,18]. In the Ordejón–Mauri functional [16–18,22], the localization of the wave functions gives rise to many Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 shallow local minima and flat regions in which the algorithm can be trapped for a long time during the R. Soc. Open Sci. 10: 230063 energy minimization. This problem is solved in the Kim functional [23] by including unoccupied states and introducing a chemical potential η, i.e. the energy separating occupied and unoccupied states. The corresponding functional is obtained by (i) an eigenspectrum shift H → H − ηS, (ii) changing 0 0 dimensions of C from m × n to m × n , where n > n, and (iii) changing the energy functional in equation (2.4) as E ~!E ~ þ hn, and energy density ρE in equation (2.6) as ρE → ρE + ηρ. It should be noted, however, that although the multiple-minima problem is solved in the Kim functional, it is sometimes hard to choose a proper value for η. It should always lie within the band gap, but the bands can move up and down during self-consistency or molecular dynamics (MD), η possibly getting into the valence or conduction bands and, as a result, converging to an erroneous solution. Care should be taken to ensure that the solution reproduces the correct number 2n of electrons. If the localization constraints on the wave functions are removed, the exact solution of equations (2.2) and (2.4) is obtained [28]. Even in this case, however, one energy minimization can demand many CG iterations. This relates to the problem of length-scale or kinetic energy ill-conditioning [42,44]. The efficiency of the CG algorithm depends on the ratio of the maximal and minimal extremal curvatures of the function minimized, which in OMM are determined by the maximal and minimal eigenvalues of the Hamiltonian. The eigenspectrum of the Hamiltonian is broad, given the large kinetic energy of high-energy eigenstates. Although such states contribute negligibly to the ground state solution, the problem becomes ill- conditioned and the convergence is slow. It is, however, possible effectively to reduce the width of the eigenspectrum by suppressing the kinetic energy contribution of high-energy states through preconditioning [45,46], by which the CG gradient matrix is multiplied by the preconditioning matrix [28]   T 1 P¼ Sþ , ð2:7Þ tT where τT is the scale for kinetic energy preconditioning and T is the kinetic energy matrix. Another approach for improving the efficiency of CG minimizations is reducing the generalized eigenvalue problem to the standard form via the Cholesky factorization [28]. Both of these approaches involve matrices that are not sparse (the preconditioning matrix or the reduced Hamiltonian) and are considered here only for cubic- scaling OMM. 3. Modular solver architecture 3.1. Solver input and output A scheme of the implemented OMM solver is shown in figure 2. At each self-consistent-field (SCF) step, the solver receives as an input the Hamiltonian and overlap matrices in the basis of strictly localized atomic orbitals, H and S, and the information on the system geometry. SIESTA uses for matrices the standard compressed sparse row format, that is the matrix information is stored in local one- dimensional arrays containing data values and column indices for individual non-zero elements of local rows as well as indices of the first non-zero elements and numbers of non-zero elements for each local row. The blocks of rows are distributed on a one-dimensional process grid (figure 4a). Here and in MS we refer to this format as pdrow to distinguish from the pdcsr format supported by DBCSR. H and S are received by the solver in the pdrow format. The density matrix ρ is the output, also in pdrow (see equation (2.5)). This matrix is used to update H for the next SCF step outside the solver. At the 6 OMM solver royalsocietypublishing.org/journal/rsos geometry H, S in SIESTA format new MD step? matrix format C† initialization reading conversion: H, S or update C† restart MD step no.1 H, S, H, S geometry in MS format update solver library guess C† (libOMM) Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 converged C† in MS format R. Soc. Open Sci. 10: 230063 ρ in MS format matrix format writing conversion: ρ C† restart ρ in SIESTA format Figure 2. Scheme of the revised OMM solver. Blocks within the solver are shown in blue. The rest of the SIESTA code is shown as a red block. Arrows indicate data flow. Hamiltonian, overlap and density matrices are denoted as H, S and ρ. The Hermitian conjugate of the coefficient matrix of expansion of the wave functions in the basis of localized atomic orbitals is denoted as Cy . The Hamiltonian and overlap matrices are converted from internal SIESTA to MatrixSwitch (MS) format for further calculation of ρ and Cy with the help of the libOMM library. The restart file for Cy can be read once at the first molecular dynamics (MD) step. end of each MD step, the solver is called again to compute the energy density matrix ρE, which, along with ρ, is later used to calculate forces (see equation (2.6)) and stresses. The scheme of the ρE calculation is analogous to that of ρ shown in figure 2. 3.2. Solver library The solver uses the libOMM library [9,27,28,32] to perform the CG minimization of the energy functional given by equation (2.4). As an input, the libOMM library requires H and S, as well as the initial guess for Cy , in one of the MS formats [9,29,30,32]. As an output, it provides the converged Cy , and ρ or ρE in the same format. The pddbc format is used for parallel calculations with dense matrices. In this case, all matrix elements are stored and algebraic operations are performed using the ScaLAPACK library [7]. The matrix is divided into two-dimensional blocks distributed on a two-dimensional or one- dimensional process grid. For parallel calculations with sparse matrices, the pdcsr format is used. The matrix is also divided into two-dimensional blocks distributed on a one-dimensional or two- dimensional process grid (figures 4b,c, respectively). However, in this case, zero blocks are not stored. The algebraic operations are performed by the DBCSR library [24–26]. At the moment, libOMM supports only equal rectangular blocks. The equations implemented in the libOMM library are compatible with all OMM flavours discussed in the previous section, including the Ordejón–Mauri and Kim functionals, with and without localization constraints. However, to make the libOMM library functional for sparse matrices, some parts to the code have been reformulated. Now block-size information is passed to the MS library during the allocation of intermediate matrices required for the CG minimization using m_allocate() (see electronic supplementary material). Also, sparsity is imposed on the gradient matrix G (with the elements ~ @ ðCm Þ ) following the sparsity pattern of the initial guess for C. Already during G Gmi ¼ @ E= i calculation [28], only matrix elements that fit into the sparsity pattern are computed in the contributions to G that are given by products of matrices (using keep_sparsity = true option of mm_multiply()). In the rest of the contributions, non-zero elements that do not fit into the sparsity pattern are omitted and no longer stored, while zero elements within the sparsity pattern are stored as zeros. The sparsity of the density (ρ) and energy density (ρE) matrices is assumed to be the same as of the overlap matrix S and only elements of these matrices that fit into the sparsity pattern are pdcsr (sparse) pddbc (dense) 7 royalsocietypublishing.org/journal/rsos m_register_pdrow (A_pdrow,...) The pointers to the arrays and the block size of A_pdrow are passed to MatrixSwitch An array of column indices in the increasing order is prepared to facilitate format conversion m_copy (A, A_pdrow) 1. An intermediate matrix distributed on the 1D process grid and keeping the block size of A_pdrow is created: Arrays with row and column indices of local non-zero blocks are prepared Non-zero elements are Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 Local non-zero blocks of the pdcsr matrix filled into local two-dimensional R. Soc. Open Sci. 10: 230063 are reserved via m_reserve_blocks ( ) arrays including all elements of the local rows Local non-zero blocks are filled in one by one via m_set_element ( ) 2. The intermediate matrix is redistributed into the final matrix A with the requested block sizes and distribution over computational nodes: via dbcsr_complete_redistribute ( ) via pdgemr2d ( ) Figure 3. A series of calls to the MatrixSwitch library required for the format conversion of matrix A from the pdrow format used in SIESTA (A_pdrow) to the pdcsr and pddbc MatrixSwitch formats (A) handled with the DBCSR and ScaLAPACK libraries, respectively. (a) (b) (c) core 1 1 2 core 2 core 3 core 4 pdrow matrix on 1D MPI grid pdcsr matrix on 1D MPI grid pdcsr matrix on 2D MPI grid Figure 4. Example of matrix format conversion from the pdrow format used in SIESTA to pdcsr MatrixSwitch format handled with DBCSR: (a) pdrow matrix distributed on the one-dimensional process grid with four CPU cores, (b) pdcsr matrix with 2 × 3 blocks distributed on the same one-dimensional process grid and (c) pdcsr matrix with 2 × 3 blocks distributed on the 2 × 2 two- dimensional process grid. Arrows indicate steps 1 and 2 of subroutine m_copy() as explained in figure 3. Small squares represent elements of the 8 × 12 matrix. Black squares are zero elements that are not stored. Red, yellow, blue and green squares correspond to elements stored on cores 1, 2, 3 and 4, respectively. computed. Additionally, the expression for the calculation of ρE has been corrected as compared with the previous libOMM version [28] in accordance with equations (2.4)–(2.6) and [17]. The Cholesky factorization and kinetic energy preconditioning are available only for dense matrices. 3.3. Cy matrix format conversion In order to incorporate the libOMM library into SIESTA within the OMM solver, the following steps are required (figure 2): (1) matrix format conversion from/to the SIESTA format to/from the MS formats and (2) initialization and update of Cy , according to the current geometry of the system. The matrix format conversion is realized using calls to MS subroutines m_register_pdrow() and m_copy() (see electronic supplementary material). The first of this subroutines has been added to the MS library and the second one has been extended to allow the conversion from/to the pdrow format to/from the pdcsr and pddbc formats. The conversion is performed as follows (figure 3). First the pointers to arrays of the pdrow matrix and its block size are passed to MS. Then a pdcsr/pddbc matrix distributed on the one-dimensional process grid with the same block size for rows as the initial pdrow matrix is filled in element by element (figure 4). The missing elements of the pddbc matrix or within non-zero blocks of the pdcsr matrix are filled with zeros. Note that to speed up the conversion and guarantee linear scaling, column and row indices of non-zero blocks of the pdcsr matrix should be passed to the DBCSR library before filling the values 8 via the call to m_reserve_blocks() (see electronic supplementary material). Once the one-dimensional- royalsocietypublishing.org/journal/rsos distributed pdcsr/pddbc matrix is ready, it can be redistributed on a two-dimensional process grid. In the case when the final matrix is distributed on the one-dimensional process grid and has the same block size for rows at the initial pdrow matrix, the last step is omitted. The conversion from pdcsr and pddbc to pdrow is implemented in a similar way. It is assumed that the row and column indices of non-zero elements of the pdrow matrix are already known. Only the values of the matrix elements are restored. 3.4. Cy matrix initialization and update The initialization of the Cy matrix in the sparse form is performed in SIESTA in the following way. It is supposed that each atom carries the number of localized wave functions equal to the atomic charge (in Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 units of elementary charge) divided by two, Qat/2. If Qat is odd, (Qat + 1)/2 localized wave functions are R. Soc. Open Sci. 10: 230063 assigned to one atom and (Qat − 1)/2 to the next one. This procedure is repeated for all the atoms in the system. Then the Cy matrix in the pdrow format with the total number of rows that corresponds to the total number of localized wave functions, NWF = Q/2, where Q is the sum of atomic charges in the system, is prepared. The local rows are assigned according to the block size bWF . By default, it equals the block size for the basis functions, bBF , multiplied by the ratio of the total number NWF of localized wave functions to the basis set size NBF: bWF = bBFNWF/NBF . For each local row, the local environment of the atom hosting the corresponding localized wave function is analysed. The row elements that correspond to atoms beyond some cut-off radius RC from the atom considered are supposed to be zero. The row elements that correspond to atoms within the cut-off radius RC are initialized by random values. This sparsity pattern is maintained during the energy functional minimization. The Cy matrix in the pdrow format is converted to the pdcsr or pddbc formats in the same manner as the Hamiltonian and overlap matrices, H and S. It should be mentioned also that the initial cut-off radius RC,ini for initialization of the Cy matrix can be set different from RC used for the energy minimization. Choosing a small initial radius RC,ini (several Å) helps to avoid convergence problems and is useful not only in calculations with sparse matrices but also with dense ones. At each new MD step, the sparsity pattern of the Cy matrix is checked again. The elements that now should be zero because the corresponding atoms got away by more than RC are set to zero and no longer stored. The elements corresponding to the atoms that got closer than RC are now stored and treated as non-zero but are assigned to zero as the initial guess. Linear extrapolation of the Cy matrix based on the information from the two previous MD steps is also possible. 3.5. Cy matrix input and output The restart file for the Cy matrix can be written at each SCF step and read at the beginning of the run. These operations are performed by calling new MS subroutines m_read() and m_write(), respectively (see electronic supplementary material). If the Cy matrix in the pdcsr or pddbc format is distributed on a two- dimensional process grid, it is first converted into a one-dimensional-distributed matrix (by analogy with the format conversion routines). Then the blocks of rows are consecutively passed to the head core and written to the file. To read the file, the reverse operations are performed. The block sizes and process grid for the Cy matrix do not need to be the same as used when writing the restart information. Upon reading, the sparsity pattern of the Cy matrix is corrected according to the current system geometry. 3.6. SIESTA input parameters The input parameters for SIESTA corresponding to the revised OMM solver are described in table 1. To use the OMM solver, SolutionMethod should be set to BLOMM (OMM with block matrices). 4. Tests 4.1. Computational details The test calculations have been carried out for single-layer boron nitride (BN) under periodic boundary conditions. Supercells of BN from 12 × 12 to 96 × 96 with up to 18 400 atoms are considered. The lattice Table 1. Principal input parameters for the revised OMM solver in SIESTA (SolutionMethod BLOMM) and their 9 default values. royalsocietypublishing.org/journal/rsos input parameter default value description OMM.UseSparse true whether to use sparse matrices OMM.UseKimFunctional true whether to use the Kim [23] (or Ordejón–Mauri [16–18,22]) functional OMM.Use2D true whether to distribute matrices on a two- dimensional process grid OMM.ReadCoeffs false whether to read the initial localized wave functions (LWFs), i.e. the Cy Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 matrix, from the restart file R. Soc. Open Sci. 10: 230063 ( .WF_COEFFS_BLOMM) OMM.WriteCoeffs false whether to write the LWFs (Cy matrix) to the restart file OMM.RelTol 10−9 the tolerance for the energy convergence in conjugate-gradient (CG) iterations. When 2(En − En−1)/(En + En−1), where En is the energy at CG iteration n, becomes smaller than this tolerance, CG iterations are stopped OMM.BlockSizeC bWF = bBFNWF/NBF the block size for LWFs (rows of the Cy matrix). By default, equals the block size for the basis functions bBF (input parameter BlockSize) multiplied by the ratio of the total number NWF of LWFs to the basis set size NBF OMM.Eta 0 eV the chemical potential for the Kim functional OMM.RcLWF 9.5 Bohr the cut-off radius RC for LWFs determining the sparsity pattern of the Cy matrix OMM.RcLWFInit 0 Bohr the initial cut-off radius RC,ini for LWFs. It is the same as OMM.RcLWF if set to 0 OMM.Extrapolate false whether to estimate LWFs at the next molecular dynamics (MD) step by the linear extrapolation of the results of two last MD steps only for the cubic-scaling OMM OMM.Precon -1 the number of self-consistent-field (SCF) steps for which to apply the preconditioning [28]. If negative, the preconditioning is applied at all SCF steps OMM.PreconFirstStep OMM.Precon OMM.Precon for the first MD step OMM.TPreconScale 10 Ry the scale τT for the kinetic energy preconditioning (see equation (2.7)) OMM.Cholesky false whether to apply the Cholesky factorization [28] 102 10 diagonalization dense, O.-M., preconditioning royalsocietypublishing.org/journal/rsos dense, O.-M., Cholesky 10 dense, O.-M. wall time (hours) 1 sparse, O.-M., RC of 10–1 sparse, O.-M., RC = 4 Å sparse, Kim, RC = 4 Å, η = –5.5 eV 10–2 103 104 Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 no. of atoms, N R. Soc. Open Sci. 10: 230063 Figure 5. Wall time (in hours) for four MD steps for a single boron nitride (BN) layer computed using different approaches versus number N of atoms in the system: (black squares) diagonalization, (green triangles up) OMM with dense matrices (using ScaLAPACK) and preconditioning using a kinetic-energy scale τT = 10 Ry, (dark green triangles down) OMM with dense matrices with Cholesky factorization (open green triangles up) plain OMM with dense matrices, (magenta diamonds) OMM with sparse matrices (using DBCSR) without wave function localization (wave function cut-off radius RC → ∞), (red circles) Ordejón–Mauri functional with RC = 4 Å and (open red circles) Kim functional with RC = 4 Å and chemical potential η = −5.5 eV. In all the cases without wave function localization, the Ordejón–Mauri functional is considered. The calculations are performed on 96 CPU cores. A double-zeta polarized (DZP) basis set is used. The block size is bWF = 6 for the wave functions and bBF = 13 for the basis functions. constant is set at 2.48 Å. The height of the simulation cell is 20 Å. The calculations have been performed at the single G point. The local density approximation [47], norm-conserving Troullier–Martins [48] pseudopotentials and standard built-in double-zeta polarized (DZP) basis set [49] are used. The atomic orbitals are set to zero beyond the cut-off determined by the energy shift of 10 meV (cut-off radii 2.5–4.5 Å). The real-space grid is equivalent to the plane-wave cut-off energy of 100 Ry. The linear mixing scheme with a mixing parameter of 0.1 is applied to converge the ground state. The tolerance is 10−4 for the density matrix and 10−3 eV for the matrix elements of the Hamiltonian. To test performance of different approaches in MD simulations, several MD steps starting from the converged ground state have been computed (the ground state is converged previously with the same method as used for MD). The microcanonical ensemble with an initial temperature of 300 K is considered. The Verlet algorithm [50] with a time step of 1 fs is used. The Pulay mixing scheme [51] with a mixing parameter of 0.2 is applied during the MD simulations. The matrices involved in the calculations consist of equal blocks. For the DZP basis set, each boron and nitrogen atom has 13 basis functions, and hosts three or two wave functions depending on whether the unoccupied states are included into consideration or not, respectively. Therefore, the block size for the wave functions is usually chosen to be bWF = 6 and for the basis functions bBF = 13. The matrices are distributed on a two-dimensional process grid. The cut-off radius for localized wave functions in typical calculations with sparse matrices is RC = 4 Å. The chemical potential for the Kim functional is η = −5.5 eV. CG iterations are performed until the difference of energies at consecutive CG iterations divided by the average energy at these iterations reaches 10−9. The tests with the preconditioning for dense matrices have been carried out using the scale for the kinetic energy of τT = 10 Ry [28]. 4.2. Results To compare the performance of diagonalization and OMM with dense and sparse matrices, we have performed test MD simulations for single-layer BN in different sizes. Figure 5 demonstrates that the approaches in which the wave functions are not confined in space have much worse scaling with system size than the methods with localized wave functions within a cut-off radius RC. The scaling of the former approaches is close to cubic for large systems (exceeding 1000 atoms in our calculations). It should be noted, however, that for small systems (within 1000 atoms) the scaling is sub-cubic. The reason is that for such systems the solver contribution to the total time plotted in figure 5 is comparable to the contributions of other parts of the code that have linear scaling with system size. Among the methods using dense matrices, OMM with applied preconditioning or Cholesky factorization, which improve (a) 11 3 royalsocietypublishing.org/journal/rsos wall time (hours) 2 1 0 0.5 1.0 1.5 2.0 no. of atoms, N (104) Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 (b) R. Soc. Open Sci. 10: 230063 1 H update (DHSCF in SIESTA) 10–1 solver library (libOMM) H, S and ρ format conversion time fraction C†, restart writing 10–2 10–3 103 104 no. of atoms, N Figure 6. (a) Wall time (in hours) for three MD steps with 12 SCF iterations each for a single BN layer described using the Kim functional versus number N of atoms (in 104 atoms). A linear fit is shown by the solid line. (b) Relative contributions of the OMM solver subroutines to the total time versus N: (magenta circles) Hamiltonian update on a three-dimensional grid (DHSCF subroutine in SIESTA [13]), (open purple diamonds) solver library libOMM, (red triangles down) format conversion of H, S and ρ and (open grey triangles right) writing restart for localized wave functions (Cy ). The calculations are performed on 192 CPU cores for a DZP basis, RC = 4 Å, η = −5.5 eV, bWF = 6 and bBF = 13. convergence, shows a slightly better scaling compared with diagonalization or plain OMM. Also OMM using the DBCSR library with no localization of wave functions (RC → ∞) clearly has a better scaling than OMM using ScaLAPACK. This is explained by the fact that the former, although having a dense coefficient matrix, still exploits the sparsity of the Hamiltonian and overlap. In the range of system sizes considered, OMM with kinetic energy preconditioning is the fastest among the approaches without wave function localization, followed by OMM with the Cholesky factorization, diagonalization and plain OMM (figure 5). The crossover between preconditioned dense OMM and the linear-scaling methods takes place for the system with about 1200 atoms. For the plain dense OMM and for diagonalization, the crossovers with linear-scaling methods occur earlier, at about 300 and 700 atoms, respectively. Our timings for single-layer BN have confirmed that the Ordejón–Mauri and Kim approaches in which the wave functions are localized within a cut-off radius RC show linear scaling with system size (figure 6a). The computational times corresponding to different parts of the solver (matrix conversion, libOMM library, initialization and update of the coefficient matrix, reading and writing of restart for localized wave functions) and other parts of the SIESTA code such as the subroutine for the Hamiltonian update called after the density matrix change at each SCF step (DHSCF), all do change linearly upon increasing the system size. As a result, relative contributions of different parts of the code do not depend on the system size (figure 6b). This is different from the cubic-scaling methods, in which the solver very early takes most of the computing time upon increasing the system size, since the rest of the code has linear scaling. It should also be noted that, for the system considered, the solver takes only 40–50% of the computational time, comparable, for example, to the subroutine for the Hamiltonian update (DHSCF in SIESTA). Most of this time (a) 1.4 12 royalsocietypublishing.org/journal/rsos 1.3 relative time 1.2 1.1 1.0 0 10 20 30 40 Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 block size for basis set, bBF R. Soc. Open Sci. 10: 230063 (b) 2.0 1.8 relative time 1.6 1.4 1.2 1.0 0 5 10 15 20 25 block size for wavefunctions, bWF Figure 7. Relative time for the solver library libOMM during four MD steps for a single BN layer with a 60 × 60 supercell (7200 atoms) using the Kim functional versus block size: (a) the block size for the basis set, bBF, is changed and the block size for the wave functions is kept as bWF = 6 and (b) bWF is changed and bBF = 13. The relative time is given with respect to the result for bWF = 6 and bBF = 13. In (a), the dashed line is shown to guide the eye for the data obtained for bBF divisible by the number of basis functions per atom. The calculations are performed on 192 CPU cores for a DZP basis, RC = 4 Å, η = −5.5 eV. corresponds to the minimization of the energy functional given by equation (2.4) performed by the solver library libOMM. The matrix format conversion takes only 0.5–1.0% of the total time. Writing of the restart files for localized wave functions takes up to 0.3% of the time, and initialization and update of the coefficient matrix take a negligible time within 0.01%. The dependence of computational time on block size for the Kim functional with DBCSR are presented in figure 7. In the case of the double-zeta polarized (DZP) basis set, each boron and nitrogen atom hosts three localized wave functions and 13 basis functions. Accordingly the computational time drops significantly at block-size values bBF for the basis functions divisible by 13 (figure 7a). For such block sizes, the computational time grows upon increasing the block size (note that the growth continues beyond the block sizes shown in figure 7a) and has the minimum at bBF = 13. The wave function block-size bWF dependence reaches the minimum at bWF = 6 −10. At small bWF , a fast growth of the computation time is observed. It can be attributed to an increase in the number of non-empty blocks considered upon decreasing the block size. At large bWF , the computational time also grows but at a slower rate. This dependence can be explained by increasing the number of matrix elements that are stored and explicitly considered in matrix operations. Therefore, we find optimal block sizes both for the wave functions and basis functions of the order of 10. Furthermore, chemical considerations can be exploited when dividing matrices into blocks. Still, the optimal choice of block sizes for complex systems is not straightforward and requires further investigation [26]. The CPU scaling of the libOMM solver library in calculations with sparse matrices using DBCSR is shown in figure 8a. A similar CPU scaling is observed for systems of different size (figure 8a), with different block and basis set sizes. The computational time decreases by a factor of about 2.5 upon (a) 1.0 13 royalsocietypublishing.org/journal/rsos 480 0.8 cost (CPU • hours) 288 240 576 192 0.6 432 384 288 576 192 0.4 480 144 96 240 384 288 0.2 192 144 96 0 2 4 6 8 10 12 Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 wall time (s) R. Soc. Open Sci. 10: 230063 (b) 2 7200 atoms, ρE 10 368 atoms, ρE 16 200 atoms, ρE relative cost 1 7200 atoms, CG + ρ 10 368 atoms, CG + ρ 16 200 atoms, CG + ρ 1 2 relative time Figure 8. (a) Computational cost (in CPU·hours) versus wall time (in s), for one call to the solver library libOMM including one conjugate gradient (CG) iteration and calculation of the density matrix ρ, for different supercells of BN using the Kim functional: (black squares) 60 × 60, (red circles) 72 × 72 and (blue triangles) 90 × 90 (7200, 10 368 and 16 200 atoms, respectively). The number of CPU cores used is indicated. (b) Relative computational cost versus relative time for the calls to the solver library including one conjugate gradient iteration and calculation of the density matrix ρ (closed symbols) or calculation of the energy density ρE (open symbols) for different supercells of BN. Relative values are given with respect to the results for 192 CPU cores. The calculations are for a DZP basis, RC = 4 Å, η = −5.5 eV, bWF = 6 and bBF = 13. doubling the computational cost. Such a speed-up is observed for CG energy minimization and subsequent calculation of ρ. It should be noted, however, that calls to libOMM for calculation of ρE involving only two matrix multiplication operations show much better CPU scaling. This can be appreciated from a twice steeper slope of computational cost versus computational time as compared with the calls for energy minimization and calculation of the density matrix (figure 8b). It can, therefore, be expected that the solver parallelization might be further improved via proper code refactoring. The use of OpenMP, GPUs and the library for small matrix multiplication (LIBXSMM) [52] are known to lead to a superior DBCSR performance [25,26], which also requires investigation. 4.3. Recommendations for orbital minimization method solver use The new modular implementation of the OMM solver makes it easier to disentangle technical problems in e.g. parallelization from drawbacks of the OMM itself. Here, we present the first implementation of the solver using external libraries that represents the starting point for further performance improvement and method polishing. Ways to improve the solver performance were mentioned in the previous subsection. We briefly discuss now the drawbacks of the OMM and how they can be addressed. One of the most important methodological problems of the OMM approach is in the minimization, which can require a large number of CG iterations. As shown in figure 9, the first SCF iteration from scratch is rather costly both for the linear and cubic-scaling OMM. For the linear-scaling methods, the first SCF iteration can include thousands of CG steps, followed by tens of SCF iterations with hundreds of CG steps each. After that each SCF step needs just a few CG iterations, becoming very fast. It should sparse (Ordejon–Mauri) 14 sparse (Kim) royalsocietypublishing.org/journal/rsos 103 dense no. of CG iterations dense (preconditioning) dense (Cholesky) 102 10 1 0 10 20 30 40 50 60 SCF step Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 R. Soc. Open Sci. 10: 230063 Figure 9. Number of conjugate-gradient (CG) iterations versus number of self-consistent-field (SCF) steps for the calculation of the ground state of a 12 × 12 BN supercell (288 atoms) from scratch using OMM with sparse matrices: (closed black squares) Ordejón– Mauri and (open black squares) Kim functionals, and OMM with dense matrices using the Ordejón–Mauri functional: (open magenta diamonds) plain, (closed magenta diamonds) preconditioned with kinetic-energy scale of τT = 10 Ry and (blue crosses) with the Cholesky factorization. The calculations are performed on 96 CPU cores for a DZP basis, RC = 4 Å, η = −5.5 eV, bWF = 6 and bBF = 13. Linear mixing with a mixing parameter of 0.1 is used. 1 O.-M. energy (eV atom−1) O.-M. force (eV Å−1) 10–2 deviation 10–4 10–6 Kim energy (eV atom−1) Kim force (eV Å−1) 10–8 0 5 10 15 20 cut-off radius RC (Å) Figure 10. Deviations of energy (in eV atom−1, closed symbols) and force (in eV Å−1, open symbols) for the 60 × 60 supercell of boron nitride (7200 atoms) with atoms displaced by 0.05 Å from their equilibrium positions from the results for the infinite cut-off radius for the wave functions RC → ∞ versus cut-off radius RC (in Å): (black squares) Ordejón–Mauri (O.-M.) and (red diamonds) Kim methods. A DZP basis set is used. The chemical potential for the Kim method is η = −5.5 eV. The block size is bWF = 6 for the localized wave functions and bBF = 13 for the basis functions. be noted that except for the very first SCF iterations, the linear-scaling and plain cubic-scaling OMM require roughly the same numbers of CG steps. However, kinetic energy preconditioning or Cholesky factorization significantly reduce the number of CG iterations required, with a considerable computational-time reduction (see also figure 5). Therefore, it is always recommended to use either of both ways to deal with kinetic energy ill-conditioning in dense OMM. The extension of these approaches to sparse matrices is not straightforward and requires further investigation. Also starting from scratch, one can get into regions in parameter space where the energy functional does not have a minimum in the CG line minimization. To avoid this situation, we recommend using a small cut-off radius RC,ini for the initial guess of wave functions both for linear- and cubic-scaling OMM. It is also recommended to preconverge the ground state using a small linear-mixing parameter. Starting from as low as 0.01 can be required for very large systems. It can then be gradually increased to normal values of 0.1–0.2. After getting close to the ground state, the use of other mixing schemes is possible. If the geometry of the system is far from the optimal one, a reduced step for geometry optimization may also be needed when starting. In figure 10, we address the accuracy of force and energy calculations with the Ordejón–Mauri and Kim functionals for BN. The deviation from the results for the wave functions without localization (RC → ∞) is plotted for different cut-off radii RC. It is seen that for both of the functionals, the accuracy 15 improves upon increasing the cut-off radius in a similar manner. The deviations of the energy and forces within 0.01 eV atom−1 and 0.02 eV Å−1 are achieved already for the cut-off radius of RC = 4 Å. These royalsocietypublishing.org/journal/rsos results confirm that for insulating systems with a substantial band gap, it is sufficient to consider cut- off radii of several Å [17,18]. The Ordejón–Mauri and Kim functionals were designed for insulating systems with a substantial gap. For metals, a smearing function needs to be introduced. However, this is not easy since the information on individual Kohn–Sham eigenstates is missing in OMM. An idea for combining OMM with another method resolving eigenstates close to the Fermi level was proposed in [28] but still requires exploration. Note that modelling of metallic systems requires a much more significant computational effort than modelling of insulators [38,39,41]. As for magnetic systems, the OMM calculations can be performed taking into account spin polarization. At each SCF step, the coefficient matrices for spin up and spin down are found Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 sequentially. All the observations for non-spin-polarized systems discussed above still hold in this case. R. Soc. Open Sci. 10: 230063 5. Conclusion We have demonstrated how modularization simplifies the implementation of new solvers in electronic structure codes by revising the OMM solver in the SIESTA code [11–15]. Matrix algebra operations and parallelization are efficiently handled via external libraries. In particular, the implementation benefits from two ESL [9,31] libraries: libOMM [9,27,28,32] and MS [9,29,30,32]. The libOMM library is used to perform the minimization of the energy functional, while the MS library serves as an interface to low-level algebraic routines facilitating switching between different matrix formats. These libraries have been extended to make possible not only cubic-scaling but also linear-scaling OMM calculations for insulating systems with a substantial band gap. Now the energy functional minimization in libOMM can be carried out for sparse matrices with the DBCSR library [24–26], in addition to dense matrices using ScaLAPACK [7]. To facilitate incorporating libOMM into electronic structure codes based on atomic orbitals, MS has been also supplemented with subroutines for matrix format conversion and matrix reading and writing. The solver library libOMM can be easily further developed in the MS language for the implementation of new solvers. The extended MS and libOMM libraries available through ESL [9,31] can be used for implementation of linear- and cubic-scaling OMM approaches in other codes. The libraries can be used with different types of local basis sets. The only condition for achieving the linear-scaling behaviour is that either the basis functions go to zero beyond some cut-off radius or the elements of the input matrices are filtered with respect to some tolerance to ensure that the matrices are sparse. Note that implementation of custom conversion routines is needed if the matrix format is different from the MS or SIESTA formats. To test the performance of the new OMM and traditional diagonalization solvers available in SIESTA, large-scale calculations have been performed for a BN layer. When sparse matrices and localized wave functions are used, linear scaling with system size is achieved in practice, as expected. Matrix conversion, reading and writing of restart files, as well as initialization and update of the localized wave functions take a small fraction of the computational time. For the linear-scaling methods that fraction does not depend on system size. The cubic-scaling OMM with kinetic energy preconditioning performs best for small systems, even better than diagonalization. For plain OMM, diagonalization and cubic-scaling OMM with kinetic energy preconditioning, the crossovers with linear-scaling methods are observed at about 300, 700 and 1200 atoms, respectively. The best performance for the linear-scaling OMM with sparse matrices is achieved when the wave functions and basis functions are divided into blocks of sizes around 10, taking into account the chemical structure. The OMM solver is MPI-parallelized. When using the DBCSR library [24–26] for algebraic operations with sparse matrices, the computational time decreases by a factor of 2.5 upon doubling the computational cost. It is expected that CPU scaling can be further improved via refactoring some operations in the libOMM library, using OpenMP and GPUs, etc. To perform OMM calculations from scratch, it is recommended to start using a small linear-mixing parameter (down to 0.01), a small step for geometry optimization, and cut-off radii for the wave functions of a few Å. For the cubic-scaling OMM, the convergence becomes much faster with kinetic energy preconditioning or Cholesky factorization. The extension of these approaches to sparse matrices demands further investigation. Data accessibility. The data and relevant code for this research work are stored in GitLab: https://gitlab.com/irina_ 16 lebedeva/siesta/-/tree/orderN (SIESTA), https://gitlab.com/ElectronicStructureLibrary/omm-bundle (omm- royalsocietypublishing.org/journal/rsos bundle) and have been archived within the Zenodo repository: https://doi.org/10.5281/zenodo.7781100 [53] (SIESTA), https://doi.org/10.5281/zenodo.7781174 [54] (MatrixSwitch and libOMM). The raw data for tests have been archived within the Mendeley Data repository: https://doi.org/10.17632/c8kz58bg5z.1 [55]. Supplementary material is available online [56]. Authors’ contributions. I.V.L.: investigation, software, validation, visualization, writing—original draft; A.G.: conceptualization, software, writing—review and editing; E.A.: conceptualization, project administration, resources, supervision, writing—review and editing; P.O.: conceptualization, methodology, project administration, writing— review and editing. All authors gave final approval for publication and agreed to be held accountable for the work performed therein. Conflict of interest declaration. We declare we have no competing interests. Funding. The authors acknowledge the European Union MaX Center of Excellence (EU-H2020 grant no. 824143), the Partnership for Advanced Computing in Europe (PRACE) for awarding us access to computational resources in Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 Joliot-Curie at GENCI@CEA, France (EU-H2020 grant no. 2019215186), computational resources at Pirineus and the R. Soc. Open Sci. 10: 230063 technical support provided by Consorci de Serveis Universitaris de Catalunya (RES grant nos. FI-2022-1-0023 and FI-2022-2-0035) as well as technical and human support provided by IZO-SGI SGIker of the University of the Basque Country (UPV/EHU) and European funding (ERDF and ESF). ICN2 is supported by the Severo Ochoa programme from Spanish MINECO (grant no. CEX2021-001214-S) and by Generalitat de Catalunya (CERCA programme). ICMAB is supported by the Severo Ochoa programme from Spanish MICIU (grant no. CEX2019- 000917-S). P.O. acknowledges support by Spanish MICIU, AEI and EU FEDER (grant no. PGC2018-096955-B-C43). A.G. acknowledges support by Spanish MICIU, AEI and EU FEDER (grant no. PGC2018-096955-B-C44). P.O. and A.G. acknowledge support from Generalitat de Catalunya (grant no. 2021SGR01519). E.A. acknowledges funding from Spanish MICINN through grant no. PID2019-107338RB- C61/AEI/10.13039/501100011033, as well as a María de Maeztu award to Nanogune, grant no. CEX2020-001038-M funded by MCIN/AEI/10.13039/501100011033. Acknowledgements. We also thank Dr David López–Durán for useful discussions and providing the introduction to the MatrixSwitch library. References 1. Martin RM. 2004 Electronic structure: basic 12. Ordejón P, Artacho E, Soler JM. 1996 Self- 19. Galli G. 1996 Linear scaling methods for theory and practical methods. Cambridge, UK: consistent order-N density-functional electronic structure calculations and quantum Cambridge University Press. calculations for very large systems. Phys. Rev. B molecular dynamics simulations. Curr. Opin. 2. Mardirossian N, Head-Gordon M. 2017 Thirty 53, R10441–R10444. (doi:10.1103/PhysRevB.53. Solid State Mater. Sci. 1, 864–874. (doi:10.1016/ years of density functional theory in R10441) S1359-0286(96)80114-8) computational chemistry: an overview and 13. Soler JM, Artacho E, Gale JD, Garcıa A, Junquera J, 20. Goedecker S. 1999 Linear scaling electronic extensive assessment of 200 density functionals. Ordejón P, Sánchez-Portal D. 2002 The SIESTA structure methods. Rev. Mod. Phys. 71, Mol. Phys. 115, 2315–2372. (doi:10.1080/ method for ab initio order-N materials simulation. 1085–1123. (doi:10.1103/RevModPhys.71.1085) 00268976.2017.1333644) J. Phys.: Condens. Matter 14, 2745–2779. (doi:10. 21. Bowler DR, Miyazaki T. 2012 O(N) methods in 3. Krylov A et al. 2018 Perspective: computational 1088/0953-8984/14/11/302) electronic structure calculations. Rep. Prog. chemistry software and its advancement as 14. Sánchez-Portal D, Ordejón P, Artacho E, Soler Phys. 75, 036503. (doi:10.1088/0034-4885/75/3/ illustrated through three grand challenge cases JM. 1997 Density-functional method for very 036503) for molecular science. J. Chem. Phys. 149, large systems with LCAO basis sets. 22. Mauri F, Galli G, Car R. 1993 Orbital formulation 180901. (doi:10.1063/1.5052551) Int. J. Quantum Chem. 65, 453–461. (doi:10. for electronic-structure calculations with linear 4. Molecular Sciences Software Institute (MolSSI). 1002/(SICI)1097-461X(1997)65:5<453::AID- system-size scaling. Phys. Rev. B 47, 2016. See https://molssi.org/. QUA9>3.0.CO;2-V) 9973–9976. (doi:10.1103/PhysRevB.47.9973) 5. BLAS technical forum. Since 1979. See http:// 15. García A et al. 2020 Siesta: recent developments 23. Kim J, Mauri F, Galli G. 1995 Total-energy www.netlib.org/blas/blast-forum. and applications. J. Chem. Phys. 152, 204108. global optimizations using nonorthogonal 6. Anderson E et al. 1999 Lapack users’ guide. (doi:10.1063/5.0005077) localized orbitals. Phys. Rev. B 52, 1640–1648. Philadelphia, PA: SIAM. (doi:10.1137/1. 16. Ordejón P, Drabold DA, Grumbach MP, (doi:10.1103/PhysRevB.52.1640) 9780898719604) Martin RM. 1993 Unconstrained minimization 24. Distributed block compressed sparse row (DBCSR) 7. ScaLAPACK. Since 1992. See http://www.netlib. approach for electronic computations that scales library. See https://github.com/cp2k/dbcsr. org/scalapack. linearly with system size. Phys. Rev. B 48, 25. Borst̆ nik U, VandeVondele J, Weber V, Hutter J. 8. EMPI forum. Since 1991. See https://www.mpi- 14 646–14 649. (doi:10.1103/PhysRevB.48. 2014 Sparse matrix multiplication: the forum.org. 14646) distributed block-compressed sparse row library. 9. Oliveira MJT et al. 2020 The CECAM electronic 17. Ordejón P, Drabold DA, Martin RM, Grumbach MP. Parallel Comput. 40, 47–58. (doi:10.1016/j. structure library and the modular software 1995 Linear system-size scaling methods for parco.2014.03.012) development paradigm. J. Chem. Phys. 153, electronic-structure calculations. Phys. Rev. B 51, 26. Kühne TD et al. 2020 CP2K: an electronic 024117. (doi:10.1063/5.0012901) 1456–1476. (doi:10.1103/PhysRevB.51.1456) structure and molecular dynamics software 10. European Centre of Excellence MAX (MAterials 18. Mauri F, Galli G. 1994 Electronic-structure package—quickstep: efficient and accurate design at the eXascale). 2019. See http://www. calculations and molecular-dynamics electronic structure calculations. J. Chem. Phys. max-centre.eu/. simulations with linear system-size scaling. 152, 194103. (doi:10.1063/5.0007045) 11. SIESTA code. See https://gitlab.com/siesta- Phys. Rev. B 50, 4316–4326. (doi:10.1103/ 27. LibOMM library. See https://esl.cecam.org/ project/siesta. PhysRevB.50.4316) software/libomm/. 28. Corsetti F. 2014 The orbital minimization binding molecular dynamics with localized 48. Troullier N, Martins JL. 1991 Efficient 17 method for electronic structure calculations with orbitals. Phys. Rev. B 51, 9455–9464. (doi:10. pseudopotentials for plane-wave calculations. finite-range atomic basis sets. Comput. Phys. 1103/PhysRevB.51.9455) Phys. Rev. B 43, 1993–2006. (doi:10.1103/ royalsocietypublishing.org/journal/rsos Commun. 185, 873–883. (doi:10.1016/j.cpc. 40. Corkill JL, Ho K-M. 1996 Electronic occupation PhysRevB.43.1993) 2013.12.008) functions for density-matrix tight-binding 49. Junquera J, Paz Ó, Sánchez-Portal D, Artacho E. 29. MatrixSwitch (MS) library. See https://esl.cecam. methods. Phys. Rev. B 54, 5340–5345. (doi:10. 2001 Numerical atomic orbitals for linear- org/software/matrixswitch/. 1103/PhysRevB.54.5340) scaling calculations. Phys. Rev. B 64, 235111. 30. MatrixSwitch using DBCSR. See https://e-cam. 41. Mohr S, Eixarch M, Amsler M, Mantsinen MJ, (doi:10.1103/PhysRevB.64.235111) readthedocs.io/en/latest/Electronic-Structure- Genovese L. 2018 Linear scaling DFT calculations 50. Allen MP, Tildesley DJ. 1987 Computer Modules/modules/MatrixSwitchDBCSR/readme. for large tungsten systems using an optimized simulation of liquids. Oxford, UK: Oxford html. local basis. Nucl. Mater. Energy 15, 64–70. University Press. 31. Electronic Structure Library (ESL). See http://esl. (doi:10.1016/j.nme.2018.01.002) 51. Kresse G, Furthmüller J. 1996 Efficiency of ab- cecam.org. 42. Payne MC, Teter MP, Allan DC, Arias TA, initio total energy calculations for metals and 32. omm-bundle. See https://gitlab.com/ Joannopoulos JD. 1992 Iterative minimization semiconductors using a plane-wave basis set. ElectronicStructureLibrary/omm-bundle. techniques for ab initio total-energy Comput. Mater. Sci. 6, 15–50. (doi:10.1016/ 33. Tsuchida E. 2007 Augmented orbital calculations: molecular dynamics and conjugate 0927-0256(96)00008-0) Downloaded from https://royalsocietypublishing.org/ on 26 September 2023 minimization method for linear scaling gradients. Rev. Mod. Phys. 64, 1045–1097. 52. Library for specialized dense and sparse matrix R. Soc. Open Sci. 10: 230063 electronic structure calculations. J. Phys. Soc. (doi:10.1103/RevModPhys.64.1045) operations, and deep learning primitives Jpn. 76, 034708. (doi:10.1143/JPSJ.76.034708) 43. Lin L, Yang C, Meza J, Lu J, Ying L, Weinan E. (LIBXSMM). See https://github.com/hfp/ 34. Bowler DR, Miyazaki T. 2010 Calculations for 2011 SelInv—an algorithm for selected libxsmm, https://libxsmm.readthedocs.io/en/ millions of atoms with density functional inversion of a sparse symmetric matrix. ACM latest/. theory: linear scaling shows its potential. Trans. Math. Softw. 37, 40. (doi:10.1145/ 53. Lebedeva IV, García A, Artacho E, Ordejón P. J. Phys.: Condens. Matter 22, 074207. (doi:10. 1916461.1916464) 2023 Modular implementation of linear and 1088/0953-8984/22/7/074207) 44. Bowler DR, Gillan MJ. 1998 Length-scale ill cubic-scaling orbital minimization methods in 35. Hohenberg P, Kohn W. 1964 Inhomogeneous conditioning in linear-scaling DFT. Comput. Phys. SIESTA. Zenodo (doi:10.5281/zenodo.7781100). electron gas. Phys. Rev. 136, B864–B871. Commun. 112, 103–111. (doi:10.1016/S0010- 54. Lebedeva IV, García A, Artacho E, Ordejón P. (doi:10.1103/PhysRev.136.B864) 4655(98)00061-7) 2023 Extension of libOMM and MatrixSwitch 36. Kohn W, Sham LJ. 1965 Self-consistent 45. Gan CK, Haynes PD, Payne MC. 2001 libraries for modular implementation of linear equations including exchange and correlation Preconditioned conjugate gradient method for and cubic-scaling orbital minimization methods effects. Phys. Rev. 140, A1133–A1138. (doi:10. the sparse generalized eigenvalue problem in in electronic structure codes using atomic 1103/PhysRev.140.A1133) electronic structure calculations. Comput. Phys. orbitals. Zenodo. (doi:10.5281/zenodo.7781174) 37. Artacho E, del Bosch LM. 1991 Nonorthogonal Commun. 134, 33–40. (doi:10.1016/S0010- 55. Lebedeva IV. 2022 Performance of linear and basis sets in quantum mechanics: 4655(00)00188-0) cubic-scaling solvers based on modular representations and second quantization. Phys. 46. Mostofi AA, Haynes PD, Skylaris C-K, Payne MC. implementation of the orbital minimization Rev. A 43, 5770–5777. (doi:10.1103/PhysRevA. 2003 Preconditioned iterative minimization for method in the SIESTA code. Mendeley Data. 43.5770) linear-scaling electronic structure calculations. (doi:10.17632/c8kz58bg5z.1) 38. Goedecker S, Colombo L. 1994 Efficient linear J. Chem. Phys. 119, 8842–8848. (doi:10.1063/1. 56. Lebedeva IV, García A, Artacho E, Ordejón P. scaling algorithm for tight-binding molecular 1613633) 2023 Modular implementation of the linear dynamics. Phys. Rev. Lett. 73, 122–125. (doi:10. 47. Perdew JP, Zunger A. 1981 Self-interaction and cubic-scaling orbital minimization methods 1103/PhysRevLett.73.122) correction to density-functional approximations in electronic structure codes using atomic 39. Goedecker S, Teter M. 1995 Tight-binding for many-electron systems. Phys. Rev. B 23, orbitals. Figshare. (doi:10.6084/m9.figshare.c. electronic-structure calculations and tight- 5048–5079. (doi:10.1103/PhysRevB.23.5048) 6601948)