A method for accelerating phonon spectrum calculation

Abstract

The present application relates to lattice dynamics, in particular to a kind of phonon spectrum acceleration calculation method, constructs q Point grid, for each q Wave vector corresponding to the point in the grid, the Dyson equation is solved by self-consistent iteration, in each step of iteration calculation of Dyson equation, the matrix of transpose Khatri-Rao product in Sternheimer equation related to the current q Point is decomposed by ISDF, and the matrix dimension is compressed from high dimension to low rank representation;The number of Sternheimer equation to be solved is reduced by introducing ACP operator and Chebyshev node interpolation;The solution result of Sternheimer equation is used in the self-consistent calculation of Dyson equation, and the condition number of iteration matrix is improved by introducing ACP operator, and the convergence speed of self-consistent calculation is accelerated;The technical scheme provided by the present application can effectively overcome the defect that phonon spectrum cannot be calculated quickly and accurately.

Description

Technical Field

[0001] This invention relates to lattice dynamics, and more specifically to a method for accelerating the calculation of phonon spectra. Background Technology

[0002] The lattice dynamics of solids largely determine their various physical properties, including infrared spectroscopy, inelastic neutron scattering, specific heat, thermal conductivity, and effects related to electro-acoustic coupling. Therefore, phonon properties are a core component of first-principles studies, with the phonon spectrum being particularly important as a key quantity for characterizing and understanding the thermal and electrical properties, structural stability, and superconductivity of materials. Currently, within the plane-wave basis set framework, the self-consistent linear response method based on density functional perturbation theory (DFPT), i.e., the self-consistent iterative solution of the Dyson equation, is widely considered one of the most accurate approaches to calculating phonon spectra.

[0003] However, existing technologies have the following significant drawbacks: 1) Extremely high computational complexity: The traditional DFPT method requires a large number of iterative calculations when solving the Dyson equation in a self-consistent iterative manner. Its computational cost increases to the fourth power with the number of atoms in the system, making it extremely difficult or even impossible to calculate the phonon spectrum of large or complex systems. 2) Poor convergence: In the iterative solution process, especially for systems with complex electronic structures, the traditional DFPT method may converge slowly or be unstable, requiring a lot of computational resources and being inefficient. 3) High resource consumption: The high computational complexity and poor convergence mentioned above directly lead to a huge demand for computing hardware (especially memory and CPU), which limits the widespread application of the DFPT method in industrial-grade material design.

[0004] The root of the aforementioned technical problems lies in the fact that the solution of the Dyson equation and the Sternheimer equation in the DFPT method involves the direct processing and repeated calculation of the high-dimensional electron density response function, which constitutes the main bottleneck of the entire calculation process. Summary of the Invention

[0005] (a) Technical problems to be solved

[0006] In view of the above-mentioned shortcomings of the existing technology, the present invention provides a method for accelerating phonon spectrum calculation, which can effectively overcome the defects of the existing technology in that it cannot calculate phonon spectrum quickly and accurately.

[0007] (II) Technical Solution

[0008] To achieve the above objectives, the present invention provides the following technical solution: A method for accelerating phonon spectrum calculations includes the following steps: S1. Input the optimized crystal structure and perform standard density functional theory (DFT) self-consistent calculations to obtain the ground state electronic wave function and energy. S2. Construct a q-point grid. For the wave vector corresponding to each q-point in the grid, perform a self-consistent iterative solution to the Dyson equation. In each iterative calculation of the Dyson equation, perform ISDF decomposition on the matrix formed by the transpose Khatri-Rao product of the Sternheimer equation related to the current q-point, compressing the matrix dimension from high dimension to low rank representation. S21. Use the weighted K-Means algorithm to calculate the expansion coefficients of the matrix, i.e. the interpolation points of the ISDF; S22. Calculate the auxiliary basis function, i.e. the interpolation vector of the ISDF, and finally decompose the matrix into a linear combination of interpolation points and interpolation vectors; S3. Reduce the number of Sternheimer equations to be solved by interpolating Chebyshev nodes, thereby speeding up the solution of Sternheimer equations; S4. Apply the solution results of the Sternheimer equation to the self-consistent calculation of the Dyson equation. At the same time, introduce the ACP operator to improve the condition number of the iteration matrix and check the convergence condition. If it does not converge, return to S2. If it converges, output the response of the electron density to the atomic displacement perturbation. S5. Calculate the force constant matrix of the grid at point q; S6. Calculate the phonon frequency at any q point by interpolation based on the force constant matrix of the q-point grid, and plot the phonon spectrum on the high symmetry path of the crystal.

