Efficient computation of the density matrix with error control on distributed computer systems
Anastasia Kruchinina Elias Rudberg Emanuel H. Rubensson


  The recursive polynomial expansion for construction of a density matrix approximation with rigorous error control [J. Chem. Phys. 128, 074106 (2008)] is implemented in the quantum chemistry program Ergo [SoftwareX 7, 107 (2018)] using the Chunks and Tasks matrix library [Parallel Comput. 57, 87 (2016)]. The expansion is based on second-order polynomials and accelerated by the scale-and-fold technique [J. Chem. Theory Comput. 7, 1233 (2011)]. We evaluate the performance of the implementation by computing the density matrix from the Fock matrix in the large-scale self-consistent field calculations. 
  We demonstrate that the amount of communicated data per worker process tends to a constant with increasing system size and number of computer nodes such that the amount of work per worker process is fixed.

Quantumness of Relative Incompatibility 
    Manish Kumar Shukla^1, Rounak Mundra^2, Arun K Pati^1,  Indranil Chakrabarty ^1 ^3, Junde Wu ^4 
    Received May 29, 2019; revised July 5, 2019; accepted
July 10, 2019. Date of publication July 31, 2019 
 date of current version
August 21, 2019
====================================================================================================================================================






§ INTRODUCTION



The density matrix D is constructed from the Fock/Kohn-Sham matrix F in each step of self-consistent field (SCF) calculations emerging in the Hartree–Fock method <cit.> and Kohn–Sham density functional theory <cit.>.
In orthogonal basis set, the matrix F has eigenvalues λ_i and corresponding eigenvectors y_i:

    Fy_i = λ_i y_i.

Each eigenpair (eigenvalue and eigenvector) corresponds to one molecular orbital. 
Let the eigenvalues of the matrix F be arranged
in ascending order

    λ_1≤λ_2≤…≤λ_homo
     <
    λ_lumo≤…≤λ_N-1≤λ_N,

where λ_homo is the eigenvalue corresponding to the
highest occupied molecular orbital (homo eigenvalue), and λ_lumo is the eigenvalue corresponding to the lowest unoccupied molecular orbital (lumo eigenvalue). We assume that the homo-lumo gap ξ := λ_lumo - λ_homo is  non-zero. The number of occupied orbitals we denote with .
The subspace spanned by the eigenvectors of the matrix F, that correspond to the occupied molecular orbitals, is referred to as the occupied
  invariant subspace of F.
At zero temperature the density matrix is the matrix for orthogonal
  projection onto the occupied invariant subspace of F:

    D = ∑_i=1^ y_iy_i^T.

Computation of the density matrix using (<ref>)  requires diagonalization of the matrix F which scales cubically with the matrix size, limiting the calculations to rather small systems. Alternative density matrix methods see the density matrix construction as an evaluation of a matrix function

    D    = θ(μ I - F),
    θ(x)    = 
    
    0     if  x < 0,
    
    1     otherwise,

where μ is located in the homo-lumo gap. 
Recursive polynomial expansions, also known as purification methods, utilize low-order polynomials f_i, i=0,1,2,… to approximate the shifted and scaled step function θ(μ-x):

    D ≈ X_n = f_n(f_n-1(… f_0(F)…)).

Recursive polynomial expansions are implemented in a number of electronic structure codes including CP2K <cit.>, Ergo <cit.>, Honpas <cit.>, and LATTE <cit.>. 
The initial polynomial f_0 usually maps the eigenvalues of F into the interval [0, 1] in reverse order.
The choice of low-order polynomials f_i, i≥1 depends on the   method, but typically one or two matrix multiplications per iteration are required. Niklasson <cit.> proposed an efficient recursive expansion algorithm based on the second-order spectral projection polynomials x^2 and 2x-x^2. We will refer to this algorithm as the SP2 algorithm. In the original SP2 algorithm, the polynomial in each iteration is chosen based on the matrix trace, which may result in slow convergence in some cases. The convergence can be improved by selecting polynomials based on the homo and lumo eigenvalues <cit.>. This will also remove the explicit dependence on the system size. In the scheme proposed in <cit.>, homo and lumo eigenvalue estimates are obtained at negligible cost using Frobenius norms and traces of intermediate matrices X_i-X^2_i. The homo and lumo estimates for the current matrix F are obtained at the end of the recursive expansion, but they can be propagated to the following SCF cycle and used there in the recursive expansion. In the first SCF cycles, where the eigenvalue bounds are not available yet or not accurate enough, we use the trace-correcting recursive expansion mentioned above <cit.>. Knowledge of homo and lumo eigenvalue estimates allows for further improvement of the SP2 scheme convergence by making use of the scale-and-fold technique proposed in <cit.>. The SP2 scheme, accelerated using the scale-and-fold technique, we refer to as the SP2ACC recursive expansion. We use the parameterless stopping criteria for recursive expansions developed in <cit.>, which automatically stop the iterations when the level of numerical errors is reached and the final density matrix approximation cannot be substantially improved.


