Linear scale electronic structure calculation method and system based on dyeing superposition state and terminal

Abstract

The application provides a linear scale electronic structure calculation method and system based on a dyeing superposition state and a terminal, and is based on the dyeing superposition state. The method is used for solving the sparse operator or matrix function with spatial locality or limited correlation length, which is involved in the Kohn-Sham self-consistent process under the local basis representation, such as the density matrix, the inverse square root of the overlap matrix, and the operator required by the orthogonal representation transformation. Under the premise of maintaining the linear scale complexity, the application breaks through the dependence of the traditional algorithm on the sparse matrix-sparse matrix multiplication (SpMSpM), and reconstructs the core calculation into the regular sparse matrix-dense matrix multiplication (SpMM). Through the transformation of the bottom layer calculation paradigm, the regularity of the Kohn-Sham electronic structure self-consistent calculation and the parallel adaptation ability of the bottom layer hardware are significantly improved, and the inherent bottleneck of the large pre-factor and the low hardware execution efficiency caused by the limitation of the core operator (SpMSpM) in the traditional linear scale algorithm is fundamentally relieved.

Description

Technical Field

[0001] This invention relates to the field of first-principles electronic structure calculation, and in particular to a linear scaling electronic structure calculation method, system, and terminal based on colored superposition states. Background Technology

[0002] First-principles calculations are methods for predicting material properties based on quantum mechanics principles and without relying on empirical parameters. Among these, Kohn-Sham density functional theory (KSDFT) is the most widely used framework, successfully transforming the multi-electron problem into an independent-particle problem within a single-electron effective potential field. Its core is solving the Kohn-Sham equations: ;in, Represents the Kohn-Sham Hamiltonian. and These represent the single-electron orbital wavefunction and its corresponding energy. The traditional solution obtains all eigenstates by diagonalizing the Hamiltonian matrix, and then constructs the density matrix using the Fermi-Dirac distribution. ;in, Let Fermi-Dirac distribution function be used. The density matrix of the system can be expressed as a matrix function of the Hamiltonian matrix. Solving for the density matrix is ​​a crucial step in obtaining the total energy, charge distribution, and atomic forces of the system. Under the non-orthogonal local basis set representation, to perform Löwdin orthogonalization and construct the equivalent Hamiltonian under the orthogonal representation, it is also necessary to solve for the inverse square root of the overlap matrix. From a mathematical structure perspective, the density matrix, The calculations of other similar physical quantities in electronic structure can be reduced to solving a class of large-scale matrix functions. However, the computational complexity of traditional diagonalization methods increases cubically. As the scale of the system increases, it will quickly face the "curse of dimensionality", which will limit its application to systems of hundreds to thousands of atoms, making it difficult to meet the simulation needs of large-scale complex systems at the million-atom level in the fields of materials science and life science.

[0003] To overcome the aforementioned scaling bottlenecks, in recent years, the academic community has developed a series of linear scaling (O(N)) algorithms based on the exponential decay locality of the density matrix in insulators and high-temperature metal systems. Among them, the Second-Order Spectral Projection (SP2) method, which recursively and rapidly approximates the Fermi-Dirac density matrix, transforms the construction of the density matrix into an iterative process of a series of sparse matrix-sparse matrix multiplications (SpMSpM), and is one of the representative schemes for achieving linear scaling.

[0004] However, existing linear scaling techniques still face bottlenecks in practical engineering applications. First, they suffer from poor hardware adaptability. Many existing linear scaling methods heavily rely on sparse matrix-sparse matrix multiplication (SpMSpM). Such computations typically exhibit irregular memory access characteristics and scattered data access, making it difficult to fully utilize the cache levels, prefetching mechanisms, vectorized instructions, and efficient block computation capabilities of modern processor architectures (such as multi-core CPUs and GPUs). This results in low hardware utilization and large pre-factors in actual operation, making it difficult to translate "theoretical linear scaling" into efficient practical computation. Second, a general solution framework is lacking. For matrix functions defined by different scalar functions, existing methods often require different iterative processes or expansion strategies, lacking versatility and hindering the formation of a unified matrix function solution framework. Third, in large-scale systems, although the target matrix often exhibits spatial locality, sparsity, or near-sparseness, existing methods do not fully utilize these structural characteristics, particularly in balancing local compression representation, efficiency in solving local subproblems, and overall reconstruction accuracy, leaving room for further improvement.

[0005] Therefore, the current field of self-consistent computation of large-scale electronic structures requires a novel method for solving matrix functions. By employing more efficient data compression methods and local subspace solution strategies, the dependence on irregular sparse matrix multiplication can be reduced, and the hardware affinity of the computation process to modern high-performance parallel computing platforms can be improved. This will enable efficient and high-precision solutions to various large-scale matrix functions (such as density matrices, inverse square root matrices, etc.) within a unified computational framework. Summary of the Invention

[0006] Given that existing technologies for large-scale self-consistent computation of the electronic structure of the Kohn-Sham equation have technical problems such as irregular sparse matrix operations, low hardware execution efficiency, and large pre-factors in the process of solving sparse operators or matrix functions,

[0007] To achieve the above and other related objectives, this invention provides a linear scaling electronic structure calculation method based on colored superposition states. The method includes: coloring and grouping basis function indices under a local basis set representation based on the spatial locality or finite correlation length of a target sparse operator or matrix function; constructing colored superposition states as a low-dimensional subspace basis based on the grouped basis function indices; projecting the target operator or matrix function onto the low-dimensional subspace formed by the colored superposition states to form a compressed low-dimensional subspace representation; calculating the compression result of the target operator or matrix function acting on the colored superposition state matrix under the compressed low-dimensional subspace representation; extracting the non-zero matrix elements of the target operator or matrix function under the original local basis set representation from the compressed low-dimensional subspace representation and reconstructing them into a sparse matrix form of the target operator or matrix function for self-consistent calculation of the electronic structure of the Kohn-Sham equation.

[0008] In one embodiment, the target operator or matrix function includes one or more of the following: density matrix, inverse square root of overlap matrix, and operator-related matrix function required for orthogonal representation transformation.

[0009] In one embodiment of the present invention, the step of using a graph coloring strategy to color and group the basis function indices under the local basis set representation, and constructing a colored superposition state as a low-dimensional subspace basis based on the grouped basis function indices includes: coloring and grouping the basis function indices with low correlation or no correlation; superimposing the local basis functions corresponding to the basis function indices of the same color to construct a colored superposition state of the corresponding color; forming a colored superposition state matrix from multiple colored superposition states corresponding to multiple colors, and using the colored superposition state matrix as a low-dimensional subspace basis.

[0010] In one embodiment of the present invention, the graph coloring strategy includes: constructing an adjacency graph based on the truncation radius, spatial locality, finite association length, or non-zero matrix element structure of the target sparse operator or matrix function, using local basis function indices as vertices; coloring the adjacency graph such that there are no adjacent edges between local basis functions corresponding to vertices of the same color, or such that the distance between local basis functions corresponding to vertices of the same color is greater than a preset truncation radius; grouping the local basis functions based on different colors and constructing a color superposition state for each color, forming a color superposition state matrix of dimension N×Nc as a low-dimensional subspace basis, and denoting the color superposition state matrix as... Where N is the total number of local basis functions and Nc is the number of color types.

[0011] In one embodiment of the present invention, the method of constructing an adjacency graph includes: reading the atomic structure information and local basis function information of the target system; for all target sparse operators or matrix functions, the nearest neighbor atoms within the truncation radius range are processed, and a nearest neighbor atom table is established for each atom to record its nearest neighbor information; an adjacency graph is constructed with local basis function indices as vertices, and adjacency edges are established between corresponding local basis function indices based on the nearest neighbor relationship.

[0012] In one embodiment of the present invention, the low-dimensional subspace compression representation includes the compression result obtained by applying a target sparse operator or matrix function to a colored superposition state matrix. If the target operator or matrix function is denoted as F, then the low-dimensional subspace compression representation includes the compression result obtained by applying F to a colored superposition state matrix. The first dense matrix represented by and by The second dense matrix is ​​represented by the first dense matrix, and F is formed by multiplying the first dense matrix and the second dense matrix. The compressed decomposition representation in this form.