[0009] Preferably, in S2, in each iterative calculation of the Dyson equation, the matrix formed by the transpose Khatri-Rao product of the Sternheimer equation related to the current q point is decomposed using ISDF, including: The Sternheimer equation is expressed as follows: ; Where Q is a projection operator over the unoccupied state space. I is the unit operator. , Let N be the wavefunction of the i-th electronic state and its complex conjugate. e Let i be the total number of electrons, i = 1, 2, ..., N e , Indicates the occupancy space projection. Let H be the energy of the i-th electronic state, and H be the Hamiltonian. Let g be the intermediate variable for the electron density response of the i-th electronic state. j This is the source of the perturbation caused by the j-th atomic shift; In the above formula, In the form of the transposed Khatri-Rao product, The Sternheimer equations contain redundant information; ISDF decomposition is used to speed up the computation. make For element M ij Perform ISDF decomposition on the constructed matrix M: ; in, g is the perturbation source induced by the i-th electronic state at the j-th atomic shift. j Below, the matrix element at position r in real space, Let be the interpolation vector at the μ-th ISDF interpolation point at real space location r. g is the perturbation source induced by the i-th electronic state at the j-th atomic shift. j Below, the matrix element at the μ-th ISDF interpolation point at real space location r, The total number of ISDF interpolation points, μ = 1, 2, ..., N μ , The wave function of the i-th electronic state The value at the μ-th ISDF interpolation point at the real space location r. The perturbation source g induced by the displacement of the j-th atom j The value at the μth ISDF interpolation point at the real space location r; After performing ISDF decomposition on matrix M, the Sternheimer equation becomes: ; in, Let be the intermediate variable of the electron density response at the i-th electronic state and the μ-th ISDF interpolation point. Essentially, it is the projection of the electron density response to atomic displacement perturbations in a low-dimensional space, and it is the variable to be solved in the Sternheimer equation. Let be the interpolation vector for the μ-th ISDF interpolation point.

[0010] Preferably, in S21, the weighted K-Means algorithm is used to calculate the expansion coefficients of the matrix, i.e., the interpolation points of the ISDF, including: Since the transpose Khatri-Rao product in the Sternheimer equation is very large, using the ordinary K-Means algorithm is inefficient. However, K-Means clustering can be performed based on real-space grid points, and the point sampling process becomes solving an optimization problem: ; in, Represents the k-th real space lattice point Belongs to the μ-th cluster , For the k-th real space lattice point The weight values, which incorporate Kohn-Sham orbital information, For the μ-th cluster The centroid, i.e. the ISDF interpolation point, is obtained by the weighted average of all real space lattice points within the cluster: , Represents the k-th real space lattice point With the μth cluster center of mass The square of the Euclidean distance between them.

[0011] Preferably, the weighting functions include four types: product squared weighting function (PSM), product weighting function (PM), sum squared weighting function (SSM), and sum weighting function (SM). 1) Product Squared Weighting Function (PSM): ; 2) Product weight function PM: ; 3) Sum-of-Square Weights Function (SSM): ; 4) and weight function SM: ; Among them, u m Let f be the electron density response corresponding to the m-th perturbation source, where N is the total number of perturbation sources, m = 1, 2, ..., N. hxc The commutative correlation functional tensor is a fourth-order tensor.

[0012] Preferably, in S3, the number of Sternheimer equations to be solved is reduced by interpolation through Chebyshev nodes, including: By using Chebyshev node interpolation to occupy the bandgap, the number of Sternheimer equations to be solved is reduced, and the Sternheimer equations become: ; in, Let c be the energy of the Chebyshev interpolation node. For the c-th Chebyshev interpolation node and the μ-th ISDF interpolation point, the intermediate variable representing the electron density response to atomic displacement perturbations is essentially the projection of the electron density response to atomic displacement perturbations onto the energy-space dual low-dimensional space. N represents the variable to be solved in the Sternheimer equation. cLet c be the total number of Chebyshev interpolation nodes, where c = 1, 2, ..., N. c .