Matrices, obtained in electronic structure calculations,  posses the exponential decay property of matrix element magnitude <cit.>. After removal of small matrix elements (truncation), matrices are sparse, and the matrix sparsity can be preserved throughout the recursive expansion iterations. Another consequence of the exponential decay property is that non-negligible matrix elements tend to be localized in certain regions. We implement the truncation scheme proposed in <cit.>, allowing for rigorous control of the error in the occupied invariant subspace, which is a measure of the density matrix approximation error. The error is controlled in the spectral norm. However, we use a computationally cheaper upper bound of the spectral norm for selecting truncation threshold values in each iteration. For that purpose, we propose the Frobenius-infinity mixed matrix norm, which provides an upper bound of the spectral norm for Hermitian matrices, and which has the same asymptotic behavior as the spectral norm for matrices with localization properties of element magnitudes. 



The main computational cost in the recursive expansion comes from sparse matrix multiplications, which makes it attractive for parallel computations. Parallel algorithms, developed for dense matrix-matrix multiplication, take advantage of a uniform distribution of matrix elements over computational resources. A few efficient algorithms for sparse matrix-matrix multiplication use similar ideas. A uniform distribution of non-zero matrix elements can be obtained dynamically during runtime by a random permutation of matrix rows and columns (see e. g. <cit.>). The permuted matrices are divided into blocks with a sparse matrix representation of elements at the block level. Then, some algorithm, inspired by the parallel dense matrix-matrix multiplication, is used for managing the communication between computer nodes. Individual blocks are multiplied locally using a chosen algorithm for sparse matrix-matrix multiplication. The communication cost per node of the 2D algorithm, utilizing a two-dimensional process grid <cit.>, scales as (√(N)) with the matrix size N in a weak scaling regime. The parallel performance can be improved by introducing a third process grid dimension <cit.>. Compared to the 2D approach, in a weak scaling limit the communication cost per node of the 3D algorithm is constant <cit.>.

Methods involving random permutation were originally developed for the general sparse matrix-matrix multiplication.  By randomly permuting matrix rows and columns, any locality information is lost. Furthermore, permutations of matrix rows and columns may be needed for some intermediate matrices in the recursive expansion to enforce load balance. An efficient sparse matrix-matrix multiplication, utilizing the locality of non-zero elements, has been implemented in the Chunks and Tasks matrix library (CHTML) <cit.>. It has been shown that the dynamic work and data distribution over computational resources allows for efficient utilization of data locality and constant communication cost per worker process in a weak scaling regime <cit.>. 



In this work, we expand CHTML with matrix operations required in the recursive expansion, and implement the whole recursive expansion using CHTML in the quantum chemistry code Ergo <cit.>. We evaluate the weak scaling performance of our implementation using the Frobenius-infinity mixed norm for the estimation of the spectral norm in the error control scheme. We request the same accuracy of the density matrix approximation for all system sizes.
The usage of the Frobenius-infinity mixed norm allows us to keep the number of non-zero elements per row asymptotically constant with increasing system size. We show numerically that the communication cost per node is constant in a weak scaling limit. 











§ ERROR CONTROL


A common approach to preserve matrix sparsity is to remove blocks of elements with norms smaller than some predefined threshold value. Often such threshold values are provided by the user and not explicitly related to the desired accuracy of the density matrix approximation.

 In <cit.>, an upper bound of the forward error D-X_n is given by

    D-X_n≤D-P__n + P__n-X_n,

where · is any unitary-invariant matrix norm, P__n is the matrix for orthogonal projection onto the occupied subspace _n of the matrix X_n.
The first norm D-P__n measures erroneous rotations of the occupied invariant subspace or errors in eigenvectors. A scheme for keeping the error in the occupied invariant subspace smaller than a desired accuracy γ is proposed in <cit.>. Let X_i = X_i+E_i be the matrix obtained after truncation of X_i,
where the truncation is written as an explicit perturbation of X_i by the error matrix E_i. Then, it has been shown in <cit.>, that

    D-P__n≤∑_i=0^nE_i/ξ̃_i - E_i_2,