[0013] In one embodiment of the present invention, the method of calculating the compression result of the target sparse operator or matrix function acting on the colored superposition matrix adopts the block Lanczos-Krylov subspace method. The method includes: using the colored superposition matrix as multiple right-hand vectors, mapping the calculation process of the target sparse operator or matrix function acting on the colored superposition matrix to the Krylov subspace through sparse matrix-dense matrix multiplication and block Lanczos iteration between the matrix to be processed and the multiple right-hand vectors; calculating the target matrix function in the Krylov subspace; and mapping the calculation result to the compression result of the target sparse operator or matrix function acting on the colored superposition matrix.

[0014] In one embodiment of the present invention, the step of extracting the non-zero matrix elements of the target sparse operator or matrix function under the original local basis set representation from the low-dimensional subspace compressed representation and reconstructing it into a sparse matrix form of the target sparse operator or matrix function includes: determining the matrix element positions that satisfy the space truncation condition under the original local basis set representation of the target operator or matrix function based on the nearest neighbor atom table, wherein the matrix element positions include corresponding row indices and column indices; establishing a mapping relationship between the matrix element positions under the original local basis set representation and the corresponding element positions in the low-dimensional subspace compressed representation based on the color superposition state; extracting the non-zero matrix elements of the target operator or matrix function from the low-dimensional subspace compressed representation according to the mapping relationship; and writing the extracted non-zero matrix elements into a sparse matrix data structure according to their row indices and column indices under the original local basis set representation, thereby reconstructing the sparse matrix form of the target operator or matrix function under the original local basis set representation.

[0015] To achieve the above objectives, this invention also provides a linear scaling electronic structure calculation system based on colored superposition states. The system includes a colored superposition state construction module, used to color and group the basis function indices under the local basis set representation using a graph coloring strategy based on the atomic structure information, local basis function information, and spatial locality or finite correlation length of the target sparse operator or matrix function, and construct colored superposition states as a low-dimensional subspace basis based on the grouped basis function indices; a compressed representation generation module, used to project the target sparse operator or matrix function onto the low-dimensional subspace composed of the colored superposition states to form a low-dimensional subspace compressed representation of the target sparse operator or matrix function; a compressed representation solving module, used to calculate the compressed result of the target sparse operator or matrix function acting on the colored superposition state matrix under the low-dimensional subspace compressed representation; and a non-zero element extraction module, used to extract the non-zero matrix elements of the target sparse operator or matrix function under the original local basis set representation from the low-dimensional subspace compressed representation, and reconstruct them into a sparse matrix form of the target sparse operator or matrix function for use in Kohn-Sham... Self-consistent calculation of electronic structure.

[0016] To achieve the above and other related objectives, the present invention provides an electronic terminal, comprising: one or more memories and one or more processors; the one or more memories being used to store a computer program; and the one or more processors being connected to the memories and used to run the computer program to execute the method.

[0017] As described above, this invention is a method, system, and terminal for calculating linear scaling electronic structures based on stained superposition states. Compared with existing technologies, this invention has at least the following advantages: This invention uses stained superposition states to represent the sparse operators or matrix functions that need to be solved in the self-consistent calculation of Kohn-Sham electronic structures as a low-dimensional subspace compressed form. Under this compressed representation, the compressed result of the target operator or matrix function acting on the stained superposition state matrix is ​​calculated. Then, based on the stained mapping relationship, the non-zero matrix elements under the original local basis set representation are extracted, and the sparse matrix form is reconstructed. Therefore, the irregular sparse matrix-sparse matrix multiplications relied upon in traditional linear scaling electronic structure calculations can be reduced or replaced. Compared with the traditional linear scaling electronic structure calculation method that relies on irregular sparse matrix-sparse matrix multiplications, this invention can transform the relevant calculations into sparse matrix-dense matrix multiplications, improving the regularity of the calculation process and the hardware parallelism adaptability. While maintaining the complexity of linear scaling, it improves the efficiency of solving for the target operator or matrix function in large-scale Kohn-Sham electronic structure self-consistent calculations. Attached Figure Description

[0018] Figure 1 The diagram shown is a schematic representation of a linear scaling electronic structure calculation method based on dye superposition states according to an embodiment of the present invention.

[0019] Figure 2 (a) shows a schematic diagram of the projection operation of a sparse matrix function in an embodiment of the present invention.

[0020] Figure 2 (b) shows a schematic diagram of the staining results of a one-dimensional periodic atomic chain adjacency diagram according to an embodiment of the present invention.

[0021] Figure 2 (c) shows a schematic diagram of the dyed superposition state corresponding to the dyed one-dimensional periodic atomic chain adjacency diagram in one embodiment of the present invention.

[0022] Figure 3 (a) shows a schematic diagram of the convergence of total energy as the subspace dimension v increases in an embodiment of the present invention.

[0023] Figure 3 (b) Calculation shown in an embodiment of the present invention A diagram showing the efficiency comparison with the SP2 method.

[0024] Figure 3 (c) shows a schematic diagram comparing the efficiency of the method in one embodiment of the present invention with the SP2 method and the diagonalization method.

[0025] Figure 4 (a) shows a schematic diagram of charge and force accuracy in one embodiment of the present invention.

[0026] Figure 4 (b) shows a schematic diagram comparing the accuracy of charge and force in one embodiment of the present invention.

[0027] Figure 4 (c) Shows a schematic diagram comparing the charge and force accuracy of the SP2 method.

[0028] Figure 4 (d) shows a schematic diagram comparing the charge and force accuracy of the SP2 method.

[0029] Figure 5 (a) shows a schematic diagram of the radial distribution function of the MD statistical water molecule system for calculating a 10000 water molecule system in one embodiment of the present invention.

[0030] Figure 5 (b) shows a schematic diagram of the convergence of the total energy of the system with millions of water molecules in one embodiment of the present invention.

[0031] Figure 6 The diagram shown is a schematic representation of a linear scaling electronic structure calculation method based on dye superposition states according to an embodiment of the present invention.

[0032] Figure 7The diagram shown is a structural schematic of an electronic terminal according to an embodiment of the present invention. Detailed Implementation

[0033] The following specific examples illustrate the implementation of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and various details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that, unless otherwise specified, the following embodiments and features described therein can be combined with each other.

[0034] It should be noted that in the following description, reference is made to the accompanying drawings, which illustrate several embodiments of the invention. It should be understood that other embodiments may also be used, and changes in mechanical composition, structure, electrical system, and operation may be made without departing from the spirit and scope of the invention. The following detailed description should not be considered limiting, and the scope of the embodiments of the invention is defined only by the claims of the published patents. The terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the invention. Spatially related terms, such as “upper,” “lower,” “left,” “right,” “below,” “below,” “lower part,” “above,” “upper part,” etc., may be used herein to illustrate the relationship between one element or feature shown in the figures and another element or feature.

[0035] Throughout this specification, when it is said that a part is "connected" to another part, this includes not only "direct connection" but also "indirect connection" by placing other elements in between. Furthermore, when it is said that a part "includes" a certain constituent element, unless otherwise stated otherwise, this does not exclude other constituent elements, but rather means that other constituent elements may also be included.

[0036] The terms "first," "second," and "third," etc., used herein are for the purpose of describing various parts, components, regions, layers, and / or segments, but are not limiting. These terms are used only to distinguish one part, component, region, layer, or segment from others. Therefore, the "first part," "component," "region," "layer," or "segment" described below may refer to a "second part," "component," "region," "layer," or "segment" without departing from the scope of this invention.

[0037] Furthermore, as used herein, the singular forms “a,” “an,” and “the” are intended to include the plural forms as well, unless the context indicates otherwise. It should be further understood that the terms “comprising,” “including,” indicate the presence of the stated feature, operation, element, component, item, kind, and / or group, but do not preclude the presence, occurrence, or addition of one or more other features, operations, elements, components, items, kinds, and / or groups. The terms “or” and “and / or” as used herein are interpreted as inclusive, or mean any one or any combination thereof. Thus, “A, B, or C” or “A, B, and / or C” means “any one of: A; B; C; A and B; A and C; B and C; A, B, and C.” Exceptions to this definition arise only when combinations of elements, functions, or operations are inherently mutually exclusive in some manner.