[0013] Preferably, the ACP operator is introduced in S4 to improve the condition number of the iteration matrix, including: make ; in, Let g be the irreducible polarizability operator and g be the perturbation source. The response is the irreducible electron density response, and hc represents Hermitian conjugation; The irreducible polarizability operator Acting on Above: ; in, Let r be the electron density response at real space position r caused by the m-th perturbation source induced by the atomic displacement at real space position r'. For the independent particle approximation, the electron density response correlation function generated at real space positions r' and r''' by the perturbation source caused by the atomic displacements at real space positions r' and r'''. For the perturbation source caused by the atomic displacement at real space position r'', the exchange correlation functional tensor of the electron at real space position r''' is... The electron density response at the real space position r'' is generated by the perturbation source caused by the atomic displacement at the real space position r''. The self-consistent calculation of the above equation only requires the diagonal elements of the electron density response u corresponding to the m-th perturbation source, therefore: ; in, The initial electron density response generated at the real space position r without exchange correlation is the m-th perturbation source caused by the atomic displacement at the real space position r'. make , Then we have: ; For any given m-th perturbation source, the electron density response u m Both can construct an ACP operator. This makes the Dyson equation become: ; Among them, the ACP operator It is a low-rank operator, relative to the solution u of the Dyson equation. m Adaptive change.

[0014] Preferably, in step S6, the phonon frequency at any q-point is calculated by interpolation based on the force constant matrix of the q-point grid, and the phonon spectrum is plotted on the high-symmetry path of the crystal, including: S61. Select a highly symmetric path in the Brillouin zone and take a series of q points along the path; S62. Perform Fourier interpolation on the force constant matrix of the solved q-point mesh: obtain the force constant matrix in real space through inverse Fourier transform, and then perform Fourier transform on the force constant matrix in real space to obtain the dynamic matrix of the wave vector position corresponding to any q-point in reciprocal space. S63. Diagonalize the dynamic matrix to obtain the phonon frequency at each q point on the high-symmetric path. Connect the interpolated phonon frequencies at different q points to obtain the phonon spectrum, which is represented as a phonon dispersion curve.

[0015] (III) Beneficial Effects

[0016] Compared with the prior art, the phonon spectrum acceleration calculation method provided by the present invention has the following beneficial effects: 1) Significantly improved computational speed: By shifting the most time-consuming calculation of the response of electron density to atomic displacement perturbation from high-dimensional space to low-dimensional space, the amount of computation is greatly reduced, effectively improving the computational speed; 2) Controllable accuracy: The computational accuracy can be controlled by adjusting the parameters of the weighted K-Means algorithm. Users can flexibly balance computational efficiency and accuracy according to actual needs. In most cases, the accuracy required for scientific research and engineering can be achieved by using far less data than the traditional DFPT method. 3) Enhanced convergence and stability: The low-rank approximation provided by the ISDF decomposition based on the weighted K-Means algorithm plays a role in numerical "smoothing", which can achieve the convergence of the Dyson equation in fewer iterations and the convergence process is more stable. 4) Reduced memory usage: Instead of storing and manipulating the high-dimensional transposed Khatri-Rao product matrix in memory, the interpolation points and interpolation vectors of the low-dimensional ISDF are stored, enabling the present invention to calculate more complex material systems with a larger number of atoms under the same hardware conditions. Attached Figure Description

[0017] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are merely some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without any creative effort.

[0018] Figure 1This is a schematic diagram of the process of the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0020] The core of this invention lies in the creative introduction of the Adaptive Compressible Polarimetry (ACP) operator into the calculation of the phonon spectrum, and the reduction of the number of Sternheimer equations to be solved at each q-point to O(N) through Chebyshev interpolation and interpolation separable density fitting (ISDF) decomposition. e This significantly accelerates the self-consistent convergence speed of the Dyson equation and increases the computational speed by orders of magnitude.

[0021] The following describes the specific process of the phonon spectrum acceleration calculation method provided by this invention with specific examples (e.g.) Figure 1 (as shown) and technical effects.

[0022] S1. Input the optimized crystal structure and perform standard density functional theory (DFT) self-consistent calculations to obtain the ground state electronic wave function and energy.

[0023] S2. Construct a grid of q points (q1×q2×q3, where q1, q2, and q3 are all positive integers). For the wave vector corresponding to each q point in the grid, perform a self-consistent iterative solution to the Dyson equation. In each iterative calculation of the Dyson equation, perform ISDF decomposition on the matrix formed by the transpose Khatri-Rao product of the Sternheimer equation related to the current q point, compressing the matrix dimension from high dimension to low rank representation. S21. Use the weighted K-Means algorithm to calculate the expansion coefficients of the matrix, i.e. the interpolation points of the ISDF; S22. Calculate the auxiliary basis function, i.e. the interpolation vector of the ISDF, and finally decompose the matrix into a linear combination of interpolation points and interpolation vectors.

[0024] Specifically, in S2, during each iterative calculation of the Dyson equation, the matrix formed by the transpose of the Khatri-Rao product in the Sternheimer equation related to the current q point is decomposed using ISDF, including: The Sternheimer equation is expressed as follows: ; Where Q is a projection operator over the unoccupied state space. I is the unit operator. , Let N be the wavefunction of the i-th electronic state and its complex conjugate. e Let i be the total number of electrons, i = 1, 2, ..., N e , Indicates the occupancy space projection. Let H be the energy of the i-th electronic state, and H be the Hamiltonian. Let g be the intermediate variable for the electron density response of the i-th electronic state. j This is the source of the perturbation caused by the j-th atomic shift; In the above formula, In the form of the transposed Khatri-Rao product, The Sternheimer equations contain redundant information; ISDF decomposition is used to speed up the computation. make For element M ij Perform ISDF decomposition on the constructed matrix M: ; in, g is the perturbation source induced by the i-th electronic state at the j-th atomic shift. j Below, the matrix element at position r in real space, Let be the interpolation vector at the μ-th ISDF interpolation point at real space location r. g is the perturbation source induced by the i-th electronic state at the j-th atomic shift. j Below, the matrix element at the μ-th ISDF interpolation point at real space location r, The total number of ISDF interpolation points, μ = 1, 2, ..., N μ , The wave function of the i-th electronic state The value at the μ-th ISDF interpolation point at the real space location r. The perturbation source g induced by the displacement of the j-th atom j The value at the μth ISDF interpolation point at the real space location r; After performing ISDF decomposition on matrix M, the Sternheimer equation becomes: ; in, Let be the intermediate variable of the electron density response at the i-th electronic state and the μ-th ISDF interpolation point. Essentially, it is the projection of the electron density response to atomic displacement perturbations in a low-dimensional space, and it is the variable to be solved in the Sternheimer equation. Let be the interpolation vector for the μ-th ISDF interpolation point.

[0025] The traditional ACP algorithm uses the QRCP method to calculate the interpolation points of the ISDF, with a calculation scale of... This uses a sampling scheme that is much faster than the QRCP method, and the calculation scale is... The weighted K-Means algorithm accelerates the solution of the Sternheimer equations at different q points, ultimately speeding up the calculation of the phonon spectrum.

[0026] Specifically, in S21, the weighted K-Means algorithm is used to calculate the expansion coefficients of the matrix, i.e., the interpolation points of the ISDF, including: Since the transpose Khatri-Rao product in the Sternheimer equation is very large, using the ordinary K-Means algorithm is inefficient. However, K-Means clustering can be performed based on real-space grid points, and the point sampling process becomes solving an optimization problem: ; in, Represents the k-th real space lattice point Belongs to the μ-th cluster , For the k-th real space lattice point The weight values, which incorporate Kohn-Sham orbital information, For the μ-th cluster The centroid, i.e. the ISDF interpolation point, is obtained by the weighted average of all real space lattice points within the cluster: , Represents the k-th real space lattice point With the μth cluster center of mass The square of the Euclidean distance between them.

[0027] The above technical solutions rely on empirical methods for selecting weighting functions, but choosing a suitable weighting function is crucial. In this application's technical solution, the weighting functions include four types: Product Square Weighting Function (PSM), Product Weighting Function (PM), Sum Square Weighting Function (SSM), and Sum Weighting Function (SM). 1) Product Squared Weighting Function (PSM): ; 2) Product weight function PM: ; 3) Sum-of-Square Weights Function (SSM): ; 4) and weight function SM: ; Among them, um Let f be the electron density response corresponding to the m-th perturbation source, where N is the total number of perturbation sources, m = 1, 2, ..., N. hxc The commutative correlation functional tensor is a fourth-order tensor.

[0028] S3. Reduce the number of Sternheimer equations to be solved by interpolating through Chebyshev nodes, thereby speeding up the solution of Sternheimer equations.

[0029] Specifically, S3 reduces the number of Sternheimer equations to be solved through Chebyshev node interpolation, including: By using Chebyshev node interpolation to occupy the bandgap, the number of Sternheimer equations to be solved is reduced, and the Sternheimer equations become: ; in, Let c be the energy of the Chebyshev interpolation node. For the c-th Chebyshev interpolation node and the μ-th ISDF interpolation point, the intermediate variable representing the electron density response to atomic displacement perturbations is essentially the projection of the electron density response to atomic displacement perturbations onto the energy-space dual low-dimensional space. N represents the variable to be solved in the Sternheimer equation. c Let c be the total number of Chebyshev interpolation nodes, where c = 1, 2, ..., N. c .

[0030] S4. Apply the solution results of the Sternheimer equation to the self-consistent calculation of the Dyson equation. At the same time, introduce the ACP operator to improve the condition number of the iteration matrix and check the convergence condition. If it does not converge, return to S2. If it converges, output the response of the electron density to the atomic displacement perturbation.

[0031] Specifically, S4 introduces the ACP operator to improve the condition number of the iteration matrix, including: make ; in, Let g be the irreducible polarizability operator and g be the perturbation source. The response is the irreducible electron density response, and hc represents Hermitian conjugation; The irreducible polarizability operator Acting on Above: ; in, Let r be the electron density response at real space position r caused by the m-th perturbation source induced by the atomic displacement at real space position r'. For the independent particle approximation, the electron density response correlation function generated at real space positions r' and r''' by the perturbation source caused by the atomic displacements at real space positions r' and r'''. For the perturbation source caused by the atomic displacement at real space position r'', the exchange correlation functional tensor of the electron at real space position r''' is... The electron density response at the real space position r'' is generated by the perturbation source caused by the atomic displacement at the real space position r''. The self-consistent calculation of the above equation only requires the diagonal elements of the electron density response u corresponding to the m-th perturbation source, therefore: ; in, The initial electron density response generated at the real space position r without exchange correlation is the m-th perturbation source caused by the atomic displacement at the real space position r'. make , Then we have: ; For any given m-th perturbation source, the electron density response u m Both can construct an ACP operator. This makes the Dyson equation become: ; Among them, the ACP operator It is a low-rank operator, relative to the solution u of the Dyson equation. m Adaptive change.

[0032] S5. Calculate the force constant matrix of the q-point mesh.

[0033] S6. Calculate the phonon frequency (usually in cm) at any q-point by interpolation based on the force constant matrix of the q-point grid. -1 (or THz), and plot the phonon spectrum along the high symmetry path of the crystal, including: S61. Select a highly symmetric path in the Brillouin zone and take a series of q points along the path; S62. Perform Fourier interpolation on the force constant matrix of the solved q-point mesh: obtain the force constant matrix in real space through inverse Fourier transform, and then perform Fourier transform on the force constant matrix in real space to obtain the dynamic matrix of the wave vector position corresponding to any q-point in reciprocal space. S63. Diagonalize the dynamic matrix to obtain the phonon frequency at each q point on the high-symmetric path, and interpolate the phonon frequencies at different q points to obtain the phonon spectrum (plot the phonon frequency as a function of the q point position), which is represented as a phonon dispersion curve.

[0034] The following example, using the calculation of the phonon spectrum of a typical semiconductor material (such as single-crystal silicon), illustrates in detail the technical effects achieved by this invention.

[0035] Example of performance comparison: On the same computing platform, the traditional DFPT method takes a long time to calculate the complete phonon spectrum of silicon, while the technical solution provided by this invention takes less time, and the obtained phonon frequencies are in good agreement with the output results of mainstream computing software (such as Quantum ESPRESSO).

[0036] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions will not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for accelerating phonon spectrum calculation, characterized in that: Includes the following steps: S1. Input the optimized crystal structure and perform standard density functional theory (DFT) self-consistent calculations to obtain the ground state electronic wave function and energy. S2. Construct a q-point grid. For the wave vector corresponding to each q-point in the grid, perform a self-consistent iterative solution to the Dyson equation. In each iterative calculation of the Dyson equation, perform ISDF decomposition on the matrix formed by the transpose Khatri-Rao product of the Sternheimer equation related to the current q-point, compressing the matrix dimension from high dimension to low rank representation. S21. Use the weighted K-Means algorithm to calculate the expansion coefficients of the matrix, i.e. the interpolation points of the ISDF; S22. Calculate the auxiliary basis function, i.e. the interpolation vector of the ISDF, and finally decompose the matrix into a linear combination of interpolation points and interpolation vectors; S3. Reduce the number of Sternheimer equations to be solved by interpolating Chebyshev nodes, thereby speeding up the solution of Sternheimer equations; S4. Apply the solution results of the Sternheimer equation to the self-consistent calculation of the Dyson equation. At the same time, introduce the ACP operator to improve the condition number of the iteration matrix and check the convergence condition. If it does not converge, return to S2. If it converges, output the response of the electron density to the atomic displacement perturbation. S5. Calculate the force constant matrix of the grid at point q; S6. Calculate the phonon frequency at any q-point by interpolation based on the force constant matrix of the q-point grid, and plot the phonon spectrum on the high-symmetry path of the crystal, including: S61. Select a highly symmetric path in the Brillouin zone and take a series of q points along the path; S62. Perform Fourier interpolation on the force constant matrix of the solved q-point mesh: obtain the force constant matrix in real space through inverse Fourier transform, and then perform Fourier transform on the force constant matrix in real space to obtain the dynamic matrix of the wave vector position corresponding to any q-point in reciprocal space. S63. Diagonalize the dynamic matrix to obtain the phonon frequency at each q point on the high-symmetry path, and interpolate the phonon frequencies at different q points to obtain the phonon dispersion curve showing the change of phonon frequency with the position of the q-wave vector on the high-symmetry path.

2. The phonon spectrum acceleration calculation method according to claim 1, characterized in that: In S2, during each iterative calculation of the Dyson equation, the matrix formed by the transpose of the Khatri-Rao product in the Sternheimer equation related to the current q point is decomposed using ISDF, including: The Sternheimer equation is expressed as follows: ； Where Q is a projection operator over the unoccupied state space. I is the unit operator. , Let N be the wavefunction of the i-th electronic state and its complex conjugate. e Let i be the total number of electrons, i = 1, 2, ..., N e , Indicates the occupancy space projection. Let H be the energy of the i-th electronic state, and H be the Hamiltonian. Let g be the intermediate variable for the electron density response of the i-th electronic state. j This is the source of the perturbation caused by the j-th atomic shift; In the above formula, In the form of the transposed Khatri-Rao product, The Sternheimer equations contain redundant information; ISDF decomposition is used to speed up the computation. make For element M ij Perform ISDF decomposition on the constructed matrix M: ； in, g is the perturbation source induced by the i-th electronic state at the j-th atomic shift. j Below, the matrix element at position r in real space, Let be the interpolation vector at the μ-th ISDF interpolation point at real space location r. g is the perturbation source induced by the i-th electronic state at the j-th atomic shift. j Below, the matrix element at the μ-th ISDF interpolation point at real space location r, The total number of ISDF interpolation points, μ = 1, 2, ..., N μ , The wave function of the i-th electronic state The value at the μ-th ISDF interpolation point at the real space location r. The perturbation source g induced by the displacement of the j-th atom j The value at the μth ISDF interpolation point at the real space location r; After performing ISDF decomposition on matrix M, the Sternheimer equation becomes: ； in, Let be the intermediate variable of the electron density response at the i-th electronic state and the μ-th ISDF interpolation point. Essentially, it is the projection of the electron density response to atomic displacement perturbations in a low-dimensional space, and it is the variable to be solved in the Sternheimer equation. Let be the interpolation vector for the μ-th ISDF interpolation point.

3. The phonon spectrum acceleration calculation method according to claim 2, characterized in that: In S21, the weighted K-Means algorithm is used to calculate the expansion coefficients of the matrix, i.e., the interpolation points of the ISDF, including: Since the transpose Khatri-Rao product in the Sternheimer equation is very large, using the ordinary K-Means algorithm is inefficient. However, K-Means clustering can be performed based on real-space grid points, and the point sampling process becomes solving an optimization problem: ； in, Represents the k-th real space lattice point Belongs to the μ-th cluster , For the k-th real space lattice point The weight values, which incorporate Kohn-Sham orbital information, For the μ-th cluster The centroid, i.e. the ISDF interpolation point, is obtained by the weighted average of all real space lattice points within the cluster: , Represents the k-th real space lattice point With the μ-th cluster center of mass The square of the Euclidean distance between them.

4. The phonon spectrum acceleration calculation method according to claim 3, characterized in that: Weighting functions include four types: product squared weighting function (PSM), product weighting function (PM), sum squared weighting function (SSM), and sum weighting function (SM). 1) Product Squared Weighting Function (PSM): ； 2) Product weight function PM: ； 3) Sum-of-Square Weights Function (SSM): ； 4) and weight function SM: ； Among them, u m Let f be the electron density response corresponding to the m-th perturbation source, where N is the total number of perturbation sources, m = 1, 2, ..., N. hxc The commutative correlation functional tensor is a fourth-order tensor.

5. The phonon spectrum acceleration calculation method according to claim 3, characterized in that: In S3, the number of Sternheimer equations to be solved is reduced through Chebyshev node interpolation, including: By using Chebyshev node interpolation to occupy the bandgap, the number of Sternheimer equations to be solved is reduced, and the Sternheimer equations become: ； in, Let c be the energy of the Chebyshev interpolation node. For the c-th Chebyshev interpolation node and the μ-th ISDF interpolation point, the intermediate variable representing the electron density response to atomic displacement perturbations is essentially the projection of the electron density response to atomic displacement perturbations onto the "energy-space" dual low-dimensional space. N represents the variable to be solved in the Sternheimer equation. c Let c be the total number of Chebyshev interpolation nodes, where c = 1, 2, ..., N. c .

6. The phonon spectrum acceleration calculation method according to claim 5, characterized in that: S4 introduces the ACP operator to improve the condition number of the iteration matrix, including: make ; in, Let g be the irreducible polarizability operator and g be the perturbation source. The response is the irreducible electron density response, and hc represents Hermitian conjugation; The irreducible polarizability operator Acting on Above: ； in, Let r be the electron density response at real space position r caused by the m-th perturbation source induced by the atomic displacement at real space position r'. For the independent particle approximation, the electron density response correlation function generated at real space positions r' and r''' by the perturbation source caused by the atomic displacements at real space positions r' and r'''. For the perturbation source caused by the atomic displacement at real space position r'', the exchange correlation functional tensor of the electron at real space position r''' is... The electron density response at the real space position r'' is generated by the perturbation source caused by the atomic displacement at the real space position r''. The self-consistent calculation of the above equation only requires the diagonal elements of the electron density response u corresponding to the m-th perturbation source, therefore: ； in, The initial electron density response generated at the real space position r without exchange correlation is the m-th perturbation source caused by the atomic displacement at the real space position r'. make , Then we have: ； For any given m-th perturbation source, the electron density response u m Both can construct an ACP operator. This makes the Dyson equation become: ； Among them, the ACP operator It is a low-rank operator, relative to the solution u of the Dyson equation. m Adaptive change.

Images