Material performance prediction and high-throughput screening method based on deep learning force field model

Abstract

This invention relates to the fields of computational materials science and materials informatics, and discloses a method for predicting and screening material properties based on a deep learning force field model. The method includes: inputting material structure data into a first-principles calculation module to generate an initial dataset; removing outlier samples through multiple physical plausibility checks; inputting the screened, consistent dataset into a force field model based on an isovariant graph neural network for training; updating model parameters by constructing a joint loss function using an automatic differentiation mechanism and physical conservation constraints; inputting candidate materials into the trained model to obtain energy and forces; outputting mechanically stable configurations after structural relaxation; calculating the multi-dimensional physical properties of these configurations; and weighting and ranking them using a multi-index comprehensive scoring function. This invention solves the problems of easy force divergence in traditional force fields and the high cost of first-principles calculations, achieving low-cost, high-throughput screening of unknown materials while ensuring physical consistency and accuracy.

Description

Technical Field

[0001] This invention relates to the fields of computational materials science and materials informatics, specifically to a method for predicting and screening material properties based on a deep learning force field model. Background Technology

[0002] With the rapid development of high-end manufacturing fields such as new microelectronics and new energy, the demand for new materials research and development for specific application scenarios is constantly increasing. High-throughput material screening technology based on computer simulation can pre-evaluate the performance of massive virtual structures in digital space, thereby effectively shortening the material research and development cycle. Currently, high-throughput computing mainly relies on two types of methods: first-principles calculations and traditional empirical force fields to obtain the physical quantities of materials and guide structural evolution.

[0003] However, existing technologies have inherent limitations when facing large-scale space exploration tasks involving a vast array of unknown elements. First-principles calculations, which require solving the Schrödinger equation for multi-electron systems, offer high microscopic physical accuracy but consume enormous computational resources. The computational cost often increases exponentially with the system size, making it difficult to support rapid evaluation of massive amounts of candidate material data in practical engineering. To reduce computational costs, researchers often use traditional empirical force fields for large-scale screening. However, traditional force fields heavily rely on pre-defined fixed mathematical functions and empirical parameters fitted to specific systems. When faced with novel alloys containing unknown chemical bonding characteristics or complex configurations with extreme structural distortions, these empirical parameters often fail, making it difficult for the force field model to maintain fundamental physical consistency and easily leading to atomic force state divergence during calculations. The contradiction between excessively high computational costs and weak force field generalization capabilities results in existing high-throughput screening methods not only having insufficient performance prediction accuracy when dealing with multi-source heterogeneous material structures, but also frequently failing to effectively converge the structural relaxation processes of the underlying atoms, making it difficult to output truly reliable mechanically stable configurations. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a material performance prediction and high-throughput screening method based on a deep learning force field model. The technical problem it solves is that existing high-throughput material screening methods are unable to maintain physical consistency and are prone to force state divergence when faced with large-scale unknown element combinations and extreme structural distortions. Furthermore, the computational cost of first-principles calculations is too high to support the evaluation of massive amounts of data, resulting in insufficient accuracy in material performance prediction and ineffective convergence of structural relaxation.

[0005] To address the aforementioned problems, this invention provides the following technical solution: This invention provides a method for predicting material properties and high-throughput screening based on a deep learning force field model, comprising the following steps: inputting the set of structures to be calculated into a first-principles calculation module to perform physical quantity calculations, and generating an initial material dataset after unifying physical units and energy reference zero points; performing physical rationality checks on the initial material dataset according to physical screening rules, eliminating abnormal structure samples, and outputting a consistent material dataset; inputting the consistent material dataset into a deep learning force field model based on an isovariant graph neural network for iterative training; constructing a joint loss function using an automatic differentiation mechanism combined with physical conservation constraints and updating model parameters to obtain a trained deep learning force field model; inputting the initial crystal structure of the candidate materials to be screened into the trained deep learning force field model, obtaining the total energy and atomic force vectors, using an unconstrained optimization algorithm to iteratively adjust atomic coordinates for structural relaxation, and outputting a mechanically stable configuration that reaches the mechanical convergence threshold; calculating the multi-dimensional physical properties of the mechanically stable configuration; constructing a multi-index comprehensive scoring function to weight the multi-dimensional physical properties to obtain a final performance evaluation score, and ranking and outputting the candidate materials based on the final performance evaluation score.

[0006] The specific innovative principle of this invention is explained as follows: During the construction of the initial material dataset, structural data from a publicly available material database is extracted, and a virtual structure generation algorithm is invoked. Virtual structures of novel materials are generated through element substitution, random perturbation of unit cell parameters, and the construction of vacancies and defects. These are then aggregated to obtain a set of structures to be calculated. This set of structures is then input into a first-principles calculation module to perform iterative electronic self-consistent field calculations, obtaining the total system energy, the force vectors of all atoms, and the stress tensor. To achieve consistency of multi-source data on a thermodynamic baseline, the total system energy is subtracted from the sum of the products of the reference baseline energy and the corresponding number of atoms of all chemical elements in the isolated single-atom state. This yields the standardized system energy. The virtual structures of the novel materials and associated physical quantities are then stored together after unifying their units.

[0007] In the data physical plausibility testing phase, for each sample in the initial material dataset, the absolute value of the energy difference between two adjacent iterations and the maximum force norm of all atoms in the system are judged to meet the convergence requirements. By calculating the shortest Euclidean distance between atomic pairs in the system, and combining the correlation mechanism between the anti-collapse lower limit coefficient, the anti-dispersion upper limit coefficient, and the empirical covalent radius of two atoms, samples with non-physical overlap or structural dispersion are identified and removed. Furthermore, the imaginary frequency branches on the high-symmetry path of the phonon spectrum are detected by combining the phonon dispersion relation, and imaginary frequency samples are removed. At the same time, the interquartile range method is used to remove high-energy or high-force samples that exceed 1.5 times the normal distribution range, ensuring that the underlying data input to the model conforms to the microscopic physical laws.

[0008] In the deep learning force field model training phase, spatial topological connections between atoms are established based on a preset cutoff radius, constructing a spatial adjacency graph representation of the material configuration. A multilayer perceptron is used to introduce chemical prior feature vectors containing atomic number, electronegativity, covalent radius, and number of valence electrons to perform node feature initialization. A multi-channel Gaussian radial basis function combined with a cosine envelope function is used to construct initial edge feature vectors with smooth boundary transitions. In the forward inference of the network hidden layers, a message passing mechanism is constructed to calculate scalar interaction messages between the center node, neighboring nodes, and spatial distance features. The three-dimensional spatial coordinate vectors of atoms are dynamically updated by multiplying one-dimensional scalar weights with the relative displacement vectors of atoms, maintaining the equivariance of translation and rotation of the Euclidean group mathematically. Finally, a summation function is used to perform an aggregation operation that satisfies permutation invariance. Combined with an energy readout multilayer perceptron, the local energy contribution values ​​of atoms are obtained, and a scalar summation is performed to predict the total energy of the system.

[0009] To ensure the model output conforms to conservative field laws, an automatic differentiation mechanism is used to obtain the negative values ​​of the partial derivatives of the total energy of the material system with respect to the spatial coordinates of each atom in the energy backpropagation path, thus analytically deriving the predicted atomic force matrix. The chain rule is then used to calculate the derivative of the total energy with respect to the strain tensor, analytically deriving the unit cell stress tensor. Subsequently, the predicted total energy, force matrix, and stress tensor are integrated to construct a joint mean squared error loss function containing energy and spatial gradient error terms. During the training phase, various weight parameters are dynamically adjusted, gradient clipping is performed to limit gradient overshoot, and an adaptive moment estimation optimization algorithm is invoked to iteratively update the model parameters.

[0010] During the structural relaxation stage of candidate materials, the stress tensor and the cell strain gradient transformed based on the inverse matrix of lattice basis vectors are integrated into a high-dimensional global gradient vector. Using the gradient difference and position difference vectors from historical iterations, a quasi-Newton method with limited memory is employed to approximate the inverse Hessian matrix of the objective function to construct the search direction. An inexact line search algorithm satisfying the strong Wolf criterion is used to determine the evolution step size to update the state vector. A mechanically stable configuration is output when the maximum value of the Euclidean norm of the atomic three-dimensional force vector, the maximum value of the absolute value of the elements of the cell stress tensor matrix, and the absolute difference of the total energy of the system simultaneously satisfy the multidimensional convergence threshold.

[0011] Before calculating the physical properties, a reference element chemical potential is introduced to calculate the average atomic generation energy for thermodynamic screening. Interatomic spacing and cell volume are verified for geometric stability screening. The dynamic matrix is ​​solved using the model's internal automatic differentiation mechanism, and eigenvalue decomposition is performed. Eigenvalue verification is then used to perform dynamic stability screening. For configurations that pass the initial screening, a three-dimensional deformation tensor is applied to obtain perturbation configurations. The first derivative of stress with respect to strain is predicted to construct the elastic stiffness tensor matrix. The Voigt-Royce-Hill approximation model is used to derive the macroscopic bulk modulus of the polycrystalline material, and the phonon vibrational frequency characteristics in the Brillouin zone of reciprocal space are extracted. The specific heat capacity at constant volume is calculated using Bose-Einstein statistical distribution theory.

[0012] In the final performance evaluation stage, a range standardization operation is performed on the positive and negative indicators of the multi-dimensional physical properties, mapping both positive and negative indicators to values ​​within a defined dimensionless scoring range. A comprehensive score is obtained by linearly summing the preset weighting coefficients and the dimensionless scores, and a Pareto front screening is performed to extract the set of configurations that are not comprehensively superior to other materials across all evaluation indicators. Based on a descending order of the final performance evaluation scores, the tensed lattice basis vectors and coordinates are converted into standard crystallographic information files and written to local storage, completing the screening and output of material structures.

[0013] This invention provides a high-throughput screening method for material property prediction based on a deep learning force field model. It offers the following advantages: 1. This invention constructs a deep learning force field model based on an equivariant graph neural network and utilizes an automatic differentiation mechanism combined with physical conservation constraints to build a joint loss function during model training. This feature enables the model to maintain the equivariance of spatial geometry in translation and rotation during derivation and to analytically determine the atomic forces and unit cell stresses that conform to conservative field laws. Compared to traditional empirical force fields and costly first-principles calculations, this method reduces the computational cost of evaluating massive amounts of data while ensuring consistency of macroscopic physical laws and accuracy of microscopic mechanics, and solves the problem of material relaxation divergence under extreme structural distortions.

[0014] 2. This invention introduces multiple physical plausibility detection mechanisms in the underlying data processing stage. It determines non-physical overlap by calculating the shortest Euclidean distance between atoms and eliminates imaginary frequencies and anomalous samples by combining phonon dispersion relations and interquartile range methods. Utilizing the constraints of these physical rules, the system can proactively intercept divergent data generated by large-scale unknown element combinations before training, avoiding interference from anomalous noise on neural network parameter updates, thereby improving the predictive stability and generalization ability of the force field model when facing unknown crystal configurations.