[0038] This invention provides a linear scaling electronic structure calculation method based on stained superposition states. This method is used to solve for target operators or matrix functions involved in the Kohn-Sham self-consistent calculation process under local basis set representations, such as the density matrix, the inverse square root of the overlap matrix, operators required for orthogonal representation transformations, orthogonalized Hamiltonian operators, and the inverse function of the overlap matrix. These target operators or matrix functions possess spatial locality or finite correlation length under local basis set representations; therefore, a low-dimensional subspace compression representation can be constructed using stained superposition states, and numerical solutions and non-zero matrix element extraction can be performed under this compression representation.

[0039] This invention employs a graph coloring strategy to color and group the low-association or unassociation basis function indices in the operators or matrix functions, and constructs a colored superposition state as a low-dimensional subspace basis based on the same-colored basis function indices, thereby achieving a low-dimensional subspace compressed representation of the original high-dimensional matrix function solution problem. The compressed result of the target operator or matrix function acting on the colored superposition state matrix is ​​then calculated under this low-dimensional subspace compressed representation. Finally, the non-zero matrix elements of the target operator or matrix function are extracted from the compressed representation of the colored superposition state. Through this process, this invention transforms the sparse-sparse matrix multiplication (SpMSpM) in traditional linear scaling computation into a more efficient sparse-dense matrix multiplication (SpMM), achieving linear scaling complexity while improving computational regularity and parallel adaptability to modern computing hardware (such as multi-core CPUs / GPUs), thus forming a hardware-efficient linear scaling solution mode. The following detailed description, with reference to the accompanying drawings, illustrates embodiments of the invention to facilitate implementation by those skilled in the art. This invention can be embodied in various forms and is not limited to the embodiments described herein.

[0040] like Figure 1 The linear scaling electronic structure calculation method based on dye superposition states shown in the embodiments of the present invention includes the following steps.

[0041] Step S1: Based on the atomic structure information and local basis function information of the target system, as well as the spatial locality or finite correlation length of the target sparse operator or matrix function, a graph coloring strategy is used to color and group the basis function indices under the local basis set representation, and a colored superposition state is constructed based on the grouped basis function indices as a low-dimensional subspace basis.

[0042] In one embodiment, step S1 includes:

[0043] Step S11: Use the basis function indices with low or no correlation to color and group the data.

[0044] Step S12: Superimpose the local basis functions corresponding to the basis function indices of the same color to construct the color superposition state of the corresponding color;

[0045] Step S13: A color superposition state matrix is ​​formed by multiple color superposition states corresponding to different colors, and the color superposition state matrix is ​​used as a low-dimensional subspace basis.

[0046] In one embodiment, the graph coloring strategy includes:

[0047] Using local basis function indices as vertices, an adjacency graph is constructed based on the cutoff radius, spatial locality, finite association length, or non-zero matrix element structure of the target sparse operator or matrix function.

[0048] The adjacency graph is colored such that there are no adjacent edges between local basis functions corresponding to vertices of the same color, or the distance between local basis functions corresponding to vertices of the same color is greater than a preset cutoff radius.

[0049] The local basis functions are grouped according to different colors, and a color superposition state is constructed for each color. The resulting color superposition state matrix of dimension N×Nc serves as the basis for the low-dimensional subspace. The color superposition state matrix is ​​denoted as... Where N is the total number of local basis functions and Nc is the number of color types.

[0050] In one embodiment, the adjacency graph is constructed in the following ways:

[0051] Read the atomic structure information and local basis function information of the target system;

[0052] For all target sparse operators or matrix functions, the nearest neighbor atoms within the truncation radius are identified, and a nearest neighbor atom table is created for each atom to record its nearest neighbor information;

[0053] An adjacency graph is constructed using local basis function indices as vertices, and adjacency edges are established between corresponding local basis function indices based on the nearest neighbor relationships.

[0054] Specifically, it includes:

[0055] (1) Setting the cutoff radius of the objective matrix function:

[0056] First, the cutoff radius of the target sparse operator or matrix function needs to be set. Its physical meaning is that when the spatial distance between the atoms to which two local basis functions belong is greater than 1, the spatial distance between them is greater than 1. At that time, the matrix elements between them are negligible within a preset precision. Figure 2 (b) shows that the values ​​of the density matrix elements decay rapidly as the physical distance between two local basis functions increases. The horizontal axis represents the physical distance between local basis functions, and the vertical axis represents the absolute value of the density matrix elements. The value is determined based on the electronic locality characteristics of the target system. In actual calculations, it can be gradually increased for different systems. The convergence test determines the optimal cutoff radius.

[0057] (2) Reading atomic structure information and local basis function information:

[0058] The complete atomic structure information of the target system is read, specifically including: basic structural parameters: atomic coordinates, atom types, and lattice parameters for periodic systems; the local basis functions corresponding to each atom are also obtained. After completing the reading of atomic structure information and local basis function information, the process of constructing the adjacency graph of local basis functions is initiated.

[0059] (3) Based on the periodic boundary conditions, obtain the coordinates of all atoms in the nearest neighbor unit cells;

[0060] For periodic target systems, the unit cell range is expanded based on periodic boundary conditions to obtain the coordinates of all atoms in the nearest neighbor unit cells, ensuring that subsequent nearest neighbor searches cover the cutoff radius. All atoms within the cell (including neighboring atoms not within the unit cell).

[0061] (4) Search for the nearest neighbor atoms within the cutoff radius of all atoms in the unit cell, and build a nearest neighbor atom table for each atom to record its nearest neighbor information;

[0062] For each atom within the unit cell, at the cutoff radius Perform a nearest neighbor search within the range: Traverse all atoms within the unit cell and determine whether the spatial distance between the current atom and other atoms (including atoms extended from the nearest unit cell) is ≤ For each atom, an independent nearest neighbor atom table is established to record information about all nearest neighbor atoms that satisfy the distance condition (including the type, coordinates, and unit cell position of the nearest neighbor atom).

[0063] (5) Using each local basis function index as a vertex, construct a local basis function adjacency graph based on the nearest neighbor atom table.

[0064] Using all local basis functions of the target system as vertices, construct an undirected adjacency graph G=(V,E) based on the nearest neighbor atom table, where:

[0065] Vertex set V: each vertex A unique local basis function (That is, all local basis functions in the system are vertices of the adjacency graph).

[0066] Edge set E: For any two vertices, determine their nearest neighbor atom list based on the nearest neighbor atom list of their respective atoms: If and The two atoms belonging to the same atom are nearest neighbors (spatial distance ≤ If ), then add an undirected edge between the two vertices; if and For different local basis functions belonging to the same atom, an undirected edge is added between the corresponding two vertices by default (because there is inherent coupling between local basis functions of the same atom, which cannot be ignored). Traversing all pairs of local basis functions completes the construction of the edge set E, and finally yields the complete local basis function adjacency graph G=(V,E).

[0067] In one embodiment, coloring the adjacency graph such that there are no adjacent edges between local basis functions corresponding to vertices of the same color, or such that the distance between local basis functions corresponding to vertices of the same color is greater than a preset cutoff radius, includes:

[0068] The graph coloring algorithm is used to color the vertices of the adjacency graph. The core constraint is that the spatial distance between the atoms of the local basis functions corresponding to vertices of the same color is greater than the cutoff radius of the density matrix. Finally, a colored undirected graph that satisfies the constraint is obtained (each vertex corresponding to a local basis function in the graph is assigned a unique color label).

[0069] Specifically, this embodiment uses a greedy coloring algorithm as the preferred implementation, with the following detailed steps: ① Use non-negative integers as color identifiers (e.g., 1, 2, ..., Nc, where Nc is the final total number of colors); ② Traverse all vertices in the adjacency graph in a preset order (e.g., the order of local basis function numbers); ③ For the current vertex to be colored, retrieve the set of colors already occupied by all its adjacent vertices; ④ Assign the current vertex the available color with the smallest number outside this set. The core objective of this coloring process is to ensure that the distance between the local basis functions of vertices of the same color is greater than the distance between the corresponding vertices of the same color. "Under the constraints, minimize the total number of colors Nc as much as possible to improve subsequent compression efficiency. It should be noted that the graph coloring algorithm in this step is not limited to the greedy coloring algorithm; other coloring methods such as backtracking coloring and DSATUR algorithm (saturation-first coloring) can also be used, as long as the local basis function spacing between vertices of the same color is greater than..." The two core requirements are "minimizing the total number of colors Nc".