where ξ̃_i is a lower bound of the homo-lumo gap in iteration  i of the recursive expansion. Assume that an upper bound n_max of the total number of iterations in the recursive expansion is known. Then, if the truncation threshold values τ_i are chosen such that:

    τ_i := γ/n_max+1ξ̃_i/1+γ/n_max+1,

and the truncation is performed so that 

    E_i≤τ_i,

then the total error in the occupied subspace is bounded by γ:

    D-P__n≤γ.

The error in the occupied subspace is increasing with the number of iterations, and any lost accuracy cannot be recovered in later iterations of the recursive expansion.


As mentioned above, any unitary-invariant norm can be used in the error control scheme. Examples of unitary-invariant norms are the spectral and Frobenius norms.  
The spectral ·_2 and Frobenius ·_F norms of a Hermitian matrix A = (a_ij) of size N are defined as:

    A_2    = ρ(A),
    A_F    = √(∑_i,j=0^Na_ij^2),

where the spectral radius ρ(A) is equal to the largest by absolute value eigenvalue of A.
The spectral norm depends only on the eigenvalue distribution, but it is expensive to compute for large systems. However, the computationally cheaper Frobenius norm depends on the system size, making its usage inappropriate in large scale calculations. Thus, we use the spectral norm in the error control scheme.  If we can measure some upper bound Ŝ_i to the spectral norm of E_i, inequality (<ref>) can be satisfied by making sure that 
Ŝ_i ≤τ_i.




In the presence of numerical errors the computed density matrix approximation slightly deviates from idempotency. The second norm P__n-X_n in (<ref>) measures errors in the density matrix approximation coming from errors in eigenvalues. 
In <cit.>, we proposed a parameterless stopping criteria for the recursive expansions, which do not require any user defined tolerance. The criteria are satisfied when the numerical errors in the density matrix approximation, coming from the removal of small matrix elements and round-off errors, prevent any further improvement of the result. The idea is to observe the numerically computed convergence order during the recursive expansion and automatically detect the convergence order drop. The convergence order in iteration i is  controlled using  the idempotency error X_i-X^2_i, which can be considered as a measure of the error P__i-X_i. Note that the error P__n-X_n may be decreased to a required level of accuracy by reducing truncation in the last recursive expansion iterations. However, due to the increased computational cost, we do not use this possibility here, and choose truncation values given by (<ref>) throughout the whole recursive expansion. If the SP2ACC scheme is used, we fall back to the SP2 scheme in the final iterations when acceleration does not have a strong effect on the convergence, and use the stopping criterion developed for the SP2 scheme, as suggested in <cit.>. In <cit.>, the proposed stopping criteria use the spectral norm. In addition, it is shown numerically, that the Frobenius norm can be used instead of the spectral norm. In contrast to the error control scheme above, the dependence of the Frobenius norm on the system size is less critical in the stopping criterion, since the Frobenius norm provides a good estimate to the spectral norm close to convergence.



In the following section we consider matrix norms providing upper bounds to the spectral norm. Such norms can be used to ensure (<ref>).











§ UPPER BOUNDS FOR THE SPECTRAL NORM




In <cit.>, the Frobenius-spectral mixed norm is proposed. Let a matrix A be partitioned into submatrices:

    A = 
      [ A_11 A_12; A_21 A_22 ].

Then the Frobenius-spectral mixed norm of the matrix A is defined as:

    A_m = 
      [ A_11_F A_12_F; A_21_F A_22_F ]_2.

The Frobenius-spectral mixed norm gives an upper bound of the spectral norm <cit.>:

    A_2≤A_m.

Moreover, the Frobenius-spectral mixed norm has the same asymptotic behavior as the spectral norm, and is computationally cheaper than the spectral norm <cit.>.  

In this work, we introduce a Frobenius-infinity mixed norm, which we refer to as Frob-Inf and denote ·_M, obtained by the combination of the Frobenius and infinity matrix norms.  The infinity norm of a matrix A = (a_ij) of size N is equal to the largest sum of absolute values of elements in each matrix row:

    A_∞ = max_1≤ i ≤ N∑_j=1^N|a_ij|.

If the matrix A is divided into submatrices as in (<ref>), the Frob-Inf norm of the matrix A is defined as:

    A_M = 
      [ A_11_F A_12_F; A_21_F A_22_F ]_∞.