[0015] 3. This invention designs a multi-dimensional physical property calculation and multi-index comprehensive evaluation system. After the candidate materials complete structural relaxation, thermodynamic, geometric, and kinetic stability screening are performed sequentially, and the bulk modulus and specific heat capacity at constant volume of the materials are further derived. Combined with range standardization operations that distinguish between positive and negative indices and Pareto front screening technology, this evaluation mechanism avoids the limitations of relying on a single performance index for screening. It can directly output usable material configurations that achieve a comprehensive balance among different physical properties, improving the engineering practical value of high-throughput material screening results. Attached Figure Description

[0016] Figure 1 This is a flowchart of the present invention; Figure 2 This is a schematic diagram of the material dataset construction and data quality verification process of the present invention; Figure 3 This is a graph showing the performance index changes during the training process of the deep learning force field model of the present invention. Figure 4 This is a schematic diagram showing the convergence curves of the maximum atomic force during the structural relaxation process, comparing the method of this invention with the traditional empirical force field optimization method. Detailed Implementation

[0017] The technical solutions in 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 embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0018] See attached document Figure 1 This invention provides a method for predicting and high-throughput screening of material properties based on a deep learning force field model, comprising the following steps: S1. Constructing a multi-source materials dataset: The system acquires material structure and energy data from public material databases, and simultaneously calls a structure generation algorithm to construct virtual structures of novel materials. The extracted and generated material structures are uniformly input into the first-principles calculation module to perform physical quantity calculations, obtaining the total energy, atomic forces, and stress tensor data for each structure. After performing unified processing of physical units and energy reference zero points, an initial materials dataset is generated.

[0019] S2 performs physical calculations to verify and assess the data quality of the constructed dataset. Based on preset physical screening rules, the system performs threshold comparison analysis on the structural topology, energy convergence accuracy, and internal force conditions after relaxation of the data samples. The system removes abnormal structural samples that fail the physical plausibility test, retains samples that conform to physical laws, and outputs a consistent material dataset, which serves as the basic data source for training the deep learning force field model.

[0020] S3. A deep learning force field model based on an isotropic graph neural network is constructed and trained. The system divides the consistent material dataset into training and testing sets, which are then input into the deep learning force field model for iterative training. During the model training phase, the system transforms the three-dimensional atomic structure of the material into graph structure data containing node and edge features, and embeds prior knowledge of chemical elements and spatial isotropic rules into the network architecture. The system uses an automatic differentiation mechanism to calculate the gradient of the total energy with respect to atomic coordinates to obtain the predicted force, and constructs a joint loss function by combining multiple physical conservation constraints. The system iteratively updates the model parameters through a backpropagation algorithm until the energy mean square error and mechanical mean square error of the model on the validation set are lower than the preset convergence threshold.

[0021] S4 utilizes a trained deep learning force field model to perform structural relaxation and screening of candidate materials. The system acquires the initial crystal structures of a large number of candidate materials to be screened and inputs them in batches into the trained deep learning force field model. Based on an internally fixed potential energy function, the model outputs the total energy of the current atomic configuration of each candidate material and the force vector of each atom. According to the force directions output by the model, the system uses an unconstrained optimization algorithm to iteratively adjust the atomic coordinates within the system. During the iterative optimization process, when the absolute value of the maximum force on all atoms in the material structure system is less than the set mechanical convergence threshold, the system determines that the structure has reached the relaxation endpoint, terminates the iteration, and outputs a mechanically stable configuration.

[0022] S5 evaluates and sorts the selected materials based on their multi-physical performance. After obtaining the mechanically stable configuration, the system calculates the electronic structure characteristics, elastic modulus parameters, and thermodynamic stability indices of each candidate material through the feature calculation module. The system constructs a multi-index comprehensive scoring function to perform weighted calculations on the acquired multi-dimensional physical properties, resulting in the final performance evaluation score for each candidate material. Based on this performance evaluation score, the system sorts the candidate material set and outputs a sorted list of selected materials.

[0023] In this embodiment, to acquire multi-source data, the system executes a data acquisition and generation process, which specifically includes the following sub-steps: S101, Obtain structural data from publicly available materials databases. The system establishes a data acquisition connection with external publicly available materials databases, extracting the initial crystal structure files of the materials and the corresponding raw data of the system's total energy. Preferably, the publicly available materials databases include the MaterialsProject database and subsets of OC20 and OC22 in the OpenCatalystProject database, extracting the structure covering bulk crystals, surface structures, and adsorption structures. Simultaneously, the system parses structural appendices from published first-principles calculation literature, supplementing them to the initial basic dataset.

[0024] To expand the data range, in S102, a new structure is constructed using a virtual material generator.

[0025] The system reads the baseline crystal structure from the initial basic dataset and calls a virtual structure generation algorithm to generate new material structure samples. The lower-level features of the generation operation include element substitution, random perturbation of unit cell parameters, vacancy and defect construction, and high-throughput structure enumeration. When performing element substitution, the system replaces specific atoms in the baseline crystal structure with elements from the same group according to the principle of consistent valence electron configurations in the periodic table, to maintain the stability of the local chemical bonding properties of the system. When performing random perturbation of unit cell parameters, the system applies a deformation matrix to the three lattice vectors of the baseline crystal to simulate the lattice distortion state of the crystal under stress. The system performs a non-zero determinant check on the applied deformation matrix to prevent the lattice volume from degenerating to zero, which would lead to physical errors. Furthermore, the system removes node atoms from the baseline crystal structure according to a set ratio to generate vacancies, or shifts atoms to interstitial positions to generate point defect structures. Finally, the virtual structures generated by the above operations are summarized to form the set of structures to be calculated.

[0026] After generating the virtual structure, in step S103, first-principles physics calculations are performed. The system batch-inputs the set of structures to be calculated into density functional theory-based calculation software to perform iterative electronic self-consistent field calculations. The calculation software uses the VASP program. The system obtains the physical quantity output results after convergence for each input structure, including the total energy of the system, the force vectors of all atoms in the system, and the stress tensor of the structure. These physical quantities are directly used as the benchmark label data for subsequent neural network model training. For the selection of exchange-correlation functionals, Brillouin zone K-point grid division, and plane wave cutoff energy settings in density functional theory calculations, those skilled in the art can make conventional adjustments based on the specific material system; these are well-known techniques in the field and will not be elaborated upon here.

[0027] Because the total number of atoms and the types of elements vary across different material systems, directly comparing absolute total energy lacks physical meaning. S104, data unification and standardization are performed. Based on the principles of energy conservation and thermodynamic reference states, the system receives the physical quantity results from first-principles calculations, defines the atomic configuration of the material system as a set including atomic numbers and three-dimensional spatial position vectors, and performs unified calculations of the zero-point energy reference. The system subtracts the reference energy of each isolated atom constituting the system to obtain the formation energy reflecting the thermodynamic stability of the material structure. The formula for calculating the standardized total energy is as follows: ; In the formula, This represents the system energy after standardization. This represents the total energy of the system directly output through first-principles calculations. This represents the total number of chemical elements contained in the material system, and its value is a positive integer. The first in the system The number of atoms of a certain element. Indicates the first The reference energy of an element in its isolated monatomic state.