[0070] After coloring, a colored undirected graph containing Nc colors is obtained, and all N local basis functions in the system are uniquely assigned a certain color. This coloring scheme has significant advantages: ① Dimensionality compression: For typical target systems, the total number of colors Nc is much smaller than the total number of local basis functions N (Nc≪N); ② Scale independence: Nc is determined only by the topological connection structure of the local basis functions of the target system (i.e., the spatial distribution characteristics of the interactions between local basis functions), and is independent of the overall scale of the system (such as the number of atoms and the total number of local basis functions N). As shown in Figure 2(b), a one-dimensional periodic atomic chain is used as an example system (the interatomic spacing is unit length, the local basis set adopts the atomic orbital basis set, each atom contains only 1 s orbital, and the cutoff radius is...). =2, (Representing the spatial distance between any two atoms), after coloring, only 3 colors are needed to satisfy the constraints. Regardless of how the atomic chain length extends (i.e., how much N increases), it can guarantee any nearest neighbor basis set (corresponding to the atomic spacing ≤ 0.5). The vertices of the coloring scheme are not multicolored, which intuitively verifies the scale independence and effectiveness of the coloring scheme.

[0071] In one embodiment, the local basis functions are grouped based on different colors, and a color superposition state is constructed for each color to form a color superposition state matrix of dimension N×Nc as a low-dimensional subspace basis. The color superposition state matrix is ​​denoted as... include:

[0072] Each color in the colored undirected graph is vectorized to construct a set of color state vectors, which are then integrated to form a color state matrix of dimension N×Nc, where N is the total number of local basis functions and Nc is the number of color types.

[0073] The specific construction process is as follows: ① Define the color superposition state vector corresponding to the k-th color. (k=1,2,...,Nc), the vector dimension is consistent with the total number of local basis functions (N dimensions); ② Vector assignment rules: For all local basis function positions corresponding to the k-th color, randomly assign any value between ±1; for the remaining local basis function positions that are not the k-th color, uniformly assign a value of 0. The core purpose of adopting the "random assignment ±1" design is to avoid the coherent superposition of the truncation error of the density matrix in the subsequent calculation process, thereby improving the accuracy of the compressed representation. After vectorizing all colors, a set of color superposition states is obtained. The number of colors is Nc. ③ Color state matrix integration: Arrange all Nc color superposition states sequentially according to their color numbers to construct the color superposition state basis matrix. The matrix has dimensions N×Nc (row dimension corresponds to the total number of local basis functions N, column dimension corresponds to the total number of colors Nc). As shown in Figure 2(c), with... Figure 2 (b) Taking a one-dimensional periodic atomic chain as an example, a schematic diagram of the coloring state vectors corresponding to the three colors is shown, where “1̃” indicates that the position of the local basis function is assigned a random value (±1), and the blank position indicates that the value is assigned to 0.

[0074] Step S2: Project the target sparse operator or matrix function onto the low-dimensional subspace composed of the color superposition states to form a low-dimensional subspace compressed representation of the target sparse operator or matrix function.

[0075] In one embodiment, the target sparse operator or matrix function F is colored with a superposition matrix. By performing a projection transformation, a low-dimensional compressed representation of the target matrix can be obtained.

[0076] The specific form is as follows:

[0077] ;

[0078] This operation decomposes the original large-scale sparse density matrix of dimension N×N into an N×Nc first-dimensional dense matrix. With the Nc×N dimensional second dense matrix The product form.

[0079] Figure 2 (a) illustrates the overall projection operation process of the target matrix function. It should be noted that this projection transformation introduces numerical padding at the zero-element positions of the target matrix function, but the original distribution characteristics and numerical information of the non-zero elements are fully preserved in the compression result. Since the non-zero elements of the target matrix function exist only among the local basis functions corresponding to the nearest-neighbor atoms, and the nearest-neighbor correlations of the local basis functions have been fully recorded through the aforementioned nearest-neighbor atom table, all non-zero element positions of the target matrix function can be accurately indexed based on the nearest-neighbor atom table, thereby achieving high-precision lossless reconstruction of the target matrix function based on the compressed projection result.

[0080] Step S3: Calculate the compression result of the target sparse operator or matrix function acting on the colored superposition matrix under the low-dimensional subspace compression representation.

[0081] After obtaining the low-dimensional subspace basis of the colored superposition state, the compression result of the target sparse operator or matrix function acting on the colored superposition state matrix can be calculated under the compressed representation of the low-dimensional subspace. The calculation method may include one or more of the following: Krylov subspace method, block Lanczos method, Chebyshev polynomial expansion method, density matrix purification method, rational function approximation method, Green's function integration method, or contour integration method. The above calculation methods are used to obtain the low-dimensional subspace compressed representation of the target sparse operator or matrix function acting on the colored superposition state matrix.

[0082] In a preferred embodiment, the calculation of the compression result of the target sparse operator or matrix function acting on the colored superposition matrix under the low-dimensional subspace compression representation adopts the Krylov subspace method of block Lanczos iteration. The method includes: using the colored superposition matrix as multiple right-hand vectors, mapping the calculation process of the target sparse operator or matrix function acting on the colored superposition matrix to the Krylov subspace through sparse matrix-dense matrix multiplication between the matrix to be processed and the multiple right-hand vectors and block Lanczos iteration; calculating the target matrix function in the Krylov subspace; and mapping the calculation result to the compression result of the target sparse operator or matrix function acting on the colored superposition matrix.

[0083] To facilitate the explanation of the basic calculation process of this invention, the following description first uses the solution of the density matrix under the orthogonal local basis set representation as an example. Based on this, the method for solving the inverse square root of the overlap matrix, the operators required for orthogonal representation transformation, and the density matrix under the non-orthogonal local basis set representation will be further explained. Let H be the Hamiltonian of the target system under the orthogonal local basis set representation, and let its density matrix be... This can be expressed as the Fermi-Dirac matrix function of the Hamiltonian: ;

[0084] in Let Fermi-Dirac distribution function be used. Let T be the chemical potential and T be the temperature. Then, the compressed representation of the density matrix acting on the staining superposition matrix can be written as:

[0085] ;

[0086] In one embodiment, the problem of the action of the matrix function can be written as follows using the Green's function integral of the matrix function:

[0087] ;

[0088] Where Z is the integration variable (corresponding to the energy dimension), I is the identity matrix, and Im represents taking the complex imaginary part.

[0089] Based on the core characteristic that all chromatic vectors share the same matrix H to be processed, block Lanczos iterative calculation is used to calculate the product of Hamiltonian and current dense matrix, gradually spanning v-dimensional Krylov subspace of each chromatic vector, and obtaining v-dimensional orthogonal basis vectors for each chromatic vector;

[0090] Since all chromaticity vectors share the same Hamiltonian H, and each chromaticity vector corresponds to a linear system... Since the same Krylov subspace structure is shared throughout the integration domain, the block Lanczos iterative algorithm can be used to simultaneously construct a unified Krylov subspace for all chromatic state vectors.

[0091] The specific iterative process is as follows: by repeatedly calculating the Hamiltonian matrix H and the current dense matrix... The product of these terms gradually spans a Krylov subspace of dimension v. This process involves SpMM operations, which operate on highly regular memory access patterns, fully utilizing CPU cache prefetching mechanisms and vectorized instruction sets for efficient computation. Through a block Lanczos iterative Gram-Schmidt orthogonalization process, each colored state vector is ultimately obtained. Obtain a set of v-dimensional approximately complete orthonormal bases And integrate them into a base matrix. (Dimension N×v, column vectors are the above orthogonal basis).

[0092] Projecting the original system's Hamiltonian H onto the subspace spanned by the orthogonal basis Qi corresponding to each stained state vector yields the v-dimensional subspace Hamiltonian matrix. The projection formula is:

[0093] ;

[0094] in basis matrix The transpose (dimension v×N), projected to obtain It is a dense matrix, and v≪N (the subspace dimension is much smaller than the original system dimension).

[0095] Based on Green's function integral formula and orthogonal basis and subspace Hamiltonian It can be deduced that Representation in subspace:

[0096] ;

[0097] In the formula Hamiltonian for subspace The integral expression of the density matrix using the Green's function.

[0098] Due to the Hamiltonian of the subspace Since the dimension is only v (v≪N), its eigenvalues ​​and eigenvectors can be quickly solved by matrix diagonalization, thus efficiently calculating Nc small-scale matrices. From the density matrix, we can obtain Thus, the original problem of directly solving large-scale sparse matrices of N×N dimensions has been transformed into the diagonalization problem of Nc small-scale dense matrices of v dimensions, achieving a substantial reduction in the dimensionality of computation and significantly reducing computational complexity and hardware resource requirements.

[0099] The core advantage of this process is that by projecting into the subspace, the solution of the density matrix of the N×N dimensional large-scale sparse Hamiltonian in the original problem is transformed into the solution of Nc small-scale dense matrices, which greatly reduces the computational complexity.

[0100] Step S4: Extract non-zero matrix elements from the low-dimensional subspace compressed representation and reconstruct the sparse matrix form of the target sparse operator or matrix function.

[0101] In one embodiment, step S4 includes:

[0102] Based on the nearest neighbor atom table, determine the matrix element positions that satisfy the spatial truncation condition under the original local basis set representation of the target operator or matrix function, wherein the matrix element positions include the corresponding row index and column index;

[0103] Based on the color superposition state, establish the mapping relationship between the matrix element positions under the original local basis set representation and the corresponding element positions in the low-dimensional subspace compressed representation;

[0104] Based on the mapping relationship, extract the non-zero matrix elements of the target operator or matrix function from the low-dimensional subspace compressed representation;

[0105] The extracted non-zero matrix elements are written into the sparse matrix data structure according to their row and column indices under the original local basis set representation, thereby reconstructing the sparse matrix form of the target operator or matrix function under the original local basis set representation.

[0106] Specifically, the matrix element positions that satisfy the spatial truncation condition under the original local basis set representation are determined based on the nearest neighbor atom table. The matrix element positions include the corresponding row index and column index.

[0107] Based on the color superposition state, establish the mapping relationship between the matrix element positions under the original local basis set representation and the corresponding element positions in the low-dimensional subspace compressed representation;

[0108] Based on the mapping relationship, extract the non-zero matrix elements of the target operator or matrix function from the low-dimensional subspace compressed representation;

[0109] The extracted non-zero matrix elements are written into a sparse matrix data structure according to their row and column indices in the original local basis set representation, thereby reconstructing the sparse matrix form of the target operator or matrix function in the original local basis set representation. After obtaining the low-dimensional subspace compressed representation of the target sparse operator or matrix function, the positions of the matrix elements that need to be recovered in the original local basis set representation of the target sparse operator or matrix function are determined based on the nearest neighbor atom table, spatial truncation conditions, or the non-zero matrix element structure of the target sparse operator. The matrix element positions include row and column indices.

[0110] Based on the nearest neighbor atom table, the local basis function index pairs (i,j) that satisfy the spatial truncation condition are determined. For each matrix element to be recovered... ij The color column corresponding to column index j in the color superposition matrix C is determined based on the color of column index j. If Then matrix elements ij Corresponding compressed representation The element in the i-th row and k-th column is used to extract all non-zero elements of the target sparse operator or matrix function from the compressed representation FC. Then, the extracted non-zero matrix elements are written into a sparse matrix data structure, such as CSR, CSC, COO, or other sparse matrix formats, according to their row and column indices, thereby reconstructing the sparse matrix form of the target sparse operator or matrix function under the original local basis set representation.

[0111] The reconstructed sparse matrix form can be used for charge density updates, total energy calculations, atomic force calculations, orthogonal representation transformations, or other related steps in the self-consistent calculation of Kohn-Sham electronic structure.

[0112] In one specific embodiment, the method is applicable to self-consistent calculations of Kohn-Sham electronic structures using non-orthogonal local basis sets. Under the non-orthogonal local basis set representation, the overlap between basis functions is represented by the overlap matrix S. To perform Löwdin orthogonalization and construct the equivalent Hamiltonian under the orthogonal representation, it is necessary to solve for the inverse square root of the overlap matrix. The equivalent Hamiltonian under orthogonal representations can be expressed as: .in, Local basis set representations typically exhibit spatial locality or finite correlation lengths; therefore, they can be solved using the same steps as described above: coloring grouping, coloring superposition state construction, low-dimensional subspace compression representation, numerical solution, and non-zero matrix element extraction. Specifically, the objective matrix function is taken as:

[0113] ;

[0114] Apply the same procedure as steps S1 and S2 to the overlapping matrix. Constructing the color superposition state matrix ,

[0115] And using the Green's function integral of the matrix function:

[0116] ;

[0117] Krylov subspaces are constructed using block Lanczos iterations, and the solution is obtained by diagonalization in the low-dimensional subspace. Finally, the reconstruction yields... .

[0118] To calculate the equivalent Hamiltonian It can be reused Only need to calculate sequentially as well as You can get The sparse compressed representation can then be obtained by simply extracting the non-zero elements as described above. The sparse representation, compared to the traditional solution... The post-Löwdin orthogonalization process avoids the SpMSpM operation, proving that this method can be applied to the computation of matrix functions. And the process of orthogonalization, including the operators and effective Hamiltonians.

[0119] Under the condition of non-orthogonal basis sets, the original Green's function integral expression of the density matrix is:

[0120] ;

[0121] Substituting this into the Löwdin orthogonalization relation, it can be transformed into the standard form based on the equivalent Hamiltonian H′:

[0122] ;

[0123] Using the same graph coloring strategy as before, we construct a coloring superposition state matrix. By performing subspace compression, we can obtain the compressed projection representation of the density matrix:

[0124] ;

[0125] make Meanwhile, the left side of the integral This can be moved outside the integral, thus the above equation becomes:

[0126] ;

[0127] The simplified expression, The expression and orthogonal basis set Since the structures are the same, the block Lanczos-Krylov subspace dimensionality reduction method can be directly reused to quickly complete the dimensionality reduction under non-orthogonal basis sets. Efficient solution.

[0128] Therefore, this method can be widely applied to the calculation of electronic structure of KS equations under various local basis set representations.

[0129] In one embodiment, the selection of key parameters in this method needs to balance computational accuracy and efficiency, and the specific principles are as follows: cutoff radius The choice directly affects the number of stainers Nc and the calculation accuracy. If the value is too small, important non-zero elements of the density matrix will be lost, leading to a decrease in accuracy. An excessively large value increases the color number Nc, reducing compression efficiency. For insulator systems, The selection of [aspect name] is related to the band gap of the system; the smaller the band gap, the better. The larger the value, the more it depends on the calculated temperature for high-temperature metals. In practical applications, this can be achieved by testing small systems at different temperatures. The optimal value is determined by the energy convergence.

[0130] Furthermore, the subspace dimension *v* controls the approximation accuracy of the Krylov subspace. If *v* is too small, the subspace will fail to fully capture the spectral information of the Hamiltonian; if *v* is too large, it increases computational overhead. In this method, the optimal value of *v* is determined by the energy convergence of the small system at different *v* values.

[0131] The sparse matrix truncation threshold is used to further improve computational efficiency while maintaining accuracy. For matrix elements with very small values ​​in the density matrix, Hamiltonian matrix, and overlap matrix, these can be set to zero to increase sparsity. The truncation threshold is selected based on the premise that it does not significantly affect the total energy and the forces acting on atoms; it is typically set to a threshold value of [value missing]. Order of magnitude. Five computational examples are provided below to verify the actual computational accuracy and efficiency of this method. To verify the feasibility and effectiveness of this method, we chose to implement and test it within the computationally efficient and widely used density functional tight-binding method framework, using atomic orbital basis sets as the local basis sets. Accuracy verification is achieved by comparing the physical properties calculated using the density matrix obtained by this method with those calculated using the diagonalization method, and by comparing the computational efficiency with diagonalization and the SP2 method, verifying the speedup effect of this method relative to diagonalization and previous linear scaling algorithms. The SP2 method is implemented using the NTPoly computational library, maintaining the same sparse matrix truncation threshold as this method. All three methods are computed under the same hardware environment (Intel 16-core CPU, 384GB memory, MPI parallelism).

[0132] Example 1: Validation of the subspace dimension of the water molecule system

[0133] This embodiment aims to verify the convergence behavior of the Krylov subspace dimension v in this method, in order to determine the minimum v value required for practical calculations. This embodiment selects a cubic box system containing 5000 water molecules, with a total number of atomic orbitals N=30000, and uses a Slater-Koster type local orbital basis set. This basis set is non-orthogonal, requiring prior processing of the overlap matrix S for calculation. .