In practice, with increasing matrix size the submatrix size is kept fixed, and the number of submatrices in Frobenius-spectral and Frob-Inf norms is increasing. Thus, the Frobenius-spectral mixed norm requires computation of the spectral norm of matrices of growing size, which may be expensive for very large systems. In addition, efficient implementation of the spectral norm is difficult on distributed systems.


Under the assumption that the matrix A is Hermitian, the infinity norm gives an upper bound of the spectral norm <cit.>. Thus, we have

    A_2≤A_m≤A_M,

which proves that the Frob-Inf norm is an upper bound of the spectral norm. We note that the infinity norm can be used for estimating the spectral norm in the error control scheme. However, due to the simplicity of the implementation on distributed systems we chose to use the Frob-Inf norm instead.


We perform numerical tests in Matlab and compare values of matrix norms. In Figure <ref> we show results for a symmetric banded matrix with random elements drawn from the standard uniform distribution on the open interval ]0,   10^-3[ and bandwidth 2× 300 + 1.
The Frob-Inf norm, computed with block sizes 100 and 500, is much smaller than the Frobenius norm and has the same asymptotic behavior as the spectral norm. However, symmetric random permutation of rows and columns delays the moment when the Frob-Inf norm becomes proportional to the spectral norm compared to the non-permuted case. Note that in Figure <ref> we use much smaller block sizes than in Figure <ref>. 



It is hard to predict the behavior of the Frob-Inf norm for general matrices. However, if matrices posses the locality property and number of non-zeros per row does not grow significantly with matrix size, we can expect that, for large enough matrix sizes, the Frob-Inf norm will be proportional to the spectral norm, and our practical experiments support it.






§ DENSITY MATRIX CONSTRUCTION



In this section we describe the density matrix construction using the SP2ACC recursive expansion implemented in the quantum chemistry program Ergo <cit.> and parallelized using CHTML. Assuming that homo and lumo estimates are available, the density matrix algorithm based on the SP2ACC expansion may be divided into three steps:


  Step 1 Initialize:
     
     
  * Following Algorithm 4 in <cit.>, estimate the number of recursive expansion iterations n_max, determine the acceleration parameters α_i and the polynomial sequence p_i, i=1, 2,…, n_max, where p_i = 1 if f_i(x) = ((1-α_i) + α_ix)^2 and p_i = 0 if f_i(x) = 2α_ix-(α_ix)^2. Moreover, select the iteration n_min
      when acceleration should be switched off, i. e. α_i = 1 for i≥ n_min. 
     
  * Select truncation threshold values τ_i, i=0, 1,…, n_max defined by (<ref>).
     
  * Rescale the matrix F such that all its eigenvalues lie in the interval [0, 1] in reverse order, see Algorithm <ref>. Here λ_max≥λ_N and λ_min≤λ_1 are obtained using Gershgorin circles theorem. 
     
     
  
  Step 2 Recursively apply polynomials p_i until the stopping criterion is satisfied, see Algorithm <ref>. The constant C in line <ref> of Algorithm <ref> is used in the stopping criterion for the SP2 recursive expansion, and it is derived in <cit.>.
  
  
  
  Step 3 
  
    Compute homo and lumo eigenvalue estimates using Frobenius norms e_i and traces t_i following Algorithm 3 in <cit.>. These estimates are then propagated to the next SCF cycle.
  


In Algorithms <ref> and  <ref> the truncation is written as an explicit addition of some error matrix E_i in each iteration. In this work, we use truncation based on the Frob-Inf and Frobenius norms. The SP2 recursive expansion can be obtained using the algorithm outlined above with the acceleration parameters α_i set to 1 for all i≥ 1, and n_min = 1.

















§ IMPLEMENTATION DETAILS



Chunks and Tasks <cit.> is a task-based programming model, where the runtime system is responsible for scheduling tasks and distributing the data into the physical resources. The idea is that the user splits the work into smaller tasks and the data into smaller chunks, and registers both to the library. As soon as chunks are registered to the library, they are read-only. Upon registration, chunks and tasks are assigned a unique identifier, which is chosen by the library. Since each task has input chunks which has to be registered to the library, it will never happen that some task requests non-existent chunks. Moreover, information about the location of the chunk can be stored in the chunk identifier, allowing efficient and easy data fetching. The CHT-MPI library <cit.> implements the Chunks and Tasks programming model. The task scheduler is based on randomized work stealing. As soon as one process runs out of work, it attempts to steal work from another process chosen at random. Tasks are scheduled for execution as soon as all input chunks are available, and each task can register new tasks and chunks.  



Chunk objects may store different kinds of data, in particular they can
store chunk identifiers of other chunks, giving the possibility to create
hierarchical structures.
The Chunks and Tasks Matrix Library (CHTML) is a hierarchical matrix library written using the Chunks and Tasks programming model and implemented on top of the CHT-MPI library. In each level of the hierarchy, the matrix is divided into four blocks forming a quad-tree structure. The tree nodes on the lowest level of the hierarchy are called leaf nodes. Each non-leaf node stores chunk identifiers of corresponding non-zero submatrices. If the submatrix is zero, its chunk identifier is set to  CHUNK_ID_NULL. On the leaf level we use the block-sparse matrix representation where the leaf matrix is divided into small blocks and only elements of non-zero blocks are stored. Other representations of the leaf matrix are possible.



We implemented a set of classes, wrapping the CHTML matrix type and various matrix operations. As mentioned above, when a chunk is registered to the library, the user in return receives  a chunk identifier which can be used as an input to tasks. Inside the wrapper class we do not store any matrix elements, but only the identifier of a root chunk of a quad-tree representing the matrix.   The actual data may be distributed on different computational nodes by some previous calculation of the matrix. We provide alternative implementations of tasks for general and symmetric matrices. The symmetric matrix wrapper implementation utilizes matrix symmetry and assumes that only the upper triangle of the matrix is stored in memory.






§ NUMERICAL EXPERIMENTS



In this section, we verify the efficiency of the density matrix construction implemented using CHTML in Ergo <cit.>. We run our tests on the Rackham cluster at the UPPMAX computer center at Uppsala university, comprised of 486 nodes with two 10-core Intel Xeon V4 CPU's running at 2.2 GHz/core, each with at least 128 GB of RAM memory. We use the Chunks and Tasks matrix library linked to CHT-MPI 1.2 <cit.> compiled with GCC 8.2.0 and Open MPI 3.1.3. The performance of leaf matrix operations is improved by linking to the optimized BLAS library OpenBlas 0.2.19 (single threaded).
We set maximum leaf matrix size to 4096 and small block size to 32. In the current implementation, the block size in the Frob-Inf norm is equal to the maximum leaf matrix size, i.e. 4096.  All timing results represent elapsed wall time.

To obtain the homo-lumo eigenvalue estimates we run the recursive expansion in each test twice, and measure time and communication cost for the second run. In the first run we do not have any assumption about the eigenvalues. Thus, we use the trace-correcting recursive expansion scheme <cit.> with ad hoc truncation tolerances. At the end of the recursive expansion we obtain homo-lumo eigenvalue estimates and use them in the second run as described in section <ref>. 





 §.§ Weak scaling tests on block-diagonal matrices


In a weak scaling test the amount of work per worker process is kept fixed. In order to achieve this, the number of non-zero elements per row in all (Fock, intermediate and final density) matrices should be kept constant with increasing system size. Moreover, to get the same number of recursive expansion iterations, the homo-lumo gap should be kept constant too.
We construct block-diagonal matrices of increasing size with the desired properties mentioned above. In the created matrices we repeat the same test matrix along the diagonal. The test matrix is a Fock matrix of size 19204  obtained in the last SCF cycle of Hartree–Fock calculations using the standard Gaussian basis set 3-21G on a glutamic acid and alanine (Glu-Ala) helix with 3330 atoms. The estimated number of SP2ACC recursive expansion iterations is 24 for all cases and the actual number of iterations is 18 for the smallest matrix size and 20 for all other matrices. The maximum allowed error in the occupied subspace is set to 0.001. Truncation is based on the Frob-Inf norm. In Figure <ref> we show that the number of non-zero elements per row is constant for large enough matrices.



In our tests we increase the number of blocks K along the diagonal proportionally to the number of worker processes N_p such that N_p = 2K. Since we are interested in the communication cost per worker process and not in the actual timings, we put one worker process per core and use only one worker thread per worker. We do not use the chunk cache, it is set to 0 GB. 
In Figure <ref> we demonstrate that the amount of data received per worker process tends to a constant with increasing number of workers. We note that these results are in agreement with the results obtained in <cit.>, where only the sparse matrix-matrix multiplication is analyzed.

We note that this is an illustrative example. The ideal distribution of work and data over processes leads to separate computation of each diagonal block, without any communication of matrix elements between blocks. However, if the random permutation is used to achieve load balance, such trivial non-zero matrix structure cannot be exploited, resulting in unnecessary communication of matrix elements. In our approach, the load balance is achieved dynamically, and we are able to utilize the non-zero matrix structure reducing the communication cost. 









 §.§ Calculations on water clusters


In this section we present results of density matrix calculations for water clusters. Water cluster geometries are generated from a large molecular dynamics simulation of bulk water at standard temperature and pressure by including all water molecules within spheres of varying radii. We perform 
calculations using our test implementation of the Hartree–Fock method parallelized for distributed-memory. We use the Gaussian basis set 3-21G. Initial guess density matrices are obtained from calculations with a smaller basis set STO-2G. We perform 5 SCF cycles and save Fock matrices in the orthogonal basis in the last cycle. The recursive expansion is applied on the saved Fock matrices. 
The chunk cache size is set to 8 GB per worker. Two workers are placed on each node (using  Open MPI flags) and each worker has 9 worker threads. 
The maximum allowed error in the occupied subspace is set to 0.001. We note that the homo-lumo gap does not change significantly with increasing problem size. The estimated number of iterations in SP2 and SP2ACC recursive expansions is 35-37 and 22-23, respectively. The obtained number of recursive expansion iterations in SP2 and SP2ACC recursive expansions is 29 and 16-17, respectively.

In our tests the number of nodes is increasing proportionally to the matrix size. Number of non-zeros per row in the Fock matrices, obtained density matrices and the maximum number of non-zeros per row between all intermediate matrices X_i are presented in Figure <ref> for both SP2 and SP2ACC recursive expansions. 





In Figure <ref> we compare the execution time and the communication cost of calculations where truncation is based on the  Frobenius and Frob-Inf norms. We use the wrapper class representing a general matrix.



In <cit.>, it is shown that the execution time of sparse matrix-matrix multiplication implemented in CHTML scales as ((log(N_p))^2) where N_p is a number of worker processes.
In Figure <ref>, data indicate that the same result is obtained for the whole recursive expansion. 




We perform recursive expansions using the CHTML wrapper for general matrices where matrix symmetry is ignored, and the wrapper for symmetric matrices where only the upper triangle of the matrix is stored and manipulated. Truncation based on the Frobenius norm is used in both cases. The total execution time and the amount of data received per worker process are given in Figure <ref>. The figure shows that CHTML with the block-sparse leaf matrix is able to utilize matrix symmetry, providing a reduction of the total execution time and communication cost by almost a factor of two. 










§ CONCLUDING REMARKS


In this work, the density matrix is constructed using the SP2ACC recursive polynomial expansion, which scales linearly with system size for sufficiently sparse matrices. The matrix sparsity is enforced by removal of small matrix elements. The combination of the error control scheme <cit.> and parameterless stopping criteria <cit.> enables us to use only one user-defined tolerance for the density matrix computation: the maximum allowed error in the occupied invariant subspace. Truncation tolerances are  chosen dynamically at the beginning of the recursive expansion to ensure the desired accuracy of the density matrix approximation.  


We study the execution time and communication cost of the recursive polynomial expansion implemented using the Chunks and Tasks matrix library (CHTML) <cit.>. Matrices are expressed as quad-trees of chunk objects, where the matrix elements are stored on the lowest level using a block-sparse matrix format. The parallelism is then expressed using tasks operating on the matrix hierarchies. The sparse matrix-matrix multiplication is an important building block of the recursive expansion. Here we use the locality-aware sparse matrix-matrix multiplication implemented in CHTML <cit.>. 
We demonstrate that, if the number of recursive expansion iterations and the number of non-zero elements per row in all matrices is bounded with increase of the system size, the communication cost is asymptotically constant. Moreover, we show that the quad-tree matrix representation in CHTML allows for efficient symmetry utilization.

To the best of our knowledge this is the first implementation of the density matrix construction with rigorous error control of the density matrix approximation and bounded communication cost in a weak scaling limit. The current study focuses on the computation of the density matrix from a given Fock/Kohn-Sham matrix. Future work includes parallelization using Chunks and Tasks of the SCF procedure implemented in the Ergo code <cit.>.




§ ACKNOWLEDGMENTS


This work was supported by the Swedish strategic research programme eSSENCE.
Computational resources were provided by the Swedish National
Infrastructure for Computing (SNIC) through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) in Uppsala, Sweden.



 siam 