[0028] Reference energy Determine by one of the following two methods: [The text abruptly ends here, likely due to an incomplete sentence or a Individual atoms of each element are placed in a supercell with a preset vacuum layer thickness and obtained through independent density functional theory calculations under the same computational parameters as the original system; or by performing multiple linear regression fitting on the total calculated energy of a large number of known standard compounds. The system embeds non-empty configuration verification logic during the calculation process to avoid errors in the total number of chemical element types. or number of atoms The calculation anomaly is caused by a value of zero.

[0029] After processing, the system will output the configuration and its corresponding normalized system energy. The atomic force vectors and stress tensors are associated and stored in the database to generate an initial material dataset. As a preferred approach, the system uniformly converts the geometric distance unit to angstroms and the energy unit to electron volts, ensuring that the input feature values ​​received by the deep learning force field model in subsequent training stages are within the same dimension and order of magnitude, thus avoiding gradient anomalies in the network.

[0030] See attached document Figure 2In this embodiment, after the system completes the construction of the initial material dataset and the calculation of physical quantities, it enters the physical consistency verification and data quality assessment stage. The purpose of the assessment is to eliminate abnormal samples caused by the non-convergence of the electronic self-consistent field or non-physical deformation, so as to avoid introducing noise into subsequent model training. The system performs a traversal detection of the dataset based on multi-dimensional physical constraint criteria, specifically including the following sub-steps: S201 performs energy convergence accuracy and internal force assessment. Based on the iterative optimization principle of density functional theory, when the change in the total energy of the system is below a set threshold and the forces on each atom approach zero, it indicates that the material configuration is in a local minimum region of the potential energy surface. Based on this principle, the system analyzes the iteration log output by the calculation software and extracts the final step energy change and the atomic force matrix after structural convergence. The assessment logic is executed by the following formula: ; ; In the formula, and Representing the first The total energy of the system output from the previous and next electronic self-consistent field iteration steps. This represents the absolute value of the energy difference between two adjacent iterations. This is the set threshold for energy convergence accuracy. This represents the maximum scalar force value experienced by all atoms within the current structural system. It is a positive integer representing the total number of atoms greater than zero in the system. , , Corresponding to the number in sequence The force components of an atom in three orthogonal directions in a three-dimensional Cartesian coordinate system. This is the force convergence threshold. As a preferred approach, the system will... Set to 10 -5 Electron volts The value is set to 0.02 electron volts per angstrom. This parameter range balances the cost of high-throughput computation with the accuracy of the dataset, ensuring that there are no residual stresses within the preserved configuration that would cause spontaneous deformation.

[0031] S202, Perform structural topological stability testing. To prevent non-physical atomic overlap or chemical bond breakage during relaxation, the system reads the crystal coordinate files before and after calculation, handles the periodic boundary conditions of the crystal using the minimum mirror convention, and calculates the shortest Euclidean distance between all atomic pairs within the system. To avoid diagonal singularities in matrix operations, the system masks zero-distance comparisons of atoms with the same index during traversal. The system's criteria for determining whether the material has experienced structural collapse or discrepancies are as follows: ; In the formula, This indicates that, after incorporating periodic boundary conditions, the first... The and the first The shortest distance between atoms in three-dimensional space, where . and These represent the empirical covalent radii of the two atoms, and their values ​​are obtained by consulting the covalent radius database of the standard periodic table. and These are the lower limit coefficient for preventing landslides and the upper limit coefficient for preventing dispersion, respectively. The system iterates through the above distance data; if any exists... If the distance is less than the lower limit, the atoms are determined to be non-physically overlapping; if the shortest distance of an atom to all its nearest neighbor atoms is greater than the upper limit, the structure is determined to be discrete. Values ​​ranging from 0.5 to 0.7 are used to accommodate bond length shortening caused by high pressure or localized strain. The value is between 1.3 and 1.5 to accommodate the bond length stretching caused by thermal vibration.

[0032] S203 performs dynamic and thermodynamic stability verification. For samples that pass the initial force and topological screening, the system checks whether they are in an unstable saddle point state on the potential energy surface. The system solves for the second-order force constant matrix based on the finite displacement method or density functional perturbation theory, obtains the phonon dispersion relation, and detects whether there are imaginary frequency branches with frequencies less than zero on the high-symmetry path of the phonon spectrum. Due to the numerical calculation error introduced by the computational grid truncation effect and Fourier interpolation at the center of the Brillouin zone, the system introduces a tolerance judgment mechanism. The system excludes unstable structures with imaginary frequency absolute values ​​greater than 0.5 terahertz, retaining samples in a thermodynamically metastable state.

[0033] S204, Remove outlier data and generate a consistent dataset. The system summarizes the detection results of each sub-step and removes structures that fail to simultaneously meet the above three physical criteria from the data stream. To prevent the single extreme value cleaning rule from mistakenly deleting high-performance metastable phases, the system simultaneously performs statistical analysis on the average formation energy and maximum force distribution of single atoms in the global dataset. In this embodiment, the system uses the interquartile range method to calculate the first and third quartiles of the average formation energy of single atoms for each sample, and removes high-energy or high-force samples that exceed 1.5 times the normal distribution range. The system serializes and packages the retained structures and their associated physical property data, outputting a consistent material dataset as the underlying constraint for training the deep learning force field model.

[0034] In this embodiment, after acquiring a consistent material dataset, the system performs the construction and training of a deep learning force field model based on an isotropic graph neural network. The foundation for building this model lies in transforming the three-dimensional continuous spatial configuration of the material into a discrete mathematical expression that the network can process, and incorporating prior physicochemical knowledge. This process specifically includes the following sub-steps: S301, constructing a spatial adjacency graph representation of the material configuration. Based on the locality principle of atomic interactions in quantum mechanics, the system reads the crystal structure information and maps it to an undirected graph model. In this graph model, nodes correspond to the atoms within the system. The system establishes connections between nodes by calculating the shortest spatial distance between atoms after combining periodic boundary conditions and lattice basis vector tensors. To balance computational complexity with the interaction range, the system constructs an interaction edge between two nodes only when the spatial distance between two atoms is less than or equal to a preset cutoff radius. During this process, the system excludes connections between atoms and themselves to prevent anomalies where the distance is zero in subsequent calculations. As a preferred approach, the cutoff radius is set to a range of 5.0 Å to 8.0 Å to cover short-to-medium-range many-body interactions in the first to third coordination layers. For the algorithm implementation of constructing a neighbor list based on periodic boundary conditions and the cutoff radius, those skilled in the art can use conventional cell list methods or generalized search tree methods. Spatial nearest neighbor search is a well-known technique in the field and will not be elaborated here.

[0035] After establishing spatial topological connections, the network needs to obtain the initial feature inputs of each node. S302, perform chemical prior embedding of node features.

[0036] The system uses a multilayer perceptron to map discrete chemical element information into high-dimensional continuous features, injecting the physicochemical properties of the elements into the network. The initial embedding logic of the nodes is executed by the following formula: ; In the formula, Indicates the first The initial node feature vector of each atom, with a dimension equal to the preset hidden layer dimension. Typically, it is set to 128 dimensions or 256 dimensions to ensure the capacity for feature representation. This represents the learnable weight matrix of the embedding layer. Indicates the first The one-hot encoded vector corresponding to the atomic number of each atom. Indicates the first The physicochemical a priori eigenvectors of each atom contain the electronegativity, covalent radius, and number of valence electrons of the corresponding element retrieved from a standard library. Electronegativity characterizes the charge transfer tendency within the system, covalent radius determines steric hindrance, and the number of valence electrons dominates bond strength. (Symbols) This represents the concatenation operation of feature vectors. This corresponds to the learnable bias vector. This operation associates physical laws at the starting point of the model's data flow, improving the model's ability to generalize to combinations of unseen elements.

[0037] After node features are initialized, edge features are constructed. A single scalar distance is insufficient to express the physical law of nonlinear decay of potential energy.

[0038] S303, Perform spatial distance radial basis expansion of edge features. The system uses a multi-channel Gaussian radial basis function to transform the distances between atoms into high-dimensional vectors, and applies a cosine envelope function to ensure that the distance features smoothly transition to zero at the truncation radius boundary, avoiding abrupt changes in the energy function at the boundary. The expansion process is shown in the following formula: ; ; In the formula, Indicates the first The and the first After the atomic spacing is expanded, the first... Gaussian channel eigenvalues. It represents the spatial distance between two atoms. This is the preset total number of channels, usually set to 50. These are the Gaussian channel eigenvalues, ranging from 1 to... Positive integers. Indicates the first The center position of the Gaussian channel eigenvalues, which is between 0 and the cutoff radius. The interval is uniformly and equidistantly distributed. This represents the width control coefficient of the function, whose value is determined by the center interval of adjacent basis functions, and is set to . This is to ensure that the features of each channel overlap in space. This represents the cosine envelope function. The cutoff radius is [value]. Pi is a constant. The system will verify this before calculation to ensure... The value is greater than zero to avoid computational errors caused by a zero denominator in division operations. Through the above expansion process, the system combines the channel values ​​into an initial edge feature vector, providing high-dimensional spatial geometric input for subsequent graph convolutional layers.

[0039] In this embodiment, after acquiring initial node and edge features containing prior information, the system extracts multi-body atom interaction features spanning multiple coordination layers through graph convolution operations. In traditional deep learning frameworks, models struggle to adapt to the physical symmetry of three-dimensional space. To ensure that the predicted physical properties follow the translational and rotational equivariance of the three-dimensional Euclidean group, the system performs message passing and energy readout operations with physical consistency constraints. This process specifically includes the following sub-steps: S304, execute spatial message passing that preserves the equivariance of the Euclidean group. In microscopic physical systems, the scalar potential energy of the system does not change with the translation or overall rotation of the coordinate system (invariance), while the force vectors of atoms undergo synchronous deflection with the overall rotation of the coordinate system (preserving equivariance). Based on the above physical principles, this embodiment adopts an equivariant graph neural network architecture, independently constructing a message passing mechanism in each hidden layer of the network. This mechanism uses a multilayer perceptron to integrate the features of the central node, neighboring nodes, and spatial distance to calculate the scalar interaction messages between atoms, thereby dynamically updating the spatial coordinate representation of the atoms. The specific update logic is executed by the following formula: ; ; In the formula, Indicates the first From atoms in layered networks Transfer to atoms The interaction message vector, where This is the hierarchical index of the network hidden layer, with values ​​ranging from 0 to... integers, This represents the total number of network layers. Indicates the first The message computation layer of the multilayer perceptron contains a linear mapping layer and a nonlinear activation function. and Each represents an atom With atoms In the Layer dimension is eigenvectors. This represents the constant edge eigenvector after radial basis expansion. and Each represents an atom In the Layer and First The updated 3D spatial coordinate vector of the layer. Represents atoms The set of spatial neighbor nodes. This represents the L2 norm of a vector, which is the Euclidean distance between two atoms. This means mapping high-dimensional messages to coordinates of one-dimensional scalar weights to update the multilayer perceptron. Represents a positive real constant. As a preferred method, The value is 10 -8The technical purpose of introducing this parameter is to prevent the norm denominator from becoming zero when the coordinates of two atoms coincide due to numerical truncation errors, thus avoiding singularities in gradient calculations. By multiplying the scalar weights by the relative displacement vectors of the atoms, this coordinate update operation mathematically ensures the rotational equivariance of the spatial vectors, enabling the model to adaptively capture the dynamic geometric evolution of the material during relaxation.

[0040] S305 performs node feature aggregation and updating for multi-body interactions. After calculating the interaction messages between neighboring atoms, the network aggregates these messages to simulate the physical effect of local multi-body potential energy. The system uses a summation function to perform a permutation-invariant aggregation operation on the messages passed from neighboring nodes, and combines this with the features from the previous layer of the central node to output the updated high-dimensional node features. The feature update process is as follows: ; In the formula, This represents the aggregation operation on all incoming messages within the set of neighboring nodes. This addition operation satisfies the permutation invariance of the input node order, which is consistent with the fact that the same type of atoms are physically equivalent in crystal materials. Represents atoms In the The new feature vector output by the layer. Indicates the first The nodes of each layer update the multilayer perceptron. The system uses stacking... The above-mentioned convolutional modules enable each atom to perceive its coordination environment at a greater distance. (Number of network layers) The settings are 3 to 6 layers to obtain a sufficient receptive field without causing feature oversmoothing.

[0041] S306 executes the closed loop for local energy readout of atoms, prediction of the total system energy, and model training. After information transfer via LL layer isomorphic graph convolution, the high-dimensional feature vectors of each node have encoded their local chemical environment and geometric configuration. To achieve the final prediction, the hidden features need to be converted into macroscopic physical quantities. Based on the physical assumption in quantum mechanics that the total system energy is approximately decomposed into the sum of the local energy contributions of each atom, the system is configured with a readout module that outputs the energy contribution of each atom. The readout calculation logic is as follows: ; In the formula, The first number obtained from the prediction The local potential energy contribution of each atom. This indicates an energy readout multilayer perceptron module, which will have a dimension of final layer node features The mapping is a one-dimensional scalar. The system contributes its energy to the local energy of all atoms within the system. Perform a scalar summation operation to obtain the total energy of the material system predicted by the neural network. The technical objective of this summation design is to comply with the physical requirement in thermodynamics that energy is an extensive property, ensuring that the total energy is proportional to the number of atoms.

[0042] To construct a complete closed-loop training process for the model and ensure that those skilled in the art can reproduce it, the system not only outputs the total energy but also uses automatic differentiation techniques to calculate the total energy of the material system. The negative values ​​of the partial derivatives with respect to the underlying input coordinate matrix XX are used to obtain the predicted atomic force matrix. During the model training phase, the system extracts the true energy labels from the aforementioned physically consistent dataset. With force label A joint loss function is constructed. The relevant calculation and evaluation constraints are as follows: ; ; In the formula, For the joint loss function, and These represent the hyperparameter weights for energy loss and force loss, respectively. As a preferred approach, since energy is a global scalar and force is a local vector, the system will initially use... The value is set to be greater than The values ​​are prioritized to stabilize the global baseline of the potential energy surface; as the training algebra progresses, the system dynamically adjusts the weights, gradually increasing... The proportion of the loss function is used to drive the model to fit the local gradient of the potential energy surface (i.e., the atomic forces). The system updates the network parameters of each layer through backpropagation algorithm and optimizer until the joint loss function converges to the preset error range, and outputs the trained physical constraint isotropic graph neural network potential energy model.

[0043] After completing the forward energy prediction and obtaining the total energy of the system, extracting the mechanical state parameters of the material to drive structural evolution is the foundation for establishing molecular dynamics or configuration optimization processes. In this embodiment, the model enters the inference application stage to obtain the complete mechanical state of the material. To ensure that the prediction results comply with the conservative field laws in physics, the system performs energy and force consistency prediction based on the automatic differentiation technique of the computational graph and drives the geometric evolution of the material configuration. This process specifically includes the following sub-steps: S307, Analytical Atomic Forces Based on Automatic Differentiation. In traditional first-principles calculations, force determination relies on the derivation of wave function derivatives, while empirical force fields often introduce truncation errors due to numerical differences. To overcome these shortcomings, this embodiment employs automatic differentiation technology to directly obtain the partial derivatives of energy with respect to the spatial coordinates of each atom along the energy backpropagation path. The calculation logic is as follows: ; In the formula, The model predicts the first The three-dimensional force vector of an atom has a dimension of 1×3. This represents the total energy scalar of the system output by the model's forward propagation. Indicates the first The spatial coordinate vectors of each atom in the Cartesian coordinate system. The purpose of using automatic differential calculation to calculate the force is that, mathematically, this method guarantees that the output force field is a conservative field, meaning the curl of the force vector field is always zero. This physical consistency characteristic avoids energy drift phenomena that violate the laws of thermodynamics in dynamic simulations.

[0044] Obtaining the forces acting on microscopic atoms only guides the relative displacements of internal atoms. For crystalline materials with periodic boundary conditions, the system also needs to evaluate the stress state under macroscopic deformation of the crystal to support the optimization of the cell shape and volume. S308 executes the analysis and derivation of the cell stress tensor. Based on the principles of continuum mechanics, the system constructs coordinate system transformation logic in the computational graph, introduces a symmetric strain tensor, and represents the absolute coordinates of atoms and lattice basis vectors as the deformed shapes after being acted upon by this strain tensor. Through the connection of the above computational graph, the system uses the chain rule to calculate the derivative of the total energy with respect to this strain tensor, thereby obtaining the stress tensor. The formula for calculating the stress tensor is as follows: ; In the formula, The stress tensor of the unit cell is represented by a 3×3 symmetric matrix. This represents the unit cell volume of the current material configuration. Before calculation, the system reads the lattice basis vectors, performs a scalar triple product to calculate the current volume, and sets a volume lower limit verification mechanism. As a preferred method, the system requires... The cubic angstrom operation is used to prevent singularity errors in division operations caused by zero volume due to abnormal model predictions or initial structural collapse. This represents the symmetric strain tensor applied to the material system. By analyzing the stress tensor, the model acquires the physical ability to sense changes in the periodic boundary scale of the lattice, providing a data foundation for predicting the structural response of materials under high-pressure physical environments. After constructing a computational graph that supports energy, force, and stress output, the system enters the model training phase.

[0045] In this embodiment, after completing forward inference and backward automatic differentiation, the deep learning force field model needs to update its network parameters through a supervised learning mechanism. To ensure that the trained model reproduces the thermodynamic and dynamic properties of real physical systems, the system constructs a joint loss function that integrates physical conservation laws to guide the iterative optimization of network weights. This process specifically includes the following sub-steps: S309, Constructing a Joint Loss Function Based on Multi-Physical Quantities. Deep learning-based potential energy surface fitting faces the problem of local curvature prediction distortion; relying solely on energy labels for regression can lead to biases in the predicted direction of atomic forces. To establish physical consistency constraints, the system integrates the total energy of the system, atomic forces, and unit cell stress tensors based on the mathematical differential relationship between scalar potential energy and vector conservative force fields, constructing a weighted mean square error joint loss function. At the data level, the system extracts molecular dynamics trajectories calculated using first-principles methods such as density functional theory. To prevent misalignment of multi-source data from deviating from the optimization direction, the system executes alignment logic based on structural frames, requiring the model's predicted energy, force, and stress tensors to match the true labels under the same time step and the same microscopic geometry. After data alignment, the system executes joint supervision constraints, ensuring the model matches the first-order spatial gradient of the potential energy surface while fitting the potential energy surface function values. To clearly express the composition of each error, the calculation logic of the joint loss function is split into two parts: energy error term and spatial gradient error term, as follows: ; ; In the formula, This represents the first-order spatial gradient loss term of the potential energy surface, which is composed of force and stress. , and These represent the non-negative weighting coefficients for the energy loss term, force loss term, and stress loss term, respectively. This represents the total number of samples in a single training iteration. This represents the index of each independent sample within a batch, with a value ranging from 1 to... Positive integers. and They represent the first The predicted total energy of a sample system and the actual total energy obtained by first-principles calculation. Indicates the first The total number of atoms in each sample system. The system divides the total energy difference by... The technical objective is to eliminate the inconsistency in the absolute energy scale caused by differences in the number of atoms among different samples, and to convert the evaluation standard into the average energy error per atom. Because The number of physical atoms actually existing within the system is always a positive integer greater than or equal to 1, which eliminates the risk of calculation anomalies due to a zero denominator at the underlying logic level.

[0046] For the spatial gradient loss term Vector and tensor terms in and They represent the first In the nth sample The model predicts the force vector of an individual atom, and the actual force vector is compared. The multiplier 3 in the denominator of the formula corresponds to the three component dimensions in the three-dimensional coordinate system. and They represent the first The predicted unit cell stress tensor and the actual unit cell stress tensor for each sample. This represents the Frobenius norm of the computation matrix, which is the root of the sum of squares of the matrix's elements. The denominator multiplier of 9 corresponds to the total number of elements in the 3×3 stress tensor matrix. Through the coordinated measurement of multiple physical quantities, the system achieves joint error constraints on the microscopic atomic positions and macroscopic unit cell scale.

[0047] S310 performs dynamic weight allocation for multi-dimensional joint loss. Considering the differences in physical dimensions and numerical distribution ranges of energy, force, and stress, using fixed weight parameters can easily lead to the model getting trapped in local optima in the later stages of training. The system dynamically adjusts the weights according to different stages of model training. , and The relative proportions. In the initial stage of model training, the network needs to determine the global baseline profile of the potential energy surface. As a preferred approach, the system sets... The initial value is in the range of 1.0 to 5.0. and The range is set to 0.01 to 0.1, with global energy convergence as the primary consideration. As the training cycle progresses, the system triggers the weight scheduler based on the proportion of the current training round to the total number of rounds, gradually reducing the weight according to a preset exponential or linear function. The proportion, and simultaneously increase and The proportion. At the end of training, the system will... and The ratio is dynamically adjusted to the range of 1000:1 to 5000:1. This dynamic adjustment mechanism forces the model to focus on fitting the microscopic local curvature of the potential energy surface in the later stages, ensuring the accuracy of the atomic force gradient.

[0048] S311, Perform adaptive parameter optimization based on gradient truncation. After determining the value of the loss function and the weight configuration, the system uses the backpropagation algorithm to calculate the scalar output value of the joint mean squared error loss function. The system calculates the gradients of all learnable weight matrices and bias vectors within the network. To prevent uncontrolled parameter update step sizes and matrix calculation divergence due to gradient overshoot in deep, equivariant networks, the system performs global norm truncation before gradient updates. The system calculates the L2 norm of the global gradient in real time. If this norm exceeds a preset truncation threshold, the system scales the gradient proportionally to the threshold range to ensure a smooth optimization trajectory. As a preferred approach, this truncation threshold is set between 1.0 and 5.0, determined based on the empirical ratio of the feature dimension of the hidden layers to the number of layers. After truncation, the system calls the adaptive moment estimation optimization algorithm to iteratively update the model parameters, setting the base learning rate to 10 during this process. -4 Up to 10 -3 To avoid misjudgments due to relying solely on a single indicator or the error of a single iteration, the system establishes a multi-dimensional early stopping evaluation mechanism on an independent validation set. The system continuously monitors the mean absolute error of energy and the mean absolute error of force on the validation set until both are below the preset tolerance for multiple consecutive rounds, at which point the model is deemed to have reached convergence and training is terminated. The specific configuration of the learning rate cosine annealing decay strategy, momentum parameter settings, and nonlinear activation functions in the optimizer can be set by those skilled in the art according to conventional machine learning principles; hyperparameter tuning is a well-known technique in the field and will not be elaborated upon here.

[0049] Reference Appendix Figure 3 , attached Figure 3 This document illustrates the training time, energy training, and training loss curves of the model during training. It details the specific configurations of the learner's learning rate cosine annealing decay strategy, momentum parameter settings, and nonlinear activation function. Specifically, it shows the training time per epoch (in seconds) and its average value, the Energy MAE / Train / Validation (in electron volts) between the training and validation sets as a function of epochs, and the trend of the total loss for both datasets. Those skilled in the art can set these parameters according to conventional machine learning principles; hyperparameter tuning is a well-known technique in the field and will not be elaborated upon here.

[0050] In this embodiment, to obtain the stable ground state configuration of the material under specific thermodynamic conditions, the system executes a gradient optimization algorithm for the geometric structure based on the force and stress tensors output by the trained deep learning force field model. The aforementioned deep force field model employs a graph neural network or an equivariant neural network structure. In the forward data flow of this model, the system extracts the element type, three-dimensional absolute coordinates, and lattice basis vectors of atoms as input features. After message passing and feature aggregation of multiple layers of local physical information, it outputs physical state quantities such as the global energy scalar, microscopic force vector, and macroscopic stress tensor. This optimization process aims to search for local minima on the potential energy surface to eliminate unreasonable atomic spacing or lattice deformation. This process specifically includes the following sub-steps: S401, extracting the initial physical state and global gradient vector. The initial configuration of materials typically deviates from equilibrium, with forces existing between internal atoms. Based on the technical requirement of constructing input parameters for a standard nonlinear optimization problem, the system reduces tensors and vectors with different physical dimensions to the same mathematical space. Specifically, the system inputs the initial atomic coordinates and lattice basis vectors into a pre-trained deep force field model. Through model forward propagation and automatic differentiation, it obtains the total energy of the system, the three-dimensional force vectors of all atoms, and the stress tensor of the unit cell. To perform unified nonlinear optimization, the system integrates the force tensor and the unit cell strain gradient transformed based on the inverse matrix of the lattice basis vectors into a high-dimensional global gradient vector. The system maps the stress tensor to its derivative with respect to the lattice basis vectors using the Jacobian matrix, avoiding physical distortion caused by directly superimposing stress onto crystal coordinates. The specific integration logic is as follows: ; In the formula, Indicates the first The global gradient vector during the second optimization iteration has a dimension of (3). +9)×1. Superscript This represents the current iteration step number, and its value is a non-negative integer. Indicates the first In the nth iteration The three-dimensional force vector of an atom, with the negative sign indicating the direction of the gradient along which the potential energy decreases. Indicates the first The 3×3 stress tensor of the next iteration. Indicates the first The unit cell volume at the next iteration. The system multiplies the stress tensor by the volume, the technical purpose of which is to unify the dimensions of a single physical quantity, so that both the macroscopic stress gradient and the microscopic atomic force gradient possess the mathematical characteristics of the first derivative of energy with respect to spatial distance. To ensure dimensional consistency between matrix and vector concatenation, the system performs a tensor flattening operation here, expanding the 3×3 stress tensor into a 1×9 row vector according to row-major order, and then transposing it and incorporating it into the global gradient column vector, thereby allowing linear addition and scalar multiplication operations to be performed within the same vector space.