[0134] First, after convergence testing of the cutoff radius, the following parameters were determined: For The calculation of the cutoff radius (Approximately 11.6 Å), color number Nc = 560, subspace dimension v = 22; for the calculation of the density matrix ρ, the cutoff radius... (approximately 9.5 Å), number of colors Nc=310, subspace dimension v=22; sparse matrix truncation threshold set to 2×10⁻ 5 That is, all matrix elements smaller than this value are set to zero in the calculation. The calculation process is as follows: First, the low-dimensional subspace compression representation of the colored superposition state is used to solve the problem using the block Lanczos-Krylov subspace framework. Then construct the equivalent Hamiltonian. Solve again using the same method. Finally, the density matrix is ​​reconstructed. The total energy of the computational system is calculated by gradually increasing the subspace dimension v while keeping other parameters constant.

[0135] The calculation results of this embodiment are shown in Figure 3In (a), the convergence trend of the total energy of the system with increasing subspace dimension is shown. The horizontal axis represents the subspace dimension v, and the vertical axis represents the difference between the total energy calculated using the method of this disclosure and the total energy obtained by diagonalization. The unit is eV. When v=22, the total energy of the system has converged, and further increasing v has a negligible effect on the total energy. At this point, the subspace dimension v=22 is much smaller than the total number of atomic orbitals N=30000, indicating that the computational cost of solving the subspace density matrix is ​​extremely small. This result verifies the core dimensionality reduction idea of ​​this method: the original diagonalization problem of a 30000-dimensional sparse matrix is ​​transformed into the diagonalization problem of a 22-dimensional dense matrix, and the computational cost after dimensionality reduction is compressed from the order of N³ to the order of v³. The main computational overhead is transferred to the SpMM operation in the block Lanczos iteration process, thus achieving a balance between linear scaling and high hardware efficiency.

[0136] Example 2: Verification of the computational accuracy of the water molecule system

[0137] This embodiment aims to verify the computational accuracy of the proposed method in a medium-scale system and compare it with the computational accuracy of the SP2 method. A cubic box system containing 2000 water molecules was selected, with a total of 6000 atoms and a total of N=12000 atomic orbitals. The calculation process and parameters remained the same as in Example 1. The SP2 method was also used to calculate the same system to compare the accuracy of the two methods. The sparse matrix truncation threshold used in the SP2 calculation was also set to 2×10⁻. 5 Based on the complete diagonalization result.

[0138] The results are displayed in Figure 4 middle, Figure 4 In (a) and (b), the horizontal axis represents the atom number, and the vertical axis represents the error of the atom's Milligen charge (RMS_charge), in units of charge (e). Figure 4 In (c) and (d), the horizontal axis represents the atom number and the forces in the three directions (x, y, z), and the vertical axis represents the error (RMS_force) of the atomic forces, in eV / Å. The results show that the total energy error calculated by this method is less than 2 meV, the maximum error of the atomic forces is less than 0.0015 eV / Å, the root mean square error of the force components is 0.00014 eV / Å, and the root mean square error of the Mulliken charge is 3 × 10⁻⁻⁻⁴. 7 e. The accuracy obtained by this method is comparable to that of the SP2 method. The above results show that this method fully meets the accuracy requirements of first-principles calculations while maintaining the computational efficiency of linear scaling.

[0139] Example 3: Efficiency Testing of Water Molecule Systems at Different Scales

[0140] This embodiment aims to verify the computational efficiency and linear scaling properties of the proposed method. Multiple water molecule systems (H₂O) at different scales were tested. n n is expanded from 64 to 30000, and the corresponding number of atomic orbitals is expanded from 384 to 180000. The same parameter settings as in Example 1 are maintained in the calculation, and this method is compared with the second-order trace purification method (SP2) and the conventional diagonalization method (ED).

[0141] The results of this embodiment are shown in Figure 3 middle, Figure 3 (c) illustrates the computational scaling behavior of three methods in self-consistent field iteration of a water molecule system. In the figure, CPU Time refers to the wall time consumed by the CPU during computation. SCFIteration refers to the multiple self-consistent iterations required to update the Hamiltonian in the DFTB self-consistent computation until convergence. Number of Water refers to the number of water molecules contained in the current cubic box. The diagonalization method exhibits typical cubic scaling; although the SP2 method is theoretically linear, its computation time increases significantly for large systems due to the low hardware efficiency of sparse matrix-sparse matrix multiplication (SpMSpM). In contrast, this method not only achieves strict linear scaling but also exhibits the smallest slope in computation time as the system size increases. For a 10,000-molecule water molecule system with 60,000 orbitals, this method requires only 3.3 seconds for a single self-consistent field iteration, achieving an 11-fold speedup compared to the SP2 method's 38 seconds, and this speedup further increases with increasing system size. This efficiency improvement is due to the fact that this method transforms the irregular SpMSpM into regular sparse-dense matrix multiplication (SpMM), which fully leverages the cache prefetching and vectorized instruction performance of modern CPUs.

[0142] At the same time, both the SP2 method and the method disclosed herein require prior calculation. Therefore, for The computational efficiency was also compared; the SP2 method uses the Newton-Schultz method for iterative solution. Its calculation still uses SpMSpM for iteration, therefore the two methods are different. The computational efficiency shows a similar trend to that of the density matrix, which proves that the proposed method can fully realize linear scaling calculation for non-orthogonal systems.

[0143] Example 4: Molecular Dynamics Simulation of Liquid Water System

[0144] This embodiment aims to verify the stability and physical reliability of the proposed method in large-scale molecular dynamics simulations. This embodiment performs NVT ensemble molecular dynamics simulations on a periodic cubic box system containing 10,000 water molecules. The box has a side length of approximately 67 Å, ​​a density of 1.0 g / cm³, a temperature of 298 K, a time step of 0.5 fs, and a total simulation duration of 10 ps. In each molecular dynamics calculation step, the proposed method is used to solve for the electronic structure and calculate the forces acting on atoms, with the parameter settings identical to those in Example 2.

[0145] The results of this embodiment are shown in Figure 5 (a) In the figure, Figure 5 (a) shows the change in the total energy of the system in the last 2 ps after 10 ps of molecular dynamics simulation. The horizontal axis represents the duration of the molecular dynamics simulation (MD simulation time) in ps, and the vertical axis represents the difference between the total energy of the structure and the total energy of the reference structure at each step, averaged over a molecule. The results showed that in the last 2 ps of molecular dynamics simulation, the average total energy change per molecule was less than 1 meV, indicating that the system had reached dynamic equilibrium at this point. Figure 5 (a) The results of the statistical radial distribution function of the molecular dynamics process after system equilibrium are shown and compared with the radial distribution function obtained by experimental testing. The horizontal axis in the figure represents the distance between two atoms in Å, and the vertical axis represents the probability of the distance between two atoms being equal to the value on the horizontal axis. The points of different shapes in the figure represent the OO, OH, and HH radial distribution functions calculated by the method of this disclosure, while the dashed lines of different colors represent the experimentally measured OO, OH, and HH radial distribution functions, respectively. The figure shows that the peak positions and peak shapes of the OO, OH, and HH radial distribution functions calculated by the method of this disclosure are in high agreement with the experimental data, accurately describing the dynamic behavior characteristics of the hydrogen bond network in liquid water. This result shows that the method maintains numerical stability and physical accuracy in molecular dynamics simulations up to 10 ps and can be used for long-scale simulations of complex systems.

[0146] Example 5: Large-scale computation at the million-atom scale

[0147] This embodiment aims to verify the computational capability and stability of the proposed method on an extreme-scale system using a single node. This embodiment tests a system containing one million water molecules (1 million water molecules and 6 million atomic orbitals) on a single computing node (configuration: Intel 96-core CPU, 768GB memory). During the test, it was found that due to memory bandwidth limitations, using all 96 cores for computation intensifies data competition between cores, leading to a decrease in speedup. Therefore, in actual production calculations, a hybrid parallel strategy of 16 cores using MPI and OpenMP was adopted to balance computational efficiency and bandwidth usage. The computational parameters remained consistent with those in Example 1. The calculation results show that a single self-consistent field iteration of the one million water molecule system takes approximately 640 seconds and consumes approximately 700GB of memory.

[0148] The convergence curve calculation results of the system with respect to the number of self-consistent iteration steps are shown in the figure. Figure 5 (b) In the figure, the horizontal axis represents the number of steps in the self-consistent calculation, and the vertical axis represents the difference between the total energy at the self-consistent step and the total energy after convergence. The results show that the self-consistent calculation converges to 10⁻ within 20 steps. 6 Ha precision. This test verifies the numerical stability of the algorithm under extreme memory consumption and long-term iteration, breaking through the bottleneck of single-node computational scale and laying the foundation for subsequent multi-node parallel expansion to the scale of tens of millions of atoms. It should be noted that the above results fully demonstrate the high precision, high efficiency, high stability and high scalability of the Colored Superposition (CSS) method. CSS transforms the problem of solving Kohn-Sham quantum mechanical linear scaling problems from the dilemma of "sparse numerical linear algebra" to "hardware-adapted reduction computation mode", fundamentally alleviating the dilemma of large pre-factor and low hardware efficiency of traditional linear scaling algorithms. In addition, the Colored Superposition method described in this invention is not limited to the self-consistent solution of electronic structure of the KS equation, but can also be extended to computational scenarios of large-scale sparse operators or matrix functions in other scientific computing fields, and is expected to provide a new and efficient solution approach for this type of problem.

[0149] The following specific embodiments are provided in conjunction with the accompanying drawings:

[0150] like Figure 6 As shown, this embodiment of the invention provides a linear scaling electronic structure calculation system based on colored superposition states. The system is applied to the linear scaling electronic structure calculation method based on colored superposition states described in the above embodiment, and is suitable for linear scaling electronic structure calculation based on colored superposition states. The system includes: a colored superposition state construction module 1, a compressed representation generation module 2, a compressed representation solution module 3, and a non-zero element extraction module 4.

[0151] The coloring superposition state construction module 1 is used to color and group the basis function indices under the local basis set representation based on the atomic structure information, local basis function information, and spatial locality or finite correlation length of the target sparse operator or matrix function, using a graph coloring strategy, and constructing a coloring superposition state as a low-dimensional subspace basis based on the grouped basis function indices.

[0152] The compressed representation generation module 2 is connected to the color superposition state construction module 1 and is used to project the target sparse operator or matrix function onto a low-dimensional subspace composed of the color superposition states to form a low-dimensional subspace compressed representation of the target sparse operator or matrix function.

[0153] The compressed representation solving module 3 is connected to the compressed representation generation module 2 and is used to calculate the compression result of the target sparse operator or matrix function acting on the colored superposition matrix under the compressed representation of the low-dimensional subspace.

[0154] The non-zero element extraction module 4 is connected to the compression result calculation module 3. It is used to extract the non-zero matrix elements of the target sparse operator or matrix function in the original local basis set representation from the low-dimensional subspace compressed representation, and reconstruct them into the sparse matrix form of the target sparse operator or matrix function for Kohn-Sham electronic structure self-consistent calculation.

[0155] Since the implementation principle of the linear scaling electronic structure calculation system based on the dye superposition state has been explained in the aforementioned method embodiments, it will not be repeated here.

[0156] The linear scaling electronic structure calculation method based on stained superposition states provided in this invention can be implemented on the terminal side or the server side. For the hardware structure of the electronic terminal, please refer to [link to relevant documentation]. Figure 7 This is a schematic diagram of an optional hardware structure of an electronic terminal 1000 provided in an embodiment of the present invention. The terminal 1000 can be a mobile phone, computer device, tablet device, personal digital processing device, factory back-end processing device, etc. The terminal 1000 includes: at least one processor 1001, a memory 1002, at least one network interface 10010, and a user interface 1009. The various components in the device are coupled together through a bus system 1005. It is understood that the bus system 1005 is used to realize the connection and communication between these components. In addition to a data bus, the bus system 1005 also includes a power bus, a control bus, and a status signal bus. However, for clarity, in... Figure 7 All buses are labeled as bus systems.

[0157] The user interface 1009 may include a monitor, keyboard, mouse, trackball, clicker, button, touchpad, or touch screen.

[0158] It is understood that memory 1002 can be volatile memory or non-volatile memory, or both. Non-volatile memory can be read-only memory (ROM) or programmable read-only memory (PROM), which serves as an external cache. By way of example, but not limitation, many forms of RAM are available, such as static random access memory (SRAM) and synchronous static random access memory (SSRAM). The memories described in the embodiments of this invention are intended to include, but are not limited to, these and any other suitable categories of memory.

[0159] In this embodiment of the invention, the memory 1002 is used to store various types of data to support the operation of the terminal 1000. Examples of this data include: any executable program for operation on the terminal 1000, such as the operating system 10021 and the application program 10022; the operating system 10021 includes various system programs, such as the framework layer, core library layer, driver layer, etc., for implementing various basic services and handling hardware-based tasks. The application program 10022 may include program modules for calculating linear scaling electronic structures based on colored superposition states, solving sparse operators or matrix functions, constructing colored superposition states, calculating compression results, and extracting non-zero matrix elements. A computer program that implements the linear scaling electronic structure calculation method based on colored superposition states provided in this embodiment of the invention may be included in the application program 10022 and executed by the processor 1001.

[0160] The methods disclosed in the above embodiments of the present invention can be applied to or implemented by the processor 1001. The processor 1001 may be an integrated circuit chip with signal processing capabilities. In the implementation process, each step of the above method can be completed by the integrated logic circuit of the hardware in the processor 1001 or by instructions in the form of software. The processor 1001 may be a general-purpose processor, a digital signal processor (DSP), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The processor 1001 can implement or execute the methods, steps, and logic block diagrams disclosed in the embodiments of the present invention. The general-purpose processor 1001 may be a microprocessor or any conventional processor, etc. The steps of the accessory optimization method provided in the embodiments of the present invention can be directly reflected as being executed by a hardware decoding processor, or being executed by a combination of hardware and software modules in the decoding processor. The software module may be located in a storage medium, which is located in a memory. The processor reads the information in the memory and combines it with its hardware to complete the steps of the aforementioned method.

[0161] In an exemplary embodiment, the terminal 1000 may be used to execute the aforementioned method by one or more application-specific integrated circuits (ASICs), DSPs, programmable logic devices (PLDs), or complex programmable logic devices (CPLDs).

[0162] Those skilled in the art will understand that all or part of the steps of the above-described method embodiments can be implemented using computer program-related hardware. The aforementioned computer program can be stored in a computer-readable storage medium. When executed, the program performs the steps of the above-described method embodiments; and the aforementioned storage medium includes various media capable of storing program code, such as ROM, RAM, magnetic disks, or optical disks.

[0163] In the embodiments provided in this application, the computer-readable and writable storage medium may include read-only memory, random access memory, EEPROM, CD-ROM or other optical disc storage devices, disk storage devices or other magnetic storage devices, flash memory, USB flash drive, portable hard drive, or any other medium capable of storing desired program code in the form of instructions or data structures and accessible by a computer. Additionally, any connection may be appropriately referred to as a computer-readable medium. For example, if instructions are transmitted from a website, server, or other remote source using coaxial cable, fiber optic cable, twisted pair, digital subscriber line (DSL), or wireless technologies such as infrared, radio, and microwave, then the coaxial cable, fiber optic cable, twisted pair, DSL, or wireless technologies such as infrared, radio, and microwave are included in the definition of medium. However, it should be understood that computer-readable and writable storage media and data storage media do not include connections, carrier waves, signals, or other transient media, but are intended for non-transient, tangible storage media. The disks and optical discs used in the application include compact optical discs (CDs), laser optical discs, optical discs, digital multifunction optical discs (DVDs), floppy disks, and Blu-ray discs, where disks typically copy data magnetically, while optical discs use lasers to copy data optically.

[0164] In summary, this invention provides a method, system, and terminal for linear scaling electronic structure calculation based on colored superposition states. This method is used to solve for target operators or matrix functions with spatial locality or finite correlation lengths, such as the density matrix, the inverse square root of the overlap matrix, and operators required for orthogonal representation transformations, involved in electronic structure calculations under local basis set representations. Based on the spatial locality of the target operators or matrix functions, this method colors and groups the low-correlation or uncorrelation basis function indices under the local basis set representation, and constructs colored superposition states as low-dimensional subspace bases based on the same-colored basis function indices, thereby forming a low-dimensional subspace compressed representation of the target operators or matrix functions. Under this low-dimensional subspace compressed representation, the compression result of the target operators or matrix functions acting on the colored superposition state matrix is ​​calculated, and the non-zero matrix elements of the target operators or matrix functions under the original local basis set representation are extracted from the low-dimensional subspace compressed representation and reconstructed into a sparse matrix form for use in Kohn-Sham electronic structure self-consistent calculations. Through the above process, this invention can reduce or replace sparse matrix-sparse matrix multiplication (SpMSpM) in traditional linear scaling electronic structure calculations, transforming the relevant calculations into more regular sparse matrix-dense matrix multiplication (SpMM). While maintaining linear scaling complexity, it improves the regularity of the calculation process and its parallel adaptability to modern computing hardware (such as multi-core CPUs / GPUs), thus forming a hardware-efficient linear scaling solution model. Therefore, this invention can improve the solution efficiency of target operators or matrix functions in self-consistent calculations of large-scale electronic structures, and has good industrial application value.

[0165] The above embodiments are merely illustrative of the principles and effects of the present invention and are not intended to limit the invention. Any person skilled in the art can modify or alter the above embodiments without departing from the spirit and scope of the present invention. Therefore, all equivalent modifications or alterations made by those skilled in the art without departing from the spirit and technical concept disclosed in the present invention should still be covered by the claims of the present invention.

Claims

1. A linear scaling electronic structure calculation method based on color superposition states, characterized in that, The method includes: Based on the atomic structure information and local basis function information of the target system, as well as the spatial locality or finite correlation length of the target sparse operator or matrix function, a graph coloring strategy is used to color and group the basis function indices under the local basis set representation, and to construct colored superposition states as a low-dimensional subspace basis based on the grouped basis function indices. The graph coloring strategy includes: constructing an adjacency graph using the local basis function indices as vertices, based on the cutoff radius, spatial locality, finite correlation length, or non-zero matrix element structure of the target sparse operator or matrix function; coloring the adjacency graph such that there are no adjacent edges between local basis functions corresponding to vertices of the same color, or such that the distance between local basis functions corresponding to vertices of the same color is greater than a preset cutoff radius; grouping the local basis functions based on different colors and constructing colored superposition states for each color, forming a colored superposition state matrix of dimension N×Nc as a low-dimensional subspace basis, denoted as […]. Where N is the total number of local basis functions, and Nc is the number of color types; The target sparse operator or matrix function is projected onto a low-dimensional subspace composed of the colored superposition states to form a low-dimensional subspace compressed representation of the target sparse operator or matrix function. The compression result of the target sparse operator or matrix function acting on the colored superposition matrix is ​​calculated under the low-dimensional subspace compression representation; Extract the non-zero matrix elements of the target sparse operator or matrix function under the original local basis set representation from the low-dimensional subspace compressed representation, and reconstruct them into the sparse matrix form of the target sparse operator or matrix function for use in self-consistent calculation of Kohn-Sham electronic structure.

2. The method according to claim 1, characterized in that, The target sparse operator or matrix function includes one or more of the following: density matrix required for self-consistent electronic structure calculation, inverse square root of overlap matrix, and operator-related matrix function required for orthogonal representation transformation.

3. The linear scaling electronic structure calculation method based on dye superposition states according to claim 1, characterized in that, The step of employing a graph coloring strategy to color and group the basis function indices under the local basis set representation, and constructing colored superposition states as low-dimensional subspace bases based on the grouped basis function indices, includes: Grouping by coloring based on basis function indices with low or no correlation; By superimposing the local basis functions corresponding to the basis function indices of the same color, a color superposition state of the corresponding color is constructed. A color superposition state matrix is ​​formed by multiple color superposition states, and the color superposition state matrix is ​​used as a basis for a low-dimensional subspace.

4. The linear scaling electronic structure calculation method based on dye superposition states according to claim 1, characterized in that, Methods for constructing adjacency graphs include: Read the atomic structure information and local basis function information of the target system; For all target sparse operators or matrix functions, the nearest neighbor atoms within the truncation radius are identified, and a nearest neighbor atom table is created for each atom to record its nearest neighbor information; An adjacency graph is constructed using local basis function indices as vertices, and adjacency edges are established between corresponding local basis function indices based on nearest neighbor relationships.

5. The linear scaling electronic structure calculation method based on dye superposition states according to claim 1, characterized in that, The low-dimensional subspace compression representation includes the compression result obtained by applying a target sparse operator or matrix function to a colored superposition matrix. If the target operator or matrix function is denoted as F, then the low-dimensional subspace compression representation includes the compression result obtained by applying F to a colored superposition matrix. The first dense matrix represented by and by The second dense matrix is ​​represented by the first dense matrix, and F is formed by multiplying the first dense matrix and the second dense matrix. The compressed decomposition representation in this form.

6. The linear scaling electronic structure calculation method based on dye superposition states according to claim 1, characterized in that, The method for calculating the compression result of the target sparse operator or matrix function acting on the colored superposition matrix employs a block Lanczos-Krylov subspace approach. This approach includes: using the colored superposition matrix as multiple right-hand vectors, mapping the calculation process of the target sparse operator or matrix function acting on the colored superposition matrix to a Krylov subspace through sparse matrix-dense matrix multiplication and block Lanczos iteration between the matrix to be processed and the multiple right-hand vectors; calculating the target matrix function within the Krylov subspace; and mapping the calculation result to the compression result of the target sparse operator or matrix function acting on the colored superposition matrix.

7. The method according to claim 1, characterized in that, The step of extracting the non-zero matrix elements of the target sparse operator or matrix function in the original local basis set representation from the compressed representation of the low-dimensional subspace, and reconstructing it into the sparse matrix form of the target sparse operator or matrix function, includes: Based on the nearest neighbor atom table, determine the matrix element positions that satisfy the spatial truncation condition under the original local basis set representation of the target operator or matrix function, wherein the matrix element positions include the corresponding row index and column index; Based on the color superposition state, establish the mapping relationship between the matrix element positions under the original local basis set representation and the corresponding element positions in the low-dimensional subspace compressed representation; Based on the mapping relationship, extract the non-zero matrix elements of the target operator or matrix function from the low-dimensional subspace compressed representation; The extracted non-zero matrix elements are written into the sparse matrix data structure according to their row and column indices under the original local basis set representation, thereby reconstructing the sparse matrix form of the target operator or matrix function under the original local basis set representation.

8. A linear scaling electronic structure calculation system based on dye superposition states, characterized in that, The system includes: A coloring superposition state construction module is used to color and group basis function indices under the local basis set representation based on the atomic structure information, local basis function information, and spatial locality or finite correlation length of the target sparse operator or matrix function, using a graph coloring strategy. Based on the grouped basis function indices, a coloring superposition state is constructed as a low-dimensional subspace basis. The graph coloring strategy includes: constructing an adjacency graph using local basis function indices as vertices, based on the cutoff radius, spatial locality, finite correlation length, or non-zero matrix element structure of the target sparse operator or matrix function; coloring the adjacency graph such that there are no adjacent edges between local basis functions corresponding to vertices of the same color, or that the distance between local basis functions corresponding to vertices of the same color is greater than a preset cutoff radius; grouping the local basis functions based on different colors and constructing a coloring superposition state for each color, forming a coloring superposition state matrix of dimension N×Nc as a low-dimensional subspace basis, denoted as […]. Where N is the total number of local basis functions, and Nc is the number of color types; The compressed representation generation module is used to project the target sparse operator or matrix function onto a low-dimensional subspace composed of the colored superposition states, thereby forming a low-dimensional subspace compressed representation of the target sparse operator or matrix function. The compression representation solution module is used to calculate the compression result of the target sparse operator or matrix function acting on the colored superposition matrix under the compression representation of the low-dimensional subspace; The non-zero element extraction module is used to extract the non-zero matrix elements of the target sparse operator or matrix function under the original local basis set representation from the low-dimensional subspace compressed representation, and reconstruct them into the sparse matrix form of the target sparse operator or matrix function for use in Kohn-Sham electronic structure self-consistent calculation.

9. An electronic terminal, characterized in that, include: One or more memories and one or more processors; The one or more memories are used to store computer programs; The one or more processors are connected to the memory and are used to run the computer program to perform the method for solving linear scaling matrix functions based on color superposition subspace projection as described in any one of claims 1 to 7.

Images