[0051] After constructing the high-dimensional gradient vector, the system needs to determine the evolution path to overcome the oscillation phenomenon on the potential energy surface.

[0052] S402, search direction construction based on the quasi-Newton method. On a multidimensional potential energy surface, directly using the negative gradient direction for updates can easily lead to oscillations in the iteration trajectory. Based on the general principles of second-order optimization algorithms, the system employs a finite-memory quasi-Newton method (L-BFGS algorithm) to construct the search direction, utilizing historical iteration information to approximate the inverse Hessian matrix of the objective function. The technical purpose of this operation is to obtain the curvature information of the potential energy surface to accelerate convergence, while avoiding the computational cost of accurately calculating the second-order derivative matrix. The formula for calculating the search direction is as follows: ; In the formula, Indicates the first The global search direction vector for the next iteration. Indicates that the system utilizes the most recent The approximate inverse Hessian matrix is ​​obtained by updating the historical position difference and gradient difference vector of each step. As a preferred method, the number of steps is memorized. The value is set to an integer between 5 and 20. To prevent computational crashes caused by non-positive definiteness or singularities in the approximate inverse Hessian matrix during iteration, the system uses a double-loop recursive algorithm based on vector inner product to replace the actual matrix inversion operation in the low-level calculation, and sets curvature verification conditions. Specifically, the system calculates the inner product of the historical gradient difference vector and the position difference vector in real time. If the inner product is less than a preset tolerance lower limit (e.g., 10), the system will check the curvature. -8 The system skips the matrix update in this step to ensure the validity of the descent direction and the numerical stability of the matrix operations.

[0053] After establishing the search direction in multidimensional space, the system needs to calculate the evolution step size along that direction to complete the state update of the material's geometric configuration. In step S403, a line search and state update satisfying the Wolf criterion are executed. An excessively large step size can easily lead to atomic overlap or exceeding local minima, while an excessively small step size results in low computational efficiency. The system uses an inexact line search algorithm satisfying the strong Wolf criterion to determine the step size and simultaneously update the material configuration. The state update formula is as follows: ; In the formula, Indicates the first The system state vector at the next iteration encompasses the Cartesian coordinates of all atoms and the lattice basis vectors expanded into one-dimensional vectors in row-major order. This represents the updated system state vector. This represents the scalar step size determined by the line search algorithm. This step size must simultaneously satisfy both the sufficient descent condition and the curvature condition. The physical significance of these two conditions is to ensure that the total energy of the system monotonically decreases after each update, preventing the configuration evolution from falling into a high-energy divergence state.

[0054] After obtaining the new material configuration through state updates, the system verifies whether the current structure has reached a state of mechanical equilibrium through a decision mechanism. S404, Multi-dimensional Equilibrium State Convergence Decision. After each state update, the system determines whether the current structure has reached a state of mechanical equilibrium based on the new state vector. The atomic coordinates and unit cell dimensions are updated, and the new forces and stresses are calculated using the deep force field model. To evaluate whether the material system has reached mechanical equilibrium, a multi-dimensional convergence criterion incorporating both microscopic and macroscopic parameters is established.

[0055] Considering that relying solely on a single energy difference extremum to determine convergence is prone to misjudgment in regions with flat potential energy surfaces, the system introduces a dimensional verification based on force and stress. In this embodiment, the system transforms the extremum comparison logic into a scalar threshold verification. Specific judgment conditions require: the maximum L2 norm of the three-dimensional force vectors of all atoms in the system is less than or equal to a preset force convergence threshold; the maximum absolute value of all elements in the unit cell stress tensor matrix is ​​less than or equal to a preset stress convergence threshold; and the absolute difference in the total energy of the system between two iterations is less than or equal to a preset energy difference convergence threshold.

[0056] The force convergence threshold is set to 0.01 to 0.05 electron volts per angstrom, the stress convergence threshold is set to 0.01 to 0.1 gigapascals, and the energy difference convergence threshold is set to 10. -5 Up to 10 -4 Electron volts. The system is determined to have reached a local minimum on the potential energy surface if and only if the conditions of the above three dimensions are simultaneously satisfied. The optimization process terminates, and the stable configuration is output as the basis data for subsequent analysis. Through this joint determination mechanism, the system avoids the defect of being misjudged as convergent due to a small single step length caused by a flat potential energy surface. For the specific derivation of the double-loop recursion in the quasi-Newton algorithm and the interpolation logic of the line search, those skilled in the art can refer to standard optimization theory for configuration, which is a well-known technique in this field and will not be elaborated here.

[0057] After obtaining material configurations at local minima of the potential energy surface through the aforementioned geometric optimization, not all structures in mathematically convergent states possess physical synthesizability. In this embodiment, to eliminate configurations in high-energy metastable states or with structural defects, the system performs preliminary stability screening on the optimized candidate materials based on multi-dimensional physicochemical criteria. This screening process comprehensively considers the thermodynamic, geometric, and kinetic properties of the system, avoiding a single judgment based solely on the total energy extremum. This process specifically includes the following sub-steps: Considering the physical requirements for assessing the spontaneous formation tendency of materials under natural conditions, the system needs to establish an energy evaluation benchmark with cross-system comparability. S405 calculates the absolute formation energy of the system and performs thermodynamic screening. The single total energy only reflects the benchmark value within a specific computational simulation framework. The system introduces the chemical potential of a reference element to convert the model-predicted total system energy into the formation energy, thereby characterizing the energy change during the combination of elements in the material under standard conditions. The specific calculation formula is as follows: ; In the formula, This represents the average formation energy per atom. This represents the total system energy output by the deep force field model for the current optimized convergent configuration. To ensure the consistency of the energy baseline, the system performs energy zero-point alignment before calculation to guarantee... The reference frame is synchronized with the reference frame of the chemical potential of the element. This indicates the total number of element types contained in the system. Indicates the element type index. Indicates the first in the system The number of atoms of a certain element. Indicates the first The system calculates the reference chemical potential per atom of each element in its standard solid state, obtained from publicly available materials databases or first-principles calculations of equivalent precision. The system then generates the sum of energies divided by the total number of atoms of all elements in the system to standardize the evaluation scale for unit cells of different sizes. Since the total number of atoms is always a positive integer, this division avoids the computational risk of the denominator approaching zero. After the calculation is complete, the system will... The energy is compared with a preset thermodynamic threshold. Preferably, this generation energy threshold is set to 0.0 to 0.1 electron volts per atom. If the material is below the threshold, it is considered to have passed thermodynamic screening; otherwise, it is rejected.

[0058] After completing the thermodynamic stability verification, the system needs to investigate non-physical geometric configurations caused by the generalization blind spots of the deep learning model. S406, geometric stability screening based on interatomic spacing and volume constraints. In the forward inference of the neural network, if the input configuration exceeds the training set distribution boundary, the model easily predicts abnormal regions of the potential energy surface, leading to non-physical collapse or expansion of atoms. The system extracts the unit cell volume of the convergent configuration. Calculate the spatial Euclidean distance between any two atoms using the Cartesian coordinates of all atoms. The specific verification logic is transformed into multi-dimensional scalar threshold judgment. System Requirements Greater than or equal to the sum of the empirical covalent radii of all atoms and the overlap tolerance coefficient The product of, where The preferred value is 0.6 to 0.8. Meanwhile, the unit cell volume... Greater than the minimum volume limit obtained by the accumulation of intrinsic atomic occupants If the system finds any atomic spacing below the threshold or the total unit cell volume below the lower limit during the traversal, it determines that the configuration does not meet the geometric constraints and then performs a rejection operation.

[0059] Even if the configuration satisfies thermodynamic and geometric criteria, it is still necessary to confirm its recovery capability under displacement perturbation, which depends on the local curvature characteristics of the potential energy surface. S407, preliminary screening based on second-order spatial gradient dynamic stability. The material being at a local minimum on the potential energy surface only indicates that the first-order force gradient approaches zero. Based on lattice dynamics principles, if the potential energy of the configuration exhibits a downward-opening parabolic characteristic in a certain displacement direction, the lattice will spontaneously distort. The system utilizes the continuously differentiable activation function within the deep force field model and solves for the Jacobian matrix of forces with respect to coordinates through a backward automatic differentiation mechanism. The formula for calculating the elements of this dynamic matrix is ​​as follows: ; In the formula, Represents the corresponding number in the dynamic matrix Atoms in Axial displacement and the first Atoms in Matrix elements of axial displacement. and These represent the physical static masses of the two atoms, and the system ensures that they are real numbers greater than zero during the preprocessing stage to guarantee the safety of division. Indicates the first Atoms in Predicting the direction of force. Indicates the first Atoms in The coordinate components of the direction. The system constructs a complete 3D model. ×3 After the matrix is ​​generated, to eliminate the asymmetric errors introduced by numerical calculation, matrix symmetry processing is performed, i.e., the mean of the original matrix and its transpose is taken. Subsequently, the system performs eigenvalue decomposition on this symmetric matrix to obtain the eigenvalues ​​of the squares of the corresponding system lattice phonon vibration frequencies. The system verifies whether all eigenvalues ​​are greater than a preset frequency tolerance threshold. As a preferred approach, this threshold is set to the square of -0.01 to -0.05 terahertz. If there are negative eigenvalues ​​below this threshold, it indicates that the lattice has imaginary frequency vibrations, and the system determines that it is dynamically unstable and terminates the candidate path. Through the above-mentioned combined filtering mechanism covering thermal, geometric, and dynamic aspects, the system achieves preliminary screening of effective material configurations. For the specific linear algebra algorithm used in eigenvalue decomposition, those skilled in the art can use a conventional numerical calculation framework configuration, which is well-known in the field and will not be elaborated here.

[0060] After the initial stability screening described above, the system retains material configurations with physical synthesis feasibility. To evaluate the service performance of these candidate materials in practical applications, in this embodiment, the system performs multi-dimensional physical property calculations on the candidate configurations based on a deep force field model and statistical physics theory. This calculation process analyzes the microscopic mechanical response and macroscopic thermodynamic statistical behavior of the materials to output the evaluation indicators required for specific business scenarios. This process specifically includes the following sub-steps: Assessing a material's resistance to external deformation is fundamental to determining its structural reliability. S501 utilizes strain perturbation-based elastic tensor calculations. The elastic tensor reflects the linear response of lattice stress to applied strain. The system applies multiple independent three-dimensional deformation tensors to the lattice basis vectors of convergent configurations, constructing a series of perturbation configurations. In this calculation logic, the system inputs tensor features containing the three-dimensional absolute coordinates of atoms and lattice basis vectors into a pre-trained graph neural network force field model. As a preferred approach, this pre-trained model uses samples and labels provided by first-principles calculations to perform parameter fitting by minimizing the joint mean square error loss function of global energy and atomic forces. During forward propagation within the model, the input features are mapped to high-dimensional node features through an embedding layer. Subsequently, a spatial adjacency graph is constructed based on a preset cutoff radius in multiple message passing layers, performing aggregation and updating of local atomic spacing and bond angle information. After feature extraction, the system outputs a global energy scalar using a global pooling layer. To obtain the macroscopic stress, the system invokes the backpropagation automatic differentiation mechanism of the computation graph to calculate the first derivative of the total energy of the system with respect to the strain tensor, and divides it by the cell volume of the current configuration. This achieves dimensional transformation and outputs the physically conserved cell stress tensor. The derivation formulas for the elastic stiffness tensor elements are as follows: ; In the formula, Represents the elements of the elastic stiffness tensor matrix. Subscript and Following the standard Voigt notation, the values ​​are integers from 1 to 6, corresponding to the three normal strain and three shear strain components in three-dimensional space, respectively. This indicates that the stress tensor output by the model is in Components in direction. Indicates that the system applies to Dimensionless deformation rate in the direction. As a preferred method, the applied deformation rate... The value ranges from -0.02 to 0.02, chosen to ensure that the material deformation remains within the applicability range of the linear Hooke's law, avoiding triggering nonlinear plastic deformation that could lead to the failure of derivative calculations. The system uses the finite difference method to extract the derivative of stress with respect to the deformation rate. Within the underlying algorithm framework, the system sets a perturbation increment. The absolute value is greater than the preset lower limit of machine precision (e.g., 10). -5 To avoid numerical divergence when the denominator approaches zero in difference operations, a complete 6×6 elastic tensor matrix is ​​constructed accordingly.

[0061] After obtaining the elastic tensor in the microscopic direction, it is necessary to derive the macroscopic mechanical properties of polycrystalline materials to match the synthesis environment. S502, Derivation of the macroscopic bulk modulus based on the polycrystalline approximation. The macroscopic mechanical properties of polycrystalline materials are mostly isotropic. The system uses the Voigt-Royce-Hill approximation model to calculate the polycrystalline bulk modulus. The calculation formulas for the polycrystalline modulus boundary are as follows: ; ; In the formula, This represents the upper limit of the Voith bulk modulus assuming uniform strain within the polycrystalline material. This represents the lower limit of the Reuss bulk modulus assuming uniform internal stress in a polycrystalline material. The elements of the elastic compliance tensor matrix are represented by their subscripts. and The definition is consistent with that of the stiffness tensor. Mathematically, this compliance tensor matrix is ​​the inverse of the corresponding elastic stiffness tensor matrix. To prevent singularities in the microscopic elastic tensor matrix from causing inversion operations to crash, the system determines the determinant of the 6×6 stiffness tensor matrix before calculation. If the absolute value of this determinant is lower than a preset numerical tolerance limit (e.g., 10), the system will proceed accordingly. -6 If the system determines that the matrix inversion condition is not met, it terminates the modulus derivation and marks the configuration as mechanically unsteady. The system then performs further checks. The absolute value of the sum of the polynomials in the formula, if less than the lower tolerance limit, is also considered a division singularity and triggers blocking logic. After obtaining the upper and lower limit moduli, the system directly calculates... and The arithmetic mean of the macroscopic bulk modulus This scalar is used to characterize the compressive strength of a material under macroscopic conditions.

[0062] Besides the mechanical properties at room temperature, the thermodynamic response of a material under different temperature conditions also determines the boundaries of its service environment. S503, based on the calculation of specific heat capacity at constant volume using phonon vibration frequencies. Specific heat capacity reflects a material's ability to absorb heat, leading to a temperature increase. The system utilizes the phonon vibration frequency characteristics extracted from the aforementioned dynamic matrix eigenvalues, combined with Bose-Einstein statistical distribution theory, to calculate the specific heat capacity at constant volume at a specific temperature. The formula for calculating this thermal property is as follows: ; In the formula, This indicates the specific heat capacity at constant volume of the material. This represents the wave vector component in reciprocal space. The index represents the phonon dispersive band. This represents the Boltzmann constant. This represents the reduced Planck constant. This indicates the absolute temperature of the target's service environment. This represents the lattice phonon angular frequency corresponding to the wave vector and energy band, and its value is derived from the square root of the aforementioned eigenvalues ​​of the dynamic matrix. To avoid the divergence anomaly of the denominator being zero due to the temperature approaching absolute zero in the formula, the system limits the absolute temperature of the target during the input data verification stage. Greater than or equal to 1 Kelvin. The lattice phonon angular frequency is the frequency at which the lattice phonon frequency is measured when lattice dynamics calculations involve the acoustic branch of the Brillouin zone center (Gamma point) in reciprocal space. The exponent polynomial in the formula will also approach zero, causing it to approach zero simultaneously. To address the computational anomalies caused by this physical limit, the system sets a lower frequency cutoff limit during the traversal accumulation. As a preferred approach, when... For frequencies below 0.01 terahertz, the system replaces this term with the analytical limit value of the corresponding Taylor expansion to maintain the numerical stability of the thermodynamic statistical results. The system performs an ergodic summation on all phonon modes within the Brillouin zone to obtain the macroscopic thermodynamic index. The specific linear algebraic algorithm for the inverse operation of the elastic-compliant tensor matrix, and the mesh generation logic for the high-symmetry point paths in the reciprocal space Brillouin zone, can be configured using conventional materials calculation frameworks by those skilled in the art; these are well-known techniques in the field and will not be elaborated upon here.

[0063] Based on the multidimensional physical property tensor output by the graph neural network in the aforementioned steps, since the indicators of different evaluation dimensions have inherent differences in physical dimensions, numerical range, and optimization direction, directly comparing absolute values ​​is difficult to objectively reflect the comprehensive performance of the material. In this embodiment, to match the needs of specific industrial application scenarios, the system constructs a comprehensive scoring function, mapping multidimensional discrete physical properties to a unified scalar evaluation space, and based on this, performs candidate sequence sorting and outputs material structure information files. This process specifically includes the following sub-steps: To address the data incomparability issue caused by inconsistent dimensions of multidimensional physical properties, the system needs to establish a unified evaluation benchmark. S504 involves dimensionless processing and comprehensive scoring of multidimensional physical properties. Based on the principle of multidimensional data normalization, the physical properties of materials include positive and negative indices. Higher positive index values ​​(e.g., compressive bulk modulus) better reflect the required properties for business needs; conversely, lower negative index values ​​(e.g., generation energy or material density) better reflect the required properties for business needs. The system extracts the physical property values ​​of all candidate materials in the current batch and performs direction-aware range standardization. The formula for dimensionless processing is as follows: ; In the formula, Indicates the first The candidate materials were in the first The dimensionless score for a physical property is constrained to a closed interval between 0 and 1. This represents the predicted value of the original physical properties of the material, specifically characterized as the dimensionality reduction value of the elastic tensor or the energy scalar output by the aforementioned pre-trained graph neural network force field model through forward inference and backward automatic differentiation. and These represent the current batch of all candidate materials in the [number]th [stage]. The maximum and minimum values ​​of the physical properties of the item. This represents a set of positive indicators. This represents a set of negative indicators. (Introduction) The physical reason is that when a batch of materials exhibits consistent performance in a specific property, the maximum and minimum values ​​are equal. In this case, the lower tolerance limit can prevent program divergence caused by a zero denominator in division operations. The technical purpose of this range standardization operation is to transform heterogeneous physical quantities spanning multiple orders of magnitude into the same numerical range to achieve comparability.

[0064] After completing the dimensionless mapping of each dimension, the system, considering the technical preferences of specific business scenarios, constructs a scalar-form comprehensive scoring objective function through a weighted linear combination. The formula for calculating the comprehensive score is as follows: ; In the formula, Indicates the first The overall score of each candidate material. Indicates that the system is targeting the first The weighting coefficients assigned to each physical indicator. These weighting coefficients are external input parameters; as a preferred method, their specific values ​​are pre-calculated and imported by external business systems based on the Analytic Hierarchy Process (AHP) or expert scoring systems. The system verifies all [indicators] during the preprocessing stage. All values ​​are non-negative real numbers, and their sum equals 1. The system assigns weights to specific items based on the physical causal boundaries of the target business conditions. For example, in the aerospace structural component selection scenario, structural components need to remain lightweight while bearing large loads. Therefore, the system assigns a high positive weight to the bulk modulus, which represents stiffness, while assigning a high weight to the material density and classifying it into the negative index set. By employing this weighted logic based on multi-dimensional and physical correlation, the system avoids one-sided judgments caused by relying solely on the extreme values ​​of a single indicator.

[0065] Considering the limitations of using a single weighted polynomial in a non-convex optimization space, which can mask the advantages of materials in certain dimensions, the system, in addition to scalar scoring, introduces a Pareto-dominant multi-objective ranking logic as a subordinate feature to fully support the business needs of multi-objective screening.

[0066] S505, Pareto front screening and standardized document output for candidate material sets. The system performs pairwise comparisons for any two material configurations within a multidimensional dimensionless scoring space. If the material... In all The scores on all evaluation indicators are greater than or equal to those for the material. And at least one of the indicators has a score that is strictly greater than the material score. The system determines the material. Pareto is superior to materials The system iterates through the entire candidate set and extracts the Pareto front set of configurations that are not dominated by any other material by executing a fast non-dominated sorting algorithm. This set represents the theoretically optimal family of solutions under the current multidimensional physical constraints that cannot improve one property without sacrificing another.

[0067] Based on the above multidimensional scoring and frontier screening results, the system needs to convert the tensor data within the deep model into a standard format usable in engineering, in order to achieve compatibility between the algorithm's output data and the format of the upper-level industrial physics simulation software. The system uses a comprehensive scoring system... The system employs a descending order sorting rule, or extracts candidate objects from the Pareto front set, and encapsulates the crystallographic information file. The system then reduces the tensorized lattice basis vectors to the lengths of the three lattice constants. Angle with three unit cells Simultaneously, the Cartesian coordinates of the atoms are converted to fractional coordinates relative to the lattice basis vectors, and combined with the element type identifier, are written to the local storage medium according to the text topology rules of standard crystallographic information files. The output file embeds metadata such as generation energy, bulk modulus, and comprehensive score, which can be accessed by downstream experimental personnel or imported into finite element simulation software. The specific string serialization and concatenation logic and memory flow management used in the underlying file input / output system can be developed and configured using conventional programming framework interfaces by those skilled in the art; these are well-known technologies in the field and will not be elaborated upon here.

[0068] To support the aforementioned material generation and measurement tasks based on graph neural networks and physics engines, this embodiment provides an electronic device. This electronic device provides computational power support for tensor operations through a heterogeneous computing architecture. Specifically, this hardware architecture includes the configuration and coordination logic of a processor, memory, and communication interfaces.

[0069] For the tensor data generated during the forward inference and automatic backward differentiation operations of the model, the electronic device is equipped with at least one processor. As a preferred approach, the processor employs a heterogeneous architecture combining a central processing unit (CPU) and a graphics processing unit (GPU) or a tensor processing unit (TPU). The CPU is responsible for executing control flow instructions, file system read / write operations, and preprocessing of unstructured data, while the GPU handles the message passing layer computations and automatic differentiation of the elastic tensors within the graph neural network. During the batch inference phase of the graph model, to avoid memory overflow, the processor's internal memory management unit dynamically calculates the maximum batch graph size for a single forward propagation based on the current hardware state. The formula for calculating this hardware memory constraint is as follows: ; In the formula, This indicates the upper limit of the number of candidate material configuration diagrams processed in a single parallel operation. The outer bracket symbol indicates a floor operation, and the output result is a non-negative integer. This indicates the total available free video memory of the current graphics processor, in megabytes. This indicates the fixed memory usage required to load the parameters of the pre-trained graph neural network force field model. This represents the underlying Computing Unified Device Architecture (CUDA) context memory reserved by the system. Its value is determined based on the hardware driver version, for example, it is 512 megabytes. This represents the average video memory size occupied by a single material configuration after constructing the spatial adjacency graph. In actual engineering deployments, the average video memory size... The value is pre-estimated by the system based on the average number of atomic nodes, node feature dimensions, and edge index size of the current batch of materials, combined with the floating-point precision of the underlying hardware. The system verifies this value before the underlying code is executed. The value is greater than a preset graph memory threshold (e.g., 10⁻³ × 10⁻³ megabytes) to avoid numerical divergence anomalies caused by the denominator approaching zero during division operations. By introducing the above computational constraints, the system divides the overall inference task into sub-batches adapted to the current hardware computing power, ensuring memory safety during the model inference process.

[0070] Electronic devices include memory that is communicatively connected to a processor. This memory consists of volatile random access memory (RAM) and non-volatile storage media (such as solid-state drives with a non-volatile memory host controller interface specification). The non-volatile storage media is used to persistently store computer program instructions, raw crystallography libraries of materials, and the runtime environment of deep learning frameworks. When the device powers on and starts a task, the processor loads the aforementioned program instructions into the volatile RAM via the system bus for execution. To match the throughput of the graphics processor, as a preferred approach, the volatile memory adopts the high-bandwidth memory (HBM) standard. Node feature matrices containing atomic three-dimensional coordinates and lattice basis vectors, as well as edge index tensors constructed based on a preset physical cutoff radius, reside in the volatile memory for the processor to access during matrix multiplication and graph convolution. By locking data in memory at a physically closer level, the system reduces the time overhead of data retrieval.

[0071] To enable data interaction between the processor, memory, and external environment, electronic devices are equipped with communication interfaces to receive element constraint ratio instructions from external business systems and output selected Pareto front material standard files to peripherals. The processor, memory, and communication interface are electrically connected and data is transferred via the PCIe (Peripheral Component Interconnect) bus. During cross-device data copying (e.g., copying tensors from CPU memory to GPU memory), the system's underlying layer sends instructions to the operating system kernel to perform page-locked memory allocation. The technical purpose of this operation is to prevent virtual memory swapping during this process, ensuring that the multi-source tensor data remains aligned with its physical memory address across time dimensions, preventing dangling pointer errors caused by memory relocation. The specific electrical specifications of the PCIe bus and the underlying thread scheduling mechanism of the operating system can be configured using conventional computer motherboards and kernel architectures, which are well-known technologies in the field and will not be elaborated upon here.

[0072] According to one embodiment of the present invention, a computer-readable storage medium is provided. Based on the aforementioned material generation and measurement tasks using graph neural networks and physics engines, this storage medium is used to persistently store computer program instructions and data dependencies supporting the program's execution. When the aforementioned computer program instructions are read and executed by the processor of an electronic device, the multi-dimensional physical property measurement, material selection, and data output logic of the aforementioned embodiments are implemented.

[0073] In this embodiment, the computer-readable storage medium is a non-volatile physical carrier. Preferably, the storage medium uses a magnetic storage medium, an optical disk, or a solid-state electronic flash memory chip. The storage medium internally has independent storage sectors, used to store the operating system boot code, the dynamic link library of the deep learning inference framework, the parameter file of the pre-trained graph neural network force field model, and the crystallographic standard file serialization module. The parameter file stores the parameters of the embedding layer for processing node features, the message passing layer for aggregating topological information, and the multilayer perceptron network parameters for outputting the prediction tensor. This structured file organization provides a defined addressing path for computing power invocation.

[0074] Based on the principle of data integrity verification, to ensure the numerical fidelity of model parameters and prevent bit flipping caused by hardware bad blocks or level fluctuations in the storage medium, the system pre-installs verification metadata in the header area of ​​the model file on the storage medium. When the processor reads tensor weights from the storage medium into volatile memory via the input / output bus, the system executes the data frame integrity verification logic. The calculation formula is as follows: ; In the formula, This represents the scalar of the cumulative error during a single file loading process. This represents the total number of layers in the neural network model parameter tensor stored in the storage medium, and is a positive integer. and They represent the first The number of rows and columns of the layer weight tensor in the corresponding dimension. This indicates the number of times the data was actually read from the storage medium into memory. The coordinates in the layer weight tensor are The parameter elements. The system iterates through all parameters of the model, parameter elements. The verification range includes the weight matrices of each layer and synchronously covers the corresponding bias vectors. Indicates the relationship with the first The sum of the reference absolute values ​​corresponding to the layer tensors is calculated by the system when the model training is completed and written to the storage medium, and is hard-coded in the file's metadata area. This verification logic pre-calculates the absolute values ​​of the parameter elements actually read to prevent positive and negative value deviations caused by individual bit flips from canceling each other out during accumulation. At the same time, compared with executing a cryptographic hash algorithm, this linear accumulation-based verification mechanism reduces the computational overhead and startup latency of the model during memory loading.

[0075] In the data loading and verification logic, the system compares the cumulative error scalar. With the preset fault tolerance threshold parameter As a preferred method, The value is set based on the truncation error caused by the floating-point bit width of the hardware memory; for example, a value of 10... -6 .when Greater than or equal to If the system determines that the read data frame has been physically damaged or has a transmission error, it will block subsequent forward inference computation of the graph neural network to prevent numerical divergence of the output elastic tensor or energy scalar due to model parameter distortion.

[0076] Once the data loading passes verification, the program instructions schedule and calculate the unified device architecture interface, mapping the graph data structure into concurrent thread blocks adapted to the processor, driving the generation of material structures and the physical engine calculations. Regarding the log-type file allocation table structure used at the bottom layer of the file system, the direct memory access control mechanism of block devices, and the page fault interrupt response logic of the operating system kernel, those skilled in the art can use conventional computer storage specifications and kernel technologies for development and configuration; these are well-known technologies in the field and will not be elaborated upon here.

[0077] In one embodiment, Figure 4 This is a schematic diagram illustrating the changes in the mechanical convergence state of the material system recorded by the system of this invention in response to a specific scenario involving sudden extreme lattice distortion (corresponding to a high-energy potential impact on the structure) during high-throughput virtual material generation. The diagram is a two-dimensional coordinate graph, with the horizontal axis representing the relaxation optimization time (or number of iterations). The vertical axis represents the maximum scalar force value of all atoms in the system. .

[0078] Reference Figure 4 The diagram contains the following elements: a horizontal line segment, marked as the mechanical convergence threshold. This line segment represents a preset compliance threshold for the structure to reach a local minimum on the potential energy surface and a state of mechanical equilibrium, i.e., the aforementioned threshold. (e.g., 0.02 electron volts per angstrom). A dashed line, labeled for conventional empirical force field optimization methods. This curve represents the maximum force in a system not employing the method of this invention. The process of change. Before time T1 (when dealing with a known steady-state crystal), the value of the dashed line is below the mechanical convergence threshold. After time T1, with the introduction of structural distortion and unknown elements, the force value rises rapidly and remains above the convergence threshold for a very long period. Furthermore, due to the lack of real quantum mechanical many-body interaction constraints, the force curve exhibits violent oscillations and numerical divergence, ultimately failing to recover below the threshold. A solid line marks the application of the method of this invention. This curve represents the maximum force in the system using the method of this invention. The curve's variation process. Over the vast majority of the optimization timeline, the curve's value converges quickly and stably to the mechanical convergence threshold. Below. Near time T1, due to the impact of a sudden lattice distortion, the curve's value experiences a brief spike, but the system immediately recovers using the global gradient vector output by automatic differentiation. And combined with the quasi-Newtonian search direction Performing configurational evolution causes the maximum force value to rapidly decrease monotonically and then return to a stable state, eventually stabilizing at the mechanical convergence threshold. Below.

[0079] Figure 4 The data also marks two specific time points: Time point T0: This represents the moment when the system's virtual material generator detects and constructs a high-energy unsteady virtual structure containing novel elemental combinations (e.g., replacing some elements in a sulfide-germanium ore type Li6PS5Cl and applying strain) by calling the random deformation matrix and element replacement logic. The method of this invention starts data parsing immediately after time T0, mapping discrete chemical element information to node features. and perform equivariant message passing. Computational preparation for collaborative optimization decision-making. Time point T1: Represents the extreme distortion structure generated at time T0, which, after being batch-transmitted by the data flow control engine, actually reaches the pre-trained deep learning force field model and triggers the first step of coordinate update calculation. At that moment.

[0080] The figure illustrates that, when addressing the technical problem of structural distortion impact at the same extreme, the traditional empirical force field optimization method (shown by the dashed line) suffers from force divergence and non-convergence (i.e., structural non-compliance) after the impact load arrives (after T1) due to the lack of physical spatial Euclidean group equivariance and accurate potential energy surface curvature sensing capability. In contrast, the method of this invention (shown by the solid line) embeds multi-body physical priors at time T0 and executes a method based on a joint loss function at time T1. Constrained Adaptive Gradient Truncation and L-BFGS State Update This allows the system to reliably and quickly search for the true mechanically stable configuration and the maximum force of the system when the extreme distortion impact load arrives (at time T1) and thereafter. Always maintain at the preset mechanical convergence threshold The following approach solves the technical problem that existing technologies cannot effectively handle high-throughput screening of large-scale, unknown, and suddenly distorted material structures.

Claims

1. A method for material performance prediction and high-throughput screening based on deep learning force field model, characterized in that, Includes the following steps: Input the set of structures to be calculated into the first-principles calculation module to perform physical quantity calculations, and generate an initial material dataset after unifying the physical units and the energy reference zero point. The initial material dataset is subjected to physical rationality detection based on physical screening rules, abnormal structure samples are removed, and a consistent material dataset is output. The consistent material dataset is input into a deep learning force field model based on an isotropic graph neural network for iterative training; By using an automatic differentiation mechanism combined with physical conservation constraints to construct a joint loss function and update the model parameters, a trained deep learning force field model is obtained. The initial crystal structure of the candidate material to be screened is input into the trained deep learning force field model to obtain the total energy and atomic force vector. An unconstrained optimization algorithm is used to iteratively adjust the atomic coordinates to relax the structure and output a mechanically stable configuration that reaches the mechanical convergence threshold. Calculate the multidimensional physical properties of the mechanically stable configuration; A multi-index comprehensive scoring function is constructed to calculate the final performance evaluation score by weighting the multi-dimensional physical properties, and the candidate materials are ranked and output based on the final performance evaluation score.

2. The method of material performance prediction and high-throughput screening based on deep learning force field model according to claim 1, characterized in that, The process of inputting the structure set to be calculated into the first-principles calculation module to perform physical quantity calculations, and generating an initial material dataset after unifying the physical units and energy reference zero point, specifically includes: Extract structural data from publicly available material databases, call virtual structure generation algorithms, and generate virtual structures of novel materials by element substitution, random perturbation of unit cell parameters, and construction of vacancies and defects. The sets of structures to be calculated are then summarized. The set of structures to be calculated is input into the first-principles calculation module to perform electronic self-consistent field iterative calculation to obtain the total energy of the system, the force vectors of all atoms, and the stress tensor. The system energy after standardization is obtained by subtracting the sum of the products of the reference energy of all chemical elements in the isolated single-atom state and the corresponding number of atoms from the total energy of the system. The configuration of the structure set to be calculated, the standardized system energy, the force vectors of all atoms, and the stress tensor are unified and stored together to generate the initial material dataset. 3.The method of predicting and high-throughput screening of material performance based on deep learning force field model according to claim 1, wherein, The steps of performing physical rationality checks on the initial material dataset according to physical screening rules, removing abnormal structure samples, and outputting a consistent material dataset specifically include: For each sample in the initial material dataset, determine whether the absolute value of the energy difference between two adjacent iterations is less than or equal to the set energy convergence accuracy threshold, and determine whether the maximum force scalar value of all atoms in the system is less than or equal to the set force convergence threshold. If both are less than or equal to, it is determined that the physical convergence condition is met and the corresponding sample is retained. If there is a case where it is greater than, it is determined that the physical convergence condition is not met and the sample is removed. The energy convergence accuracy threshold is set based on the default electronic step convergence accuracy of the first-principles calculation module. The set force convergence threshold is set based on the negligible reference noise level of the interatomic forces in the material. Calculate the shortest Euclidean distance between all atomic pairs in the system. If the shortest Euclidean distance is less than the product of the lower limit coefficient for preventing collapse and the sum of the empirical covalent radii of the two atoms, non-physical overlap is determined and the pair is eliminated. If the distance is greater than the product of the upper limit coefficient for preventing dispersion and the sum of the empirical covalent radii of the two atoms, structural dispersion is determined and the pair is eliminated. By combining the phonon dispersion relation, it is determined whether there are imaginary frequency branches with frequencies less than zero on the high symmetry path of the phonon spectrum. If they exist, they are judged as abnormal and the corresponding samples are removed. If they do not exist, they are judged as normal and the corresponding samples are retained. High-energy or high-force samples that exceed 1.5 times the interquartile range of the normal distribution range are directly removed using the interquartile range method, and the consistent material dataset is output.

4. The method of material performance prediction and high-throughput screening based on deep learning force field model according to claim 3, characterized in that, The steps of inputting the consistent material dataset into a deep learning force field model based on an isotropic graph neural network for iterative training specifically include: For each sample in the consistent material dataset, spatial topological connections between atoms are established based on a preset cutoff radius to construct a spatial adjacency graph representation of the material configuration; Within the deep learning force field model based on the isovariant graph neural network, a chemical prior feature vector containing atomic number, electronegativity, covalent radius and number of valence electrons is injected into the network through a multilayer perceptron to initialize the node features; A multi-channel Gaussian radial basis function is used to transform the spatial distance between atoms into a high-dimensional vector, and a cosine envelope function is applied to ensure that the distance feature smoothly transitions to zero at the truncation radius boundary, thereby constructing an initial edge feature vector for iterative training in the deep learning force field model based on the isovariant graph neural network.

5. The method of material performance prediction and high-throughput screening based on deep learning force field model according to claim 4, characterized in that, Before constructing a joint loss function and updating model parameters using an automatic differentiation mechanism combined with physical conservation constraints to obtain the trained deep learning force field model, the method further includes performing forward inference in the deep learning force field model based on the isotropic graph neural network: A message passing mechanism is constructed in the hidden layer of the network. The scalar interaction message between the central node and the neighboring nodes is calculated by combining the initial edge feature vector. The three-dimensional spatial coordinate vector of the atom is dynamically updated by multiplying the one-dimensional scalar weight with the relative displacement vector of the atom, thus maintaining the equivariance of translation and rotation of the Euclidean group. The scalar interaction messages passed by neighboring nodes are subjected to a permutation-invariant aggregation operation using a summation function, and the updated high-dimensional node features are output by combining the features of the previous layer of the central node. The energy readout multilayer perceptron module maps the features of the final layer nodes to one-dimensional scalars to obtain the local energy contribution value of the atoms. A scalar summation operation is performed on the local energy contribution values ​​of all atoms in the system to predict and output the total energy of the material system.

6. The method for predicting material properties and high-throughput screening based on a deep learning force field model according to claim 1, characterized in that, The steps of constructing a joint loss function and updating model parameters using an automatic differentiation mechanism combined with physical conservation constraints to obtain the trained deep learning force field model specifically include: The predicted atomic force matrix is ​​obtained by analyzing the negative values ​​of the partial derivatives of the total energy of the material system with respect to the spatial coordinates of each atom in the energy back propagation path through an automatic differentiation mechanism. The predicted unit cell stress tensor is obtained by calculating the derivative of the total energy of the material system with respect to the symmetric strain tensor using the chain rule. The predicted total energy of the material system, the predicted atomic force matrix, and the predicted unit cell stress tensor are integrated and aligned with the normalized system energy, the force vector of all atoms, and the stress tensor in the consistent material dataset to construct a joint mean square error loss function that includes energy error terms and spatial gradient error terms. Adjust the weight parameters of each term in the joint mean squared error loss function according to the different stages of model training, perform global norm truncation after calculating the gradient to limit the gradient from exceeding the preset truncation threshold, call the adaptive moment estimation optimization algorithm to iteratively update the model parameters, and obtain the trained deep learning force field model. The preset truncation threshold is set based on the gradient explosion critical value during the backpropagation process of the deep neural network.

7. The material property prediction and high-throughput screening method based on a deep learning force field model according to claim 1, characterized in that, The steps of inputting the initial crystal structure of the candidate material to be screened into the trained deep learning force field model to obtain the total energy and atomic force vectors, iteratively adjusting the atomic coordinates using an unconstrained optimization algorithm to perform structural relaxation, and outputting a mechanically stable configuration that reaches the mechanical convergence threshold specifically include: The initial crystal structure of the candidate material to be screened is input into the trained deep learning force field model for forward inference, and the total energy and the atomic force vector corresponding to the current configuration are extracted. The atomic force vector and the cell strain gradient transformed based on the inverse matrix of lattice basis vectors are integrated into a high-dimensional global gradient vector. The strain gradient is obtained by multiplying the predicted cell stress tensor by the cell volume of the current configuration after expanding it according to the row priority principle. The search direction is constructed by approximating the inverse Hessian matrix of the objective function using the gradient difference and position difference vector of historical iterations and the quasi-Newton method with limited memory. The evolution step size is determined and the state vector of the candidate material is updated by using an inaccurate line search algorithm that satisfies the strong Wolf criterion. Determine whether the following conditions are met simultaneously: the maximum value of the Euclidean norm of all atomic three-dimensional force vectors is less than or equal to the set force convergence threshold, the maximum value of the absolute value of the elements in the unit cell stress tensor matrix is ​​less than or equal to the set stress convergence threshold, and the absolute difference of the total energy between the two iterations is less than or equal to the set energy difference convergence threshold. If all conditions are met simultaneously, the structure relaxation is determined to be complete and the mechanically stable configuration is output. If all conditions are not met simultaneously, the structure is determined to be incomplete and the iteration continues. The stress convergence threshold is set based on the internal stress tolerance for lattice deformation to stabilize, and the energy difference convergence threshold is set based on the minimum limit of energy fluctuations in adjacent evolution steps.

8. The method for predicting material properties and high-throughput screening based on a deep learning force field model according to claim 1, characterized in that, Before calculating the multidimensional physical properties of the mechanically stable configuration, the method also includes a preliminary stability screening step for the mechanically stable configuration: The average generation energy of each atom in the current mechanically stable configuration is calculated by introducing the chemical potential of the reference element. It is then determined whether the average generation energy is lower than a set thermodynamic threshold. If it is lower, the thermodynamic screening is passed; if it is not lower, the current configuration is rejected. Determine whether the spatial Euclidean distance between any two atoms and the total volume of the unit cell meet the set minimum physical occupancy threshold. If they meet the threshold, the geometric stability screening is passed; otherwise, the current configuration is rejected. The thermodynamic threshold is set based on the thermodynamic barrier of the material against decomposition at room temperature, and the minimum physical occupancy threshold is set based on the empirical radius of the ultimate compressibility of the atom. The Jacobian matrix of the atomic force vector with respect to the atomic coordinates is solved using the backward automatic differentiation mechanism inside the deep learning force field model that has been trained, and used as the dynamic matrix. Eigenvalue decomposition is performed on the dynamic matrix to determine whether all eigenvalues ​​of the square of the corresponding phonon vibration frequency are greater than the set frequency tolerance threshold. If they are all greater than the threshold, the dynamic stability screening is passed. If there are any values ​​that are not greater than the threshold, the current configuration is rejected. The set frequency tolerance threshold is based on the allowable range of numerical error fluctuations in phonon spectrum calculation.

9. The method for predicting material properties and high-throughput screening based on a deep learning force field model according to claim 1, characterized in that, The steps for calculating the multidimensional physical properties of the mechanically stable configuration specifically include: A three-dimensional deformation tensor is applied to the lattice basis vectors of the mechanically stable configuration to obtain a perturbation configuration. The first derivative of stress with respect to deformation rate is predicted using the trained deep learning force field model to construct the elastic stiffness tensor matrix. The macroscopic bulk modulus of polycrystalline materials is derived and calculated using the elastic stiffness tensor matrix and its inverse matrix, employing the Voigt-Royce-Hill approximation model. Phonon vibration frequency characteristics within the Brillouin zone of the reciprocal space of the mechanically stable configuration are extracted. Combined with the Bose-Einstein statistical distribution theory and a preset lower frequency cutoff limit, the isotropic specific heat capacity at the target absolute temperature is calculated by iterative accumulation to obtain the multidimensional physical properties.

10. The method for predicting material properties and high-throughput screening based on a deep learning force field model according to claim 1, characterized in that, The step of constructing a multi-index comprehensive scoring function to calculate the final performance evaluation score by weighting the multi-dimensional physical properties, and then ranking and outputting the candidate materials based on the final performance evaluation score, specifically includes: The multidimensional physical properties are subjected to a direction-aware range standardization operation, which maps positive and negative indicators to values ​​within a set dimensionless scoring range. The multi-index comprehensive scoring function is constructed by combining the pre-assigned weight coefficients with the dimensionless scores of each dimension. The multi-index comprehensive scoring function is then used to perform a weighted linear combination accumulation operation to obtain the final performance evaluation score. Pareto front screening is performed to extract the set of configurations that are not comprehensively superior to other materials in all evaluation metrics. Based on the final performance evaluation score, the tensor lattice basis vectors and coordinates are converted into standard crystallographic information files and written to the local storage medium to complete the sorting and output.

Images