An integrated circuit simulation system and method

Abstract

The application relates to the technical field of integrated circuit reliability simulation and analog technology, in particular to an integrated circuit simulation system and method. The method comprises the following steps: circuit topological structure data of an integrated circuit layout is acquired and processed to obtain an oxide layer area grid map; stress distribution mapping generation is carried out on the integrated circuit according to the oxide layer area grid map to obtain an area stress state matrix; electric field-temperature collaborative analysis is carried out according to the area stress state matrix to obtain a random fluctuation factor matrix; time sequence damage accumulation calculation is carried out according to the random fluctuation factor matrix to obtain a time sequence damage state set; damage state feedback adjustment is carried out on the time sequence damage state set to obtain a time-varying damage evolution graph. The fragile performance index and area of the integrated circuit are identified, and a comprehensive circuit performance degradation characteristic spectrum is generated, so that a scientific basis is provided for reliability design and performance optimization of the integrated circuit.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of integrated circuit reliability simulation technology, and particularly relates to an integrated circuit simulation system and method. BACKGROUND

[0002] The reliability of integrated circuits, especially gate oxide dielectric breakdown (TDDB), is a key factor affecting the long-term stable operation and service life of chips. Existing TDDB analysis methods have made some progress in predicting the average failure time. Most of the existing methods are based on macroscopic models or mean field theory, and fail to simulate the microphysical processes of trap formation, capture, de-trapping and chain evolution in the oxide layer, especially the inherent randomness and non-uniformity of these processes, which leads to the inability to accurately reflect the differences and dynamic changes in damage accumulation in different regions. The existing methods mainly focus on "when" breakdown occurs, i.e., predicting the failure time, but often ignore "how" breakdown occurs, i.e., the inability to simulate how multiple breakdown points form the final conduction path and its specific topology after local damage accumulation reaches a certain level. Due to the inability to simulate the specific morphology and electrical characteristics of the breakdown path, the existing methods cannot real-time feedback the TDDB effect (including the increase in leakage current before breakdown and the conduction path after breakdown) to the circuit simulation model, so as to accurately evaluate the impact of these changes on the actual performance parameters (such as bias point, power consumption, timing, noise margin) of the circuit and even the final function, leading to insufficient comprehensive and accurate reliability evaluation of integrated circuits under complex working conditions.

[0003] In summary, the existing technology has significant deficiencies in comprehensively and dynamically simulating the entire process of TDDB from micro-damage accumulation to macro-circuit failure, especially considering its randomness, breakdown path formation and actual impact on circuit performance, which needs to be solved urgently. SUMMARY

[0004] Therefore, it is necessary to provide an integrated circuit simulation system and method to solve at least one of the above technical problems.

[0005] To achieve the above-mentioned purpose, an integrated circuit simulation method comprises the following steps:

[0006] Step S1: Obtain and dataize the circuit topology structure of the integrated circuit layout to obtain an oxide layer region grid map; generate a region stress state matrix according to the stress distribution mapping of the oxide layer region grid map;

[0007] Step S2: Perform electric field-temperature synergistic analysis according to the regional stress state matrix to obtain a random fluctuation factor matrix; perform time sequence damage accumulation calculation according to the random fluctuation factor matrix to obtain a time sequence damage state set; perform damage state feedback adjustment on the time sequence damage state set to obtain a time-varying damage evolution map;

[0008] Step S3: Perform breakdown point identification and regional adjacency analysis according to the time-varying damage evolution map to obtain a breakdown point space-time distribution table and a regional adjacency relationship graph; perform time sequence connectivity analysis according to the breakdown point space-time distribution table and the regional adjacency relationship graph to obtain a connected path set; construct a breakdown path topological structure according to the connected path set;

[0009] Step S4: Perform performance simulation of the integrated circuit according to the breakdown path topological structure to obtain a time sequence degradation data table; perform performance influence evaluation according to the time sequence degradation data table to obtain a circuit performance degradation characteristic spectrum.

[0010] The present application converts the integrated circuit layout into a fine oxide region grid map, and calculates the electric field and temperature distribution of each region in combination with the initial working condition, which comprehensively and spatially quantifies the initial stress state of all oxide regions in the integrated circuit. This overcomes the limitations of traditional methods that stress evaluation is too macroscopic or only focuses on specific regions, and provides necessary and accurate initial data basis for subsequent accurate simulation of the damage accumulation process of each tiny oxide region, so that the simulation can be closer to the non-uniformity of the actual chip internal stress distribution. The electric field-temperature synergistic acceleration model is combined with a random fluctuation factor with spatial correlation to dynamically simulate the nonlinear process of trap formation and damage accumulation in the oxide layer. By considering the feedback effect of damage on local electric field and temperature at each time step and adjusting the calculation in real time, this step can capture the interaction between stress and damage during the damage evolution process, more truly reflect that the oxide degradation is a dynamic, non-uniform and random process, break through the limitations of relying only on average failure time prediction, and obtain a time-varying damage evolution graph with physical basis. Not only the local breakdown point where the damage reaches the threshold is identified, but also the breakdown path of the conductive channel is simulated and determined by analyzing the connectivity of these points in space and time. By constructing the topological structure of the breakdown path and converting it into an electrical model containing equivalent resistance and capacitance changes, this step provides a specific physical and electrical description of the TDDB failure mode, so that the simulation can go from "when to break down" to "how to break down", overcoming the shortcomings of existing methods that cannot simulate the morphology and electrical characteristics of the conductive path after breakdown. The electrical model of the topological structure of the breakdown path is fed back and updated to the circuit netlist for simulation, which directly quantifies the specific impact and degradation of the TDDB breakdown path on the actual performance parameters (such as timing, bias point, power consumption) of the integrated circuit. By further combining parameter sensitivity and failure threshold for functional failure risk assessment, this step can identify the most vulnerable performance indicators and regions, and generate a comprehensive circuit performance degradation feature spectrum, thereby providing more valuable information than simple lifetime prediction, and providing a scientific basis for reliability design and performance optimization of integrated circuits.

[0011] Therefore, the present application provides an integrated circuit simulation method, which dynamically calculates the accumulation and evolution of oxide damage by introducing a fluctuation factor reflecting the randomness of microscopic physical processes, innovatively simulates and constructs the topological structure of the breakdown path and its electrical equivalent model by analyzing the spatial connectivity, and feeds back the model to the circuit simulation in real time, thereby accurately quantifying and evaluating the actual dynamic impact and degradation trend of the breakdown path on key performance indicators such as circuit bias, timing, power consumption and function, overcoming the shortcomings of existing methods in randomness simulation, path modeling and performance feedback evaluation.

[0012] Preferably, the present application also provides an integrated circuit simulation system for performing the integrated circuit simulation method as described above, the integrated circuit simulation system comprising:

[0013] a stress distribution mapping module for obtaining and digitizing the circuit topology data of the integrated circuit layout to obtain an oxide layer region grid map, and generating a stress distribution mapping of the integrated circuit according to the oxide layer region grid map to obtain a region stress state matrix;

[0014] a damage evolution calculation module for performing a field-temperature collaborative analysis according to the region stress state matrix to obtain a random fluctuation factor matrix, performing a time sequence damage accumulation calculation according to the random fluctuation factor matrix to obtain a time sequence damage state set, and performing a damage state feedback adjustment on the time sequence damage state set to obtain a time-varying damage evolution map;

[0015] a breakdown path analysis module for performing a breakdown point identification and a region adjacency analysis according to the time-varying damage evolution map to obtain a breakdown point space-time distribution table and a region adjacency relationship graph, performing a time sequence connectivity analysis according to the breakdown point space-time distribution table and the region adjacency relationship graph to obtain a connected path set, and constructing a breakdown path topology according to the connected path set;

[0016] a circuit performance evaluation module for performing a performance simulation of the integrated circuit according to the breakdown path topology to obtain a time sequence degradation data table, and performing a performance influence evaluation according to the time sequence degradation data table to obtain a circuit performance degradation characteristic spectrum.

[0017] The integrated circuit simulation system can comprehensively and dynamically simulate the whole process of the dielectric breakdown (TDDB) in the integrated circuit, from the initial stress bearing, the random accumulation evolution of the oxide layer micro-damage, the formation of the physical conduction path, to the influence on the actual performance of the circuit. The system overcomes the shortcomings of the traditional method which only focuses on the average life prediction and ignores the randomness of damage evolution, the breakdown path morphology, and the dynamic influence on the circuit performance. Through accurate spatial stress analysis, a damage accumulation model considering randomness, a breakdown path modeling based on spatial connectivity and its electrical equivalent, and real-time feedback of these information to the circuit simulation, the system can more accurately predict the reliability of the circuit under different stresses, quantitatively evaluate the performance degradation degree, and identify the potential functional failure risk, thereby providing a powerful simulation tool and more instructive data for the reliability design, process optimization, and failure analysis of the integrated circuit. BRIEF DESCRIPTION OF DRAWINGS

[0018] Fig. 1 FIG. 1 is a schematic diagram of the steps of the integrated circuit simulation method according to the present application;

[0019] Fig. 2 is a schematic diagram of an integrated circuit board.

[0020] The purposes, functional features and advantages of the present application will be further illustrated in conjunction with embodiments and with reference to the accompanying drawings. DETAILED DESCRIPTION

[0021] The technical method of the present application will be described clearly and completely below in conjunction with the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all the other embodiments obtained by those skilled in the art without creative work fall within the scope of the present application.

[0022] In addition, the accompanying drawings are only schematic diagrams of the present application, and are not necessarily drawn to scale. The same reference signs in the drawings represent the same or similar parts, and thus repeated descriptions thereof will be omitted. Some of the block diagrams shown in the drawings are functional entities, which do not necessarily have to correspond to physically or logically independent entities. The functional entities can be implemented in the form of software, or in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0023] It should be understood that although the terms "first", "second" and the like can be used herein to describe various elements, these elements should not be limited by these terms. These terms are only used to distinguish one element from another. For example, a first element could be termed a second element, and, similarly, a second element could be termed a first element without departing from the scope of the example embodiments. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0024] To achieve the above-mentioned purposes, please refer to Figs. 1-2 The present application provides an integrated circuit simulation method, comprising the following steps:

[0025] Step S1: obtaining and circuit topological structure data of the integrated circuit layout, and obtaining an oxide layer area grid map; according to the oxide layer area grid map, stress distribution mapping generation is performed on the integrated circuit to obtain an area stress state matrix;

[0026] In the embodiment of the present application, by processing integrated circuit layout data and converting it into a grid representation, all oxide layer regions are identified and their geometric information is obtained to form an oxide layer region grid map. Then, by circuit simulation, the node voltage and current of the circuit in the initial operating state are calculated to obtain initial node electrical parameters. Based on these electrical parameters and the oxide layer grid map, the initial electric field intensity and temperature distribution of each oxide layer region are calculated to generate a regional electric field distribution map and a regional temperature distribution map. At the same time, according to the material type and process information of the oxide layer region, the material database is queried to obtain the regional material property table. Combined with the regional material properties, electric field and temperature information, the initial damage state of each region is evaluated to obtain an initial damage coefficient map. Finally, the electric field intensity, temperature value and initial damage coefficient of each oxide layer region are integrated into a data structure to generate a regional stress state matrix as the initial condition for subsequent damage evolution calculation.

[0027] Step S2: Perform electric field-temperature collaborative analysis according to the regional stress state matrix to obtain a random fluctuation factor matrix; perform time sequence damage accumulation calculation according to the random fluctuation factor matrix to obtain a time sequence damage state set; perform damage state feedback adjustment on the time sequence damage state set to obtain a time-varying damage evolution map;

[0028] In the embodiment of the present application, according to the initial electric field and temperature in the regional stress state matrix, the electric field-temperature collaborative action coefficient of each oxide layer region is calculated, and the regional random reference is set based on this. According to the regional random reference, an initial random number sequence is generated, and the spatial correlation adjustment is performed on the initial random number sequence by analyzing the oxide layer region grid map to establish a spatial adjacency matrix, so that the random fluctuation factor with spatial correlation is obtained. At the same time, according to the regional stress state matrix, the electric field influence factor and the temperature influence factor are calculated, and these factors are combined with the spatially adjusted random sequence to obtain the final random fluctuation factor matrix, which reflects the damage growth trend with randomness under stress driving. Then, according to the random fluctuation factor matrix, the unit time damage increment is calculated, and a nonlinear accumulation function considering the relationship between damage degree and growth rate is constructed. By using the numerical iteration method, the nonlinear accumulation function is integrated and calculated according to the set time step to obtain the cumulative damage value of each region at different time points to form a time sequence damage state set. In this process, according to the damage value in the time sequence damage state set, the regional conductivity change is calculated through the damage-conductivity conversion, the redistribution calculation of electric field and temperature is performed, the feedback influence of damage on local stress state is evaluated, and the electric field and temperature acceleration factor is calculated according to the updated electric field and temperature, so that the acceleration effect of damage accumulation is superimposed to obtain the adjusted stress state sequence. Finally, the time sequence damage state set and the adjusted stress state sequence are mapped to space to generate a time-varying damage evolution map, which dynamically displays the evolution process of damage in space over time.

[0029] Step S3: According to the time-varying damage evolution map, the breakdown point recognition and region adjacency analysis are performed to obtain a breakdown point space-time distribution table and a region adjacency graph; according to the breakdown point space-time distribution table and the region adjacency graph, the time sequence connectivity analysis is performed to obtain a connected path set; and according to the connected path set, the breakdown path topological structure is constructed.

[0030] In the embodiment of the application, the breakdown characteristic parameters of the oxide layer material are obtained, and the damage threshold value triggering local breakdown is set in combination with the cumulative damage value of the region in the time-varying damage evolution map to obtain a breakdown threshold standard table. The time-varying damage evolution map is traversed, the region whose cumulative damage reaches or exceeds the threshold value is identified as a local breakdown point, and the occurrence time and spatial position thereof are recorded to form a breakdown point space-time distribution table. At the same time, according to the oxide layer region grid map, the spatial adjacency relationship between the regions where breakdown occurs is established to obtain a region adjacency graph. Then, according to the breakdown point space-time distribution table, the breakdown events are divided according to the time window to obtain a time window breakdown point table. The breakdown points in the time window are clustered by using the region adjacency graph to obtain a breakdown point clustering table, and the connection relationship between the clusters is determined to obtain a cluster connectivity table. The connected threshold value is calculated according to the cluster connectivity table, the cluster sequence meeting the condition is identified as a potential path to form a potential path table. The potential path table is subjected to time continuity analysis, and the paths stably or continuously existing are screened out to generate a connected path set. Finally, the paths in the connected path set are subjected to dominance evaluation to determine a dominant path to form a dominant path identification map. According to the dominant path identification map and the breakdown point space-time distribution table, the conductivity of the breakdown points constituting the path is calculated, the path is divided into segments, the segment resistance is calculated, and the voltage dependence is considered for correction, and finally the equivalent total resistance of the dominant path is calculated to form a path resistance distribution table. At the same time, the influence of the formation of the dominant path on the original oxide layer capacitance is calculated to obtain capacitance change data. The identification, spatial information, formation time, equivalent resistance and capacitance change of the dominant path and other information are integrated to generate a breakdown path topological structure, which describes the physical and electrical characteristics of the breakdown path.

[0031] Step S4: According to the breakdown path topological structure, the performance simulation of the integrated circuit is performed to obtain a time sequence degradation data table; and according to the time sequence degradation data table, the performance influence evaluation is performed to obtain a circuit performance degradation characteristic spectrum.

[0032] In this embodiment of the invention, the original circuit netlist of the integrated circuit is obtained. Based on the breakdown path topology generated in step S3, corresponding resistors and capacitors are added or modified in the original netlist to form a modified circuit netlist, introducing the electrical effects of the breakdown path into the circuit model. Then, the original circuit netlist is simulated, and various performance indicators of the circuit under non-breakdown conditions are extracted as benchmarks to obtain a benchmark performance indicator table. The modified circuit netlist and the same operating conditions are input into a circuit simulator for simulation. Changes in critical path delay, timing margin, etc., are calculated through timing simulation to obtain a timing degradation data table. The offset of voltage and current at each node is calculated through DC simulation to obtain a bias point drift table. Changes in static and dynamic power consumption of the circuit are calculated based on the bias point drift table to obtain power consumption change data. Next, combining the power consumption change data and the timing degradation data table, key performance parameters that significantly affect functionality are identified, and the sensitivity of these parameters to the underlying breakdown path electrical parameters is calculated to obtain a parameter sensitivity matrix. Failure modes are defined based on the sensitivity matrix, and parameter failure thresholds are determined according to the failure modes and design specifications. A causal relationship graph, or failure causal path graph, is constructed from underlying breakdown to functional failure. On this graph, the probability of occurrence of each node (including failure mode nodes) is calculated based on the degree of performance parameter degradation and thresholds, resulting in a node failure probability table. Finally, based on the node failure probability table, the functional failure risk of the circuit is quantitatively assessed and correlated with circuit regions to generate a functional failure risk graph. By integrating information such as timing degradation, bias drift, power consumption changes, and functional failure risk, a circuit performance degradation characteristic spectrum is generated, comprehensively demonstrating the degradation trends of various key circuit performance characteristics as TDDB evolves and their impact on reliability.

[0033] Preferably, step S1 includes the following steps:

[0034] Step S11: Calculate the initial node electrical parameters of each node of the integrated circuit in its initial operating state;

[0035] Step S12: Generate the regional electric field distribution map and the regional temperature distribution map based on the initial node electrical parameters and the oxide layer region mesh map;

[0036] Step S13: Parameterize the material properties based on the oxide layer region mesh diagram to obtain the region material property table;

[0037] Step S14: Based on the regional material properties, regional electric field distribution map, and regional temperature distribution map, perform an initial damage state assessment to obtain an initial damage coefficient map;

[0038] Step S15: Generate the regional stress state matrix based on the initial damage coefficient map, the regional electric field distribution map, and the regional temperature distribution map.

[0039] In the embodiment of the present application, a complete circuit netlist model of the integrated circuit is needed. The netlist model contains all the components (such as transistors, resistors, capacitors, etc.) and their interconnection relationships in the circuit. Then, the initial working conditions of the circuit are set, such as the power supply voltage V_DD, the ground voltage V_SS, and the initial level of the input signal. The netlist model and the initial working conditions are input into a circuit simulator. The circuit simulator performs DC operating point analysis to calculate the voltage value of each node in the circuit and the current value flowing through each branch in the steady state. For example, for a CMOS inverter, the simulator will calculate the DC voltage of all nodes such as the input terminal, the output terminal, the power rail, and the ground rail, as well as the drain current and the gate current flowing through the transistor. The calculated voltage values and current values of all nodes are summarized to form an initial node electrical parameter table, in which each row corresponds to a circuit node and records its unique identifier, voltage value, and current value.

[0040] According to the oxide region grid map, the spatial position, geometric size (including thickness d), and the circuit nodes connected by the upper and lower electrodes of each oxide region can be determined. The voltage values of these electrode nodes are found in the initial node electrical parameter table. For each oxide region, the voltage difference ΔV between the upper and lower electrode nodes is calculated. Assuming that the electric field is uniform within the oxide layer (simplified model), the average electric field strength E of the region can be calculated as E = ΔV / d. The electric field strengths of all oxide regions in the oxide region grid map are calculated, and these electric field strength values are associated with the corresponding region spatial positions to generate a region electric field distribution map. At the same time, according to the voltage and current values in the initial node electrical parameter table, the power consumption of each component in the circuit is calculated. These power consumptions are the sources of heat. The power consumption distribution and the physical layout of the circuit (spatial information derived from the oxide region grid map) as well as the external heat dissipation conditions are input into a thermal simulator to solve the steady-state heat conduction equation. The thermal simulator calculates the temperature distribution of each region of the circuit in the initial working state. The temperature values of all oxide regions in the oxide region grid map are extracted and associated with the corresponding region spatial positions to generate a region temperature distribution map.

[0041] The oxide region grid map contains not only the spatial information of the oxide regions, but also identifies the material type (e.g., SiO2, HfO2, etc.) of each region and the associated process information (e.g., oxide growth process, anneal condition, etc.). Based on this information, a pre-established material property database is queried. The database stores key material property parameters that affect TDDB, such as dielectric constant ε, bandgap width E_g, intrinsic defect density N0, defect capture cross section σ_c, defect energy level position E_t, etc., for different oxide material types and process conditions. For each oxide region in the oxide region grid map, the corresponding parameter values are extracted from the material property database based on its material type and process information. These parameter values are recorded to form a region material property table, where each row corresponds to an oxide region and records its identifier and all the associated material property parameters.

[0042] The initial damage state refers to the defect level or damage degree that the oxide has before it is put into use or when it first starts to experience stress. This initial damage is caused by process defects during manufacturing or transient stress when first powered on. Based on the material property parameters (e.g., intrinsic defect density N0) of each region in the region material property table, the electric field intensity values E of each region in the region electric field distribution map, and the temperature values T of each region in the region temperature distribution map, an initial damage assessment model is used to calculate the initial damage coefficient D0 for each region. The model can be an empirical formula based on physical mechanisms, for example, a simplified model can be expressed as D0 = g(N0, E, T), where g is a function that reflects the contribution of initial defect density and stress level to initial damage. For each oxide region in the oxide region grid map, its corresponding N0, E, T values are input into the model to calculate its initial damage coefficient D0. The initial damage coefficients calculated for all regions are associated with the corresponding region spatial positions to generate an initial damage coefficient map. The initial damage coefficient is an index that quantifies the initial damage degree of a region, for example, it can be a normalized value between 0 and 1, with 1 indicating the most severe initial damage.

[0043] The regional stress state matrix is a comprehensive data structure for storing the key stress information and initial damage state of each oxide layer region at the initial moment. The data in the initial damage coefficient map, the regional electric field distribution map and the regional temperature distribution map are integrated. Each row or each entry of the matrix corresponds to an oxide layer region in the oxide layer region grid map. For each region, the electric field intensity value E obtained from the regional electric field distribution map, the temperature value T obtained from the regional temperature distribution map, and the initial damage coefficient D0 obtained from the initial damage coefficient map, as well as the unique identifier and spatial position information of the region, are stored in the data entry corresponding to the region. The finally formed regional stress state matrix contains the key initial state information of all oxide layer regions, for example, for N oxide layer regions, the matrix contains N data entries, each entry contains {region ID, spatial coordinates, E, T, D0} fields, and provides complete initial input data for subsequent damage accumulation dynamic calculation.

[0044] Preferably, the electric field-temperature synergistic analysis according to the regional stress state matrix in step S2 comprises:

[0045] The electric field-temperature synergistic coefficient of the oxide layer region is calculated according to the regional stress state matrix to obtain the synergistic coefficient;

[0046] The regional random reference table is set according to the synergistic coefficient;

[0047] The initial random number sequence is generated according to the regional random reference table;

[0048] The spatial adjacency matrix is determined according to the oxide layer region grid map;

[0049] The spatial correlation of the initial random number sequence is adjusted according to the spatial adjacency matrix to obtain the spatial adjusted random sequence;

[0050] The electric field influence factor table is calculated according to the regional stress state matrix;

[0051] The temperature influence factor table is calculated according to the regional stress state matrix;

[0052] The spatial adjusted random sequence, the electric field influence factor table and the temperature influence factor table are subjected to.

[0053] In the embodiment of the present application, the electric field strength E_i and the temperature T_i of each oxide layer region i are extracted from the regional stress state matrix. Then, a model reflecting the synergistic acceleration effect of electric field and temperature is adopted to calculate the synergy coefficient S_synergy,i. For example, an exponential form model can be adopted: S_synergy,i = exp(A / T_i + B·E_i), where A is a parameter related to the thermal activation energy of the material, and B is a parameter related to the electric field acceleration parameter of the material. These two parameters A and B are constants determined in advance according to the material type and process characteristics of the oxide layer. For each region in the regional stress state matrix, the corresponding E_i and T_i values are substituted into the model to calculate the synergy coefficient S_synergy,i of the region. The S_synergy,i values calculated for all regions are summarized to generate a synergy coefficient table, which records the identifier of each region and its corresponding synergy coefficient.

[0054] According to the synergy coefficient S_synergy,i calculated in the previous step, a random reference B_rand,i is set for each region. This random reference quantifies the potential amplitude of random fluctuations in the damage accumulation process of the region under the current stress level. For example, the random reference B_rand,i can be set in linear or nonlinear relationship with the synergy coefficient: B_rand,i = k × S_synergy,i, where k is a scaling factor used to adjust the overall strength of the random fluctuations. This reference value B_rand,i can be interpreted as a parameter of the probability distribution used when generating subsequent random numbers, for example, if a uniform distribution is used, B_rand,i can determine the width or range of the distribution; if a normal distribution is used, B_rand,i can determine the standard deviation. The B_rand,i values set for all regions are summarized to generate a regional random reference table, which records the identifier of each region and its corresponding random reference value.

[0055] According to the B_rand,i values in the regional random reference table, an initial random number R_initial,i is generated for each oxide layer region, which is independent and subject to a specific probability distribution. For example, a random number can be drawn from a uniform distribution with a mean of 1 and a range determined by B_rand,i for each region i. Specifically, R_initial,i can be randomly selected from the interval [1-B_rand,i, 1+B_rand,i]. This initial random number sequence reflects the independent random characteristics of each region before being affected by space. The R_initial,i values generated for all regions are summarized to form an initial random number sequence, which is a sequence containing the same number of random number values as the number of oxide layer regions.

[0056] A spatial adjacency matrix A is constructed by analyzing the spatial information provided by the oxide region grid map. The matrix is an N x N square matrix, where N is the total number of oxide regions. The value of matrix element A_ij is defined as: A_ij = 1 if region i and region j are spatially adjacent (e.g., share a boundary or a vertex); otherwise A_ij = 0. The diagonal elements A_ii are usually set to 0. The determination of spatial adjacency is based on the geometric data in the region grid map, by checking the spatial relationship between the grid cells of different regions. The constructed spatial adjacency matrix A describes the spatial connectivity between all oxide regions.

[0057] Spatial correlation is introduced so that the random fluctuations of adjacent regions are not completely independent. The spatial smoothing or weighted averaging is performed on the initial random number sequence R_initial using the spatial adjacency matrix A obtained in the previous step. For example, for region i, its spatially adjusted random number R_adjusted,i can be calculated as the weighted average of its initial random number and the initial random numbers of its adjacent regions: R_adjusted,i = (1 - w) · R_initial,i + w · (∑{j,A_ij=1}R_initial,j / ∑{j,A_ij=1}1), where w is a weight factor between 0 and 1, controlling the strength of spatial correlation, and ∑_{j,A_ij=1} represents the summation over all regions j that are adjacent to region i. This operation makes the random number values of adjacent regions closer, thereby simulating the spatial correlation of microscopic physical processes. The adjusted R_adjusted,i values of all regions are collected to obtain the spatially adjusted random sequence.

[0058] The electric field intensity E_i of each oxide region i is extracted from the region stress state matrix, and its electric field influence factor F_E,i is calculated. This factor represents the acceleration effect of the electric field on the damage accumulation rate. A formula based on the electric field acceleration model is usually used, for example: F_E,i = exp(γ · E_i), where γ is the electric field acceleration parameter of the material. The calculated F_E,i values of all regions are collected to form the electric field influence factor table.

[0059] The temperature T_i of each oxide region i is extracted from the region stress state matrix, and its temperature influence factor F_T,i is calculated. This factor represents the acceleration effect of the temperature on the damage accumulation rate. A formula based on the thermal activation model is usually used, for example: F_T,i = exp(-E_a / (k_B · T_i)), where E_a is the thermal activation energy of the material, and k_B is the Boltzmann constant. The calculated F_T,i values of all regions are collected to form the temperature influence factor table.

[0060] The spatial adjustment random sequence R_adjusted calculated above, the electric field influence factor table F_E and the temperature influence factor table F_T are integrated to generate a final random fluctuation factor matrix R_fluct. For each region i, its random fluctuation factor R_fluct,i is a function of R_adjusted,i, F_E,i and F_T,i. The random fluctuation factor R_fluct,i will be used in subsequent damage accumulation calculation to dynamically modulate the damage growth rate, reflecting the microscopic random process with spatial correlation driven by stress (electric field and temperature). One integration method is: R_fluct,i = R_adjusted,i × F_E,i × F_T,i. In this way, the final fluctuation factor combines the dual effects of randomness (R_adjusted,i) and stress acceleration (F_E,i, F_T,i). The R_fluct,i values calculated for all regions are summarized to form a random fluctuation factor matrix.

[0061] Preferably, the time-dependent damage accumulation calculation according to the random fluctuation factor matrix in step S2 comprises:

[0062] calculating a unit damage increment table according to the random fluctuation factor matrix;

[0063] constructing a set of nonlinear accumulation functions according to the unit damage increment table;

[0064] performing time point iteration calculation on the set of nonlinear accumulation functions to obtain a set of time-dependent damage states.

[0065] In the embodiment of the present application, the random fluctuation factor R_fluct,i of each oxide layer region i at the current time point t is extracted from the random fluctuation factor matrix. Then, the instantaneous basic damage rate Rate_instant_base,i(t) of each region when not affected by the current damage degree is calculated. The basic damage rate reflects the damage trend under the current stress (electric field and temperature) and the influence of random fluctuation. The calculation method is: Rate_instant_base,i(t) = K_damage × R_fluct,i, where K_damage is the intrinsic damage rate constant of the oxide layer material, a predetermined positive value parameter. The Rate_instant_base,i(t) values calculated for all regions are summarized to form a unit damage increment table, which records the instantaneous basic damage rate of each region at the current time point.

[0066] Based on the physical mechanism of oxide damage accumulation, a nonlinear function f is constructed to describe the relationship between the instantaneous rate of damage accumulation dD_i / dt and the current accumulated damage value D_i(t). The function also depends on the instantaneous basic damage rate Rate_instant_base,i(t) calculated in the previous step. A commonly used nonlinear model is that the damage rate accelerates with the increase of the current damage level, which can be expressed as: dD_i / dt = f(D_i(t), Rate_instant_base,i(t)) = Rate_instant_base,i(t) x exp(a x D_i(t)), where a is a nonlinear coefficient of damage accumulation, usually a positive parameter, reflecting the acceleration effect of damage (e.g. trap formation) on subsequent damage rate. The function f defines the damage accumulation kinetics of each region. The set of nonlinear accumulation functions is the set of such functions f for all regions, with a adjusted for the material properties of different regions.

[0067] Using numerical integration method, iteration calculation is performed according to the time points in the time point sequence table. Starting from the initial time t0, the initial damage value D_i(t0) of each region is known. For each continuous time interval At_k = t_{k+1} - t_k in the time point sequence table, for each oxide layer region i: first, according to the accumulated damage value D_i(t_k) at the current time t_k and the Rate_instant_base,i(t_k) in the unit damage increment table, the actual damage accumulation instantaneous rate Rate_actual,i(t_k) = f(D_i(t_k), Rate_instant_base,i(t_k)) of region i at time t_k is calculated using the nonlinear accumulation function f. Then, the damage increment At_k in this time interval is calculated: At_k ≈ Rate_actual,i(t_k) x At_k. Next, the accumulated damage value of region i at the next time t_{k+1} is updated: D_i(t_{k+1}) = D_i(t_k) + At_k. Repeat this process until all time points in the time point sequence table are calculated. The accumulated damage values D_i(t_k) of all regions at each time point are collected to form a time sequence damage state set, which records the accumulated damage value of each region at each discrete time point.

[0068] Preferably, the damage state feedback adjustment on the time sequence damage state set in step S2 comprises:

[0069] Converting the time sequence damage state set into a region conductivity change table;

[0070] Performing electric field redistribution calculation according to the region conductivity change table to obtain an updated electric field distribution map;

[0071] According to the updated electric field distribution map, local thermal effect evaluation is performed to obtain an updated temperature distribution map;

[0072] According to the updated electric field distribution map, an electric field acceleration factor is calculated;

[0073] According to the updated temperature distribution map, a temperature acceleration factor is calculated;

[0074] According to the electric field acceleration factor and the temperature acceleration factor, acceleration effect superposition is performed on the time sequence damage state set to obtain an adjusted stress state sequence;

[0075] According to the adjusted stress state sequence, a time-varying damage evolution map is generated.

[0076] In the embodiment of the application, the cumulative damage value D_i(t) of each oxide layer region i at the current simulation time t is extracted from the time sequence damage state set. A physical model between damage degree and oxide layer conductivity is established. The model describes how the conductivity increases as the internal defects (damage) of the oxide layer accumulate. For example, a nonlinear model can be used to map the damage value D_i(t) to the conductivity σ_i(t): σ_i(t) = σ_oxide_initial + (σ_oxide_breakdown - σ_oxide_initial) × (D_i(t) / D_breakdown_char) 2 . Wherein, σ_oxide_initial is the initial conductivity of the undamaged oxide layer, σ_oxide_breakdown is the conductivity in the completely breakdown state, and D_breakdown_char is a damage characteristic value representing the breakdown characteristics of the material. For all regions in the time sequence damage state set at the current time, the D_i(t) value is substituted into the model to calculate the conductivity σ_i(t) of each region i at time t. The σ_i(t) values calculated for all regions are summarized to form a region conductivity change table, which records the identifier of each region and its conductivity value at the current time.

[0077] Using the conductivity values in the region conductivity change table obtained in the previous step, the electrical model of the circuit or local region is updated. These updated conductivity values are input as material properties into an electric field solver. The electric field solver solves the electric field distribution based on the spatial geometric information of the oxide layer region grid map and the boundary conditions (for example, the voltage of the upper and lower electrodes of the oxide layer) obtained from the circuit simulation. The solver performs numerical calculations based on partial differential equations inside, such as finite element method or finite difference method, to calculate the electric field distribution that satisfies the current continuity equation and the relationship between electric field and potential potential V and electric field E distribution, where current density J = σE, σ is the value in the region conductivity change table. The calculation result is the electric field intensity distribution of all the oxide layer regions in the circuit at the current time t. The electric field intensity values are associated with the corresponding region spatial positions to generate an updated electric field distribution map.

[0078] According to the updated electric field distribution map and the region conductivity change table, the local power dissipation of each oxide layer region is calculated. For each region i, its unit volume power dissipation P_i(t) can be calculated by Joule's law: P_i(t) = σ_i(t) x |E_i(t)|2. Where σ_i(t) is the conductivity in the region conductivity change table, |E_i(t)| is the electric field intensity amplitude in the updated electric field distribution map. These power dissipations are input as heat sources into a thermal simulator. The thermal simulator solves the heat conduction equation where κ is the thermal conductivity, P is the heat source power density, ρ is the density, and c_p is the specific heat capacity. The simulator takes into account the heat dissipation conditions of the circuit (such as contact with the package, heat sink), and calculates the temperature distribution of the circuit at the current time t. The temperature values of all oxide layer regions are extracted and associated with the corresponding region spatial positions to generate an updated temperature distribution map.

[0079] The electric field intensity E_i(t) of each oxide layer region i at the current time t is extracted from the updated electric field distribution map. An electric field acceleration model, such as an exponential model, is used to calculate the damage acceleration factor F_E_accel,i(t) corresponding to the electric field intensity. The model is usually: F_E_accel,i(t) = exp(γ_E·E_i(t)), where γ_E is the electric field acceleration parameter of the oxide layer material. The F_E_accel,i(t) values calculated for all regions are summarized to form an electric field acceleration factor list.

[0080] The temperature T_i(t) of each oxide layer region i at the current time t is extracted from the updated temperature distribution map. A temperature acceleration model, such as the Arrhenius model, is used to calculate the damage acceleration factor F_T_accel,i(t) corresponding to the temperature. The model is usually: F_T_accel,i(t) = exp(-E_a / (k_B·T_i(t))), where E_a is the thermal activation energy of the oxide layer material, and k_B is the Boltzmann constant. The F_T_accel,i(t) values calculated for all regions are summarized to form a temperature acceleration factor list.

[0081] The electric field acceleration factor F_E_accel,i(t) and the temperature acceleration factor F_T_accel,i(t) calculated above are combined with the random fluctuation characteristics inherent to the region to obtain a comprehensive acceleration factor for damage calculation at the next time step. The new comprehensive random fluctuation factor R_fluct_new,i(t) is calculated as R_fluct_new,i(t) = R_adjusted,i x F_E_accel,i(t) x F_T_accel,i(t) by combining the original spatially adjusted random sequence R_adjusted,i, which represents the random nature inherent to the region. This new factor R_fluct_new,i(t) reflects the damage growth trend of the region at the next time step under the combined action of the actual stress (electric field and temperature) at the current time t and the random nature inherent to the region. The cumulative damage value D_i(t) of all regions at the current time t, the updated electric field E_i(t), the updated temperature T_i(t), and the calculated R_fluct_new,i(t) value are recorded as stress state data at a time point. The set of stress state data at all simulation time points arranged in order constitutes the adjusted stress state sequence. This sequence reflects the dynamic changes in stress state (electric field, temperature, random fluctuation factor) during damage accumulation.

[0082] The data in the adjusted stress state sequence are organized and visualized. For each time point in the sequence, the cumulative damage values D_i(t) of all oxide layer regions and their spatial position information (derived from the oxide layer region grid map) are extracted. These data are mapped onto a two-dimensional or three-dimensional spatial representation to generate a damage distribution map at this time. The damage distribution maps at different time points are arranged in chronological order to form a dynamic atlas that shows how the damage state evolves in space over time. This atlas intuitively shows which regions accumulate damage quickly and which regions accumulate damage slowly, as well as the spatial expansion process of damage. In addition, the time-varying damage evolution atlas can also include updated electric field and temperature distribution maps to comprehensively demonstrate the stress-damage interaction process.

[0083] Preferably, the breakdown point identification and region adjacency analysis according to the time-varying damage evolution atlas in step S3 include:

[0084] Obtaining material breakdown characteristic parameters, setting damage threshold criteria according to the time-varying damage evolution atlas to obtain a breakdown threshold standard table;

[0085] Identifying local breakdown points according to the breakdown threshold standard table to obtain a breakdown point space-time distribution table;

[0086] Establishing region adjacency relationships for the oxide layer region grid map according to the breakdown point space-time distribution table to obtain a region adjacency relationship map.

[0087] In the embodiments of the present application, the key breakdown characteristics of different oxide materials are obtained from a pre-set material characteristics database. These parameters include, for example, critical defect density, critical conductivity, or average time-to-failure under certain stress, which determine the extent of damage at which the material will break down locally. According to the region material type information contained in the time-varying damage evolution map, a specific damage threshold D_threshold,i is set for each oxide region i. The threshold is a quantitative value representing the minimum accumulated damage level required for region i to break down locally under the current material and stress conditions. For example, if the material parameters indicate that a certain oxide layer tends to break down when the damage value reaches 0.7, the D_threshold,i of this region is set to 0.7. This setting process converts the intrinsic breakdown characteristics of the material into a specific damage value threshold that can be used for comparison. The D_threshold,i values of all regions are collected to form a breakdown threshold standard table, which records the identifier of each oxide region and its corresponding damage threshold.

[0088] All time points t and all oxide regions i in the time-varying damage evolution map are traversed. For each accumulated damage value D_i(t) of region i at each time point t, it is compared with the damage threshold D_threshold,i corresponding to this region i in the breakdown threshold standard table. If at a certain time point t, the accumulated damage value D_i(t) of region i first reaches or exceeds its damage threshold D_threshold,i, it is determined that a local breakdown event occurs in region i at this time point. The time and spatial location (i.e., the identifier of region i) of this local breakdown event are recorded. The time and spatial location information of all local breakdown events identified during the simulation are collected to form a breakdown point space-time distribution table. This table is an event list, each record containing an identifier of an oxide region that has broken down locally and a timestamp of the event.

[0089] From the breakdown point space-time distribution table, the identifiers of all the oxide layer regions that have ever experienced local breakdown are extracted. Then, the spatial geometry information of all the oxide layer regions provided by the oxide layer region grid map is utilized. For all the extracted regions that have experienced local breakdown, their spatial mutual position relationship is analyzed. It is determined whether any two of these regions i and j are spatially adjacent, for example, by checking whether they share a side (for a two-dimensional grid) or a face (for a three-dimensional grid) on the grid map. A graph data structure is constructed to represent this spatial adjacency relationship. In the graph structure, each oxide layer region that has experienced local breakdown is a node in the graph, and if two regions are determined to be spatially adjacent, an edge is added between their corresponding nodes. The finally formed graph structure is the region adjacency graph, which describes the spatial connection possibilities between all the potential, locally breakdown-experienced regions, providing a basis for subsequent analysis of the formation of breakdown paths.

[0090] Preferably, the temporal connectivity analysis according to the breakdown point space-time distribution table and the region adjacency graph in step S3 comprises:

[0091] According to the breakdown point space-time distribution table, time window division is performed to obtain a time window breakdown point table;

[0092] According to the time window breakdown point table, in-window breakdown point clustering is performed on the region adjacency graph to obtain a breakdown point clustering table;

[0093] Inter-cluster connection determination is performed on the breakdown point clustering table to obtain a cluster connectivity relationship table;

[0094] A connectivity threshold criterion is calculated according to the cluster connectivity relationship table;

[0095] Initial path identification is performed according to the connectivity threshold criterion to obtain a potential path table;

[0096] Path time continuity analysis is performed according to the potential path table to obtain path continuity data;

[0097] A connected path set is generated according to the path continuity data.

[0098] In the embodiment of the present application, a fixed time window length ΔT_window is first set. Then, all the breakdown events in the breakdown point space-time distribution table are traversed. According to the time stamp of each breakdown event, it is assigned to the corresponding time window. For example, events occurring in [0, ΔT_window) are assigned to window 1, events occurring in [ΔT_window, 2ΔT_window) are assigned to window 2, and so on until the end of the simulation. All the breakdown point events in each time window are organized to form a time window breakdown point table. This table is a view of the original breakdown point space-time distribution table grouped by time windows.

[0099] For each time window in the breakdown point table, extract all the oxide layer regions (i.e. breakdown points) that breakdown within this window. Then, utilize the region adjacency graph. In this graph structure, only consider the regions that breakdown within the current window as nodes. Apply a connected component algorithm (e.g. depth-first search or breadth-first search) on these nodes. If two regions that breakdown within the current window are connected by a series of adjacent regions in the region adjacency graph, they belong to the same cluster. This algorithm identifies all the connected region sets within the current time window that are composed of spatially adjacent breakdown points. Aggregate all the cluster information (including the breakdown point list contained in each cluster) identified in each time window, forming a breakdown point cluster table.

[0100] Traverse the clusters in the breakdown point cluster table, either across different time windows or within the same time window. Determine whether there is a spatial connection between any two clusters C1 and C2. This determination of connection is based on the region adjacency graph: if there exists at least one breakdown point p1 in cluster C1 and at least one breakdown point p2 in cluster C2, and according to the region adjacency graph, p1 and p2 are spatially adjacent, then it is determined that there is a connection between clusters C1 and C2. Record all the pairs of clusters that are determined to have a connection. Aggregate these connection relationships, forming a cluster connectivity table that lists all the pairs of clusters that have a connection and their connection properties (e.g. the number of pairs of adjacent breakdown points in the connection).

[0101] According to the connection information in the cluster connectivity table, set a threshold value that quantifies the "strength" or "importance" of a connection. This threshold value is used to distinguish accidental local breakdown point aggregation from persistent structures that constitute potential conductive paths. For example, the minimum number of pairs of adjacent breakdown points required to connect a pair of clusters can be calculated as the connectivity threshold, or a comprehensive calculation based on the size, density, and other attributes of the clusters. This calculation process converts the original connection data in the cluster connectivity table into a standard value for filtering or evaluating subsequent paths. This standard value is the connectivity threshold standard.

[0102] Utilize the cluster connection relationships in the cluster connectivity table that meet the connectivity threshold standard calculated in the previous step. Construct a new graph structure, where the nodes are breakdown point clusters and the edges are cluster-to-cluster connections that meet the connectivity threshold standard. On this graph, perform a path search algorithm to find sequences of connected clusters from the boundary region on one side of the oxide layer (e.g. connected to the gate) to the boundary region on the other side (e.g. connected to the substrate). These paths are composed of a series of clusters connected head-to-tail. Record all the identified sequences of clusters that meet the conditions, forming a potential path table that lists all the potential breakdown paths across the oxide layer, with each path composed of a series of cluster identifiers.

[0103] The evolution of the identified paths in the potential path table is tracked in different time windows. For a potential path identified in a time window, it is checked whether it still exists or has a corresponding evolved form in the subsequent time window. For example, if a path P_k is identified in a time window k, which is composed of a cluster sequence C_{k,1}, C_{k,2},..., C_{k,m}, it is checked whether there is a path P_{k+1} in the time window k+1 composed of a cluster sequence C_{k+1,1}, C_{k+1,2},..., C_{k+1,n}, and P_k and P_{k+1} have significant overlap or inheritance in spatial position, contained breakdown point set or cluster composition. The temporal continuity of each potential path is quantified, for example, the first appearance time, the last disappearance time, the number of time windows that exist continuously, or the degree of overlap in different time windows are recorded. These quantification results are summarized to form path continuity data.

[0104] A time continuity criterion is set, for example, requiring the number of time windows that the path exists continuously to exceed a certain minimum value N_min_windows, or the path continuity score to be higher than a certain threshold S_min_continuity. Only potential paths that meet the criterion are selected as final connected paths. These selected paths represent stable or continuously existing conductive channels formed during the simulation process. The set of these filtered paths forms a connected path set, which contains all important breakdown paths simulated by the simulation.

[0105] Preferably, the step S3 of constructing the breakdown path topology according to the connected path set comprises:

[0106] The path dominance evaluation is performed on the connected path set to obtain a path dominance ranking table;

[0107] The dominant path identification map is determined according to the path dominance ranking table;

[0108] The path resistance distribution table is calculated according to the dominant path identification map and the breakdown point spatio-temporal distribution table;

[0109] The capacitance change data is calculated according to the dominant path identification map;

[0110] The breakdown path topology is generated according to the dominant path identification map, the path resistance distribution table and the capacitance change data.

[0111] In the embodiments of the present application, for each connected path in the connected path set, the "advantage degree" of the connected path as a leakage current channel is evaluated. The advantage degree can be measured according to the characteristics of the breakdown points or regions constituting the path. For example, the average damage value when all breakdown points on the path breakdown can be calculated, or the average electric field strength or temperature borne by the regions on the path can be calculated. A path constituted by a region with high damage degree and large stress (electric field or temperature) is expected to have stronger conduction ability and higher advantage degree. The specific method is as follows: for a path P_k in the connected path set, it is composed of a series of breakdown points. The information of these breakdown points is obtained from the breakdown point space-time distribution table, and the cumulative damage value D_i_breakdown of these points at the time of breakdown is found from the time-varying damage evolution map. The advantage degree Score_k of the path P_k is calculated as the average of the damage values of these points: Score_k = (∑_ (i∈P_k) D_i_breakdown) / N_k, where N_k is the number of breakdown points included in the path P_k. The advantage degree scores of all paths in the connected path set are calculated and sorted from high to low. The sorting result forms a path advantage degree sorting table.

[0112] According to the path advantage degree sorting table, the path with the highest advantage degree is identified, which is determined as the dominant leakage current path. The dominant path refers to the path of the main conduction channel in the time range of the simulation simulation. For example, the first path in the path advantage degree sorting table is simply selected as the dominant path. A threshold value can also be set, and all paths with an advantage degree higher than a certain value are marked as dominant paths (or important paths). A data structure is created to explicitly identify which paths are determined as dominant paths. The data structure is the dominant path identification map, which contains the identification of all paths in the connected path set and marks whether it is a dominant path.

[0113] For each of the identified dominant paths (and potentially secondary paths), calculate its equivalent resistance. The resistance of a breakdown path depends on the conductivity of the oxide regions that constitute the path after breakdown, as well as the geometry of each region. Obtain the spatial locations and occurrence times of the breakdown sites that constitute the dominant path from the breakdown spatiotemporal distribution map. Look up the cumulative damage values of these sites at the time of breakdown from the time-varying damage evolution map. Convert the damage value of each breakdown site to its corresponding conductivity σ_i_breakdown using a damage-conductivity conversion model (as described in step S2.5). Assume that the path can be modeled as a series or parallel combination of conductive regions. For the series model, calculate the path resistance R_path = ∑_{i∈path} ρ_i_breakdown × L_i / A_i, where ρ_i_breakdown = 1 / σ_i_breakdown is the resistivity of region i at the time of breakdown, L_i is the effective length of the region along the current direction, and A_i is the effective cross-sectional area perpendicular to the current direction. L_i and A_i are determined from the geometry information in the oxide region grid map. Calculate the equivalent resistance values for all dominant paths. Form a path resistance distribution table that records the identification of the dominant paths and their calculated equivalent resistance values.

[0114] Evaluate the impact of dominant path formation on the original capacitance of the oxide layer. The oxide layer, when not broken down, is a dielectric layer with dielectric properties, forming a dielectric capacitance. When a conductive breakdown path is formed, the path "short-circuits" part of the oxide layer in space, changing the dielectric properties of that region and causing the overall equivalent capacitance to change. Determine the spatial location and morphology of the dominant path (derived from the sequence of regions that constitute the path) from the dominant path identification map. Using a capacitance calculation tool, calculate the equivalent capacitance C_initial between the two ends of the oxide layer (e.g., the gate electrode and the substrate) when no breakdown path is formed, based on the oxide region grid map and the material dielectric constant. Then, modify the material properties of the regions in the oxide region grid map that are occupied by the dominant path, for example, treating them as highly conductive conductor regions or setting their dielectric constant to change significantly. Recalculate the modified equivalent capacitance C_broken_down using the capacitance calculation tool. The capacitance change data can be represented as ΔC = C_broken_down - C_initial, or directly record the value of C_broken_down. These calculation results are taken as the capacitance change data.

[0115] All the key information about the dominant breakdown paths calculated and identified above are integrated to form a structured data representation. For each dominant path identified in the dominant path identification map, its equivalent resistance value is looked up from the path resistance distribution table, and its associated capacitance change information is looked up from the capacitance change data. Meanwhile, its spatial location, constituent region sequence, formation time, etc. are obtained from the breakdown point spatio-temporal distribution table and region adjacency graph, etc. underlying data. These information are organized into a hierarchical or object-oriented data structure, for example, a Path object can be defined, which contains: path unique ID, path type (dominant / secondary), list of regions / breakdown points constituting the path (containing spatial coordinates), formation time, calculated equivalent resistance value, and associated capacitance change data. All the Path objects of the dominant paths are collected to form the final breakdown path topology structure. This structure comprehensively describes the physical morphology, spatial location, and equivalent electrical characteristics (resistance and capacitance change) of the breakdown paths, and provides a direct input model for the subsequent circuit performance evaluation.

[0116] Especially important is that the path resistance distribution table is calculated according to the dominant path identification map and the breakdown point spatio-temporal distribution table, specifically:

[0117] The breakdown point conductivity table is calculated according to the breakdown point spatio-temporal distribution table;

[0118] The path segment table is divided from the dominant path identification map according to the breakdown point conductivity table;

[0119] The segment resistance table is calculated according to the path segment table;

[0120] The voltage correction coefficient table is obtained by voltage-dependent correction according to the segment resistance table;

[0121] The path resistance distribution table is generated according to the voltage correction coefficient table;

[0122] In the embodiments of the present application, all local breakdown events recorded in the breakdown point spatiotemporal distribution table are traversed. For each oxide layer region (i.e. breakdown point) i where breakdown occurs, the cumulative damage value D_i(t_breakdown,i) of this region at the time t_breakdown,i when breakdown occurs is looked up from the time-varying damage evolution map. Then, the damage value is converted to the equivalent conductivity σ_i_breakdown of the breakdown point using a pre-established nonlinear physical model between damage degree and oxide layer conductivity. For example, the model can be σ_i_breakdown = σ_initial + (σ_max - σ_initial) x (D_i(t_breakdown,i) / D_char)^β, where σ_initial is the conductivity before breakdown, σ_max is the upper limit of the conductivity when fully on, D_char is a parameter representing the breakdown damage threshold of the material, and β is a nonlinear index. The equivalent conductivities σ_i_breakdown calculated for all breakdown points are associated with the corresponding breakdown point identifiers to form a breakdown point conductivity table.

[0123] The constituent sequence of each dominant path identified in the dominant path identification map is analyzed. A dominant path is composed of a series of oxide layer regions (breakdown points) that are spatially adjacent and have undergone local breakdown. The conductivity values of each breakdown point constituting the dominant path are obtained from the breakdown point conductivity table. According to these conductivity values, the dominant path is divided into several segments. The division principle can be based on the similarity of conductivity, for example, breakdown points along the path that are continuous and have conductivity values within a certain tolerance range are grouped into a segment; or simply each independent breakdown point (or its corresponding oxide layer region) constituting the path is regarded as an independent segment. The segment information of each dominant path (including the breakdown point list contained in each segment

[0124] The equivalent resistance of each segment is calculated. A segment is composed of one or more spatially adjacent breakdown points. For each breakdown point i that constitutes the segment, its conductivity σ i breakdown is obtained from the breakdown point conductivity table, and the spatial geometry information of the region corresponding to the breakdown point is obtained from the oxide region grid map, especially the effective length L i along the path direction and the effective cross-sectional area A i perpendicular to the path direction. It is assumed that the current mainly flows through these breakdown points along the path direction. The equivalent resistance of the breakdown point can be calculated as R i = L i / (σ i breakdown · A i ). If a segment contains multiple breakdown points, and they are in series along the path direction, the total resistance of the segment R segment = ∑{i∈Segment}R i. If they are in parallel in the path cross-sectional direction, the total conductance of the segment G segment = ∑{i∈Segment}σ i breakdown · A i / L i, and the resistance R segment = 1 / G segment. According to the composition of the segment and the current flow assumption, the equivalent resistance value of each segment is calculated. The calculation results are collected to form a segment resistance table, which records the identification of each segment and its calculated equivalent resistance value.

[0125] The nonlinear dependence of the breakdown path's conductive properties on the voltage is considered. For each segment in the segment resistance table, the voltage V segment it is subjected to in the actual circuit working state needs to be evaluated or obtained. Based on the physical model or empirical data of the oxide material, a function R actual = R calculated × f(V segment ) is established to describe how the segment resistance changes with the voltage. The function f(V segment ) is the voltage correction coefficient. For example, for some conduction mechanisms, the resistance decreases as the voltage rises, f(V segment ) = exp(-α V · |V segment |), where α V is the voltage-dependent coefficient. The voltage correction coefficient C voltage segment of each segment under the voltage V segment is calculated. The voltage correction coefficients of all segments are collected to form a voltage correction coefficient table.

[0126] The calculated resistance value in the segmented resistance table is combined with the correction coefficient in the voltage correction coefficient table to obtain the actual equivalent resistance of each segment. For each dominant path in the path segment table, all segments contained therein are sequentially traversed. For each segment s, its calculated resistance R calculated,s is obtained from the segmented resistance table, and its voltage correction coefficient C voltage,s is obtained from the voltage correction coefficient table. The actual equivalent resistance R actual,s of the segment is calculated as R calculated,s x C voltage,s. Then, the actual equivalent resistances of all segments of the dominant path are summed (assuming series connection between segments) to obtain the total equivalent resistance R path_total of the dominant path. The identification of all dominant paths and the total equivalent resistance values calculated therefrom are recorded to form the final path resistance distribution table.

[0127] Preferably, step S4 comprises the following steps:

[0128] Step S41: obtaining an original circuit netlist; performing circuit equivalent updating on the original circuit netlist according to the breakdown path topology to obtain a corrected circuit netlist;

[0129] Step S42: simulating the original circuit to extract a benchmark performance index table;

[0130] Step S43: performing timing degradation simulation according to the benchmark performance index table and the corrected circuit netlist to obtain a timing degradation data table;

[0131] Step S44: calculating a bias point drift table according to the benchmark performance index table and the corrected circuit netlist;

[0132] Step S45: calculating power consumption change data of the circuit according to the bias point drift table;

[0133] Step S46: performing functional failure risk assessment according to the power consumption change data and the timing degradation data table to obtain a functional failure risk graph;

[0134] Step S47: generating a circuit performance degradation characteristic spectrum according to the functional failure risk graph.

[0135] In the embodiment of the present application, the original circuit netlist file of the integrated circuit to be simulated is obtained. The netlist file is a text description of the circuit, containing the types, parameters of all elements in the circuit (e.g. transistors, resistors, capacitors, inductors, voltage sources, current sources, etc.) and the node information of their mutual connections. The netlist is a standard input format of circuit simulation tools (e.g. SPICE format). At the same time, the breakdown path topology data generated in step S3 is obtained. The structure describes the leakage current paths formed during the TDDB evolution process, including the circuit nodes connected by each path (mapped to the nodes connected by the upper and lower electrodes on the oxide layer through the spatial position of the path), the equivalent resistance value of the path and the equivalent change to the original oxide layer capacitance. Then, the original circuit netlist is parsed and modified. Each dominant path in the breakdown path topology structure is traversed. For the path P_k, the two nodes N1 and N2 connected in the original circuit netlist are determined. The equivalent resistance value R_k and the related capacitance change ΔC_k of the path P_k are extracted from the breakdown path topology structure. In the original netlist, a new resistance element with a resistance value of R_k is added between the nodes N1 and N2. At the same time, the equivalent capacitance value between the nodes N1 and N2 (or between the ground / power supply, depending on the connection mode of the original capacitance) is modified according to ΔC_k, for example by adding a parallel capacitance element C_added = ΔC_k. Repeat this process for all dominant paths in the breakdown path topology structure until the electrical effects of all breakdown paths are equivalently added to the original netlist. After modification, save it as a new netlist file, i.e. the modified circuit netlist.

[0136] The simulation parameters of the circuit under typical working conditions are set, such as power supply voltage, input signal waveform, working frequency, environmental temperature, etc. The original circuit netlist and the set simulation parameters are input into a standard circuit simulator (e.g. a simulation tool based on the SPICE core). Perform various types of circuit simulation analysis to obtain the baseline performance indicators of the circuit when no TDDB breakdown occurs. For example, perform DC analysis to calculate the initial stable voltage and branch current (bias point) of each node in the circuit; perform transient analysis, input specific excitation signals (e.g. clock signal, input data), calculate the waveform, delay, rise / fall time, setup / hold time margin of the output signal, and instantaneous and average power consumption in the working period; perform AC analysis or noise analysis to evaluate the frequency response and noise margin. Extract key performance values from the simulation results, such as the propagation delay τ_d of a specific signal path, the clock-to-output delay τ_co, the static power consumption P_static, the dynamic power consumption P_dynamic, the DC voltage V_node of a specific node, the signal swing, etc. These extracted values are summarized and arranged to form a baseline performance indicator table. The table records the reference values of various key performance parameters of the circuit in the ideal non-breakdown state.

[0137] The same circuit operating conditions and input stimuli as the baseline performance simulation are set. The modified circuit netlist and simulation parameters are input to the circuit simulator. A transient simulation is performed, focusing on the timing characteristics of the circuit. In the modified circuit, the propagation delay, clock-to-output delay, setup time margin, and hold time margin of the same key signal paths as in the baseline performance metrics table are measured. For example, the time required for a signal on a critical path to change from input to output is measured. These measured timing values in the modified circuit are compared to the corresponding baseline values in the baseline performance metrics table. The timing degradation, e.g., the increase in delay Δτ_d = τ_d_modified - τ_d_baseline, the decrease in setup time margin Δτ_setup_margin = τ_setup_margin_baseline - τ_setup_margin_modified, is calculated. The degradation values for all key paths and timing margins are recorded to form a timing degradation data table. This table quantifies the impact of the shoot-through path on the circuit switching speed and timing margins.

[0138] The same DC operating conditions as the baseline performance simulation are set. The modified circuit netlist and DC operating conditions are input into the circuit simulator. A DC analysis is performed to calculate the steady-state DC voltages and branch currents at each node in the modified circuit. These new DC operating points reflect the additional leakage current and resistive voltage division effects introduced by the breakdown path. The DC voltage and current values at the key nodes in the modified circuit are compared with the corresponding baseline values in the baseline performance metrics table. The bias point shifts are calculated, e.g., node voltage shift AV node = V node modified - V node baseline, branch current change AI branch = I branch modified - I branch baseline. The voltage and current shift values at all key nodes are recorded to form the bias point shift table. This table quantifies the impact of the breakdown path on the static operating state of the circuit. The same circuit operating conditions as the baseline performance simulation are set. The modified circuit netlist and simulation parameters are input into the circuit simulator. A DC analysis is performed to calculate the static power consumption of the modified circuit at steady state, i.e., the DC power supplied by the power source P static modified. A transient simulation is performed to calculate the dynamic power consumption of the modified circuit over the operating period, P dynamic modified. The increase in static power consumption is mainly due to the leakage current formed by the breakdown path, which flows directly from the power source to ground or substrate. The change in dynamic power consumption is due to the change in node capacitance or signal swing / slew rate caused by the breakdown path. The calculated static and dynamic power consumptions of the modified circuit are compared with the corresponding baseline values in the baseline performance metrics table. The power consumption changes are calculated, e.g., static power consumption increase AP static = P static modified - P static baseline, dynamic power consumption change AP dynamic = P dynamic modified - P dynamic baseline. These power consumption change values are recorded to form the power consumption change data.

[0139] The bias point drift table of step S44 is also considered (although not directly listed in Claim 9, bias drift is an important factor affecting both performance and noise margin). The functional failure criteria of the circuit are set. These criteria are based on the circuit type and application requirements, such as: critical path delay exceeds the target clock period; voltage drift of a certain node exceeds the threshold range recognizable by the subsequent logic gate; the circuit output does not reach the correct logic state within the specified time; total power consumption exceeds the upper limit of the system heat dissipation or power supply capacity; noise margin decreases, leading to increased susceptibility to interference. Based on the timing degradation data table, bias point drift table, and power consumption change data, the performance of the modified circuit is evaluated. For example, check if timing degradation leads to setup / hold time violations; check if bias point drift affects the transistor operating region or logic level; check if the power consumption increase is too large. Based on these evaluation results, identify parts of the circuit that are at risk of functional failure (e.g., specific functional modules, critical paths). Quantify the risk level, for example, a risk scoring mechanism can be set up based on the proximity to the performance degradation and failure threshold. Map the evaluation results to the physical layout of the circuit (using the spatial information provided by the oxide layer region grid map), or present them in a list form. Organize the identified high-risk areas, related performance degradation indicators, and risk level information to form a functional failure risk map.

[0140] Integrate and analyze these scattered performance degradation and risk information. The circuit performance degradation characteristic spectrum is a comprehensive summary that not only indicates the amount of performance change of the circuit under the influence of breakdown, but also reveals the spatial distribution, temporal evolution trend, and potential impact on the final function and reliability of the circuit. The values of timing degradation, bias drift, and power consumption change are associated with the high-risk areas in the functional failure risk map. For example, the spectrum can include: the curve of critical path delay increasing with the evolution time of breakdown; the drift curve of voltage or current at a specific node over time; the curve of total power consumption increasing over time; and the circuit area or working mode with the highest risk of functional failure. The spectrum can be in the form of a chart, table, or structured report. Through the spectrum, designers and reliability engineers can clearly understand how TDDB breakdown specifically affects various performance indicators of the circuit, which performance is the weakest link, and how these degradations ultimately lead to functional failure, thereby providing a basis for improving design and reliability assessment.

[0141] In particular, the functional failure risk assessment based on the power consumption change data and the timing degradation data table is as follows:

[0142] Identify performance critical parameters based on the power consumption change data and the timing degradation data table;

[0143] Calculate the parameter sensitivity matrix of the critical parameters;

[0144] Define failure modes based on the parameter sensitivity matrix;

[0145] determining parameter failure thresholds from the failure mode and parameter sensitivity matrix;

[0146] constructing a failure cause-effect path graph from the parameter failure thresholds;

[0147] calculating a node failure probability table from the failure cause-effect path graph;

[0148] generating a functional failure risk graph from the node failure probability table;

[0149] In the embodiments of the present application, the values in the power consumption variation data and the timing degradation data table are analyzed to identify the performance parameters that have the most significant impact on the circuit functionality and reliability. These parameters are usually the focus of the circuit design specification, or the parameters that show the largest degradation in simulation. For example, the propagation delay of the few critical paths that have the largest delay increase is identified from the timing degradation data table; the increase in static power consumption is identified from the power consumption variation data; if the bias point shift is considered (see step S4.4), the critical node parameters that have the largest voltage or current shift are identified. This identification process can be based on a pre-defined list of critical parameters, or dynamically sorted to select the most drastic varying parameters according to the simulation results. The list of identified performance parameters is taken as the performance critical parameters.

[0150] calculating a parameter sensitivity matrix for the critical parameters. This step evaluates the sensitivity of the performance critical parameters to the underlying circuit parameter variations that lead to the TDDB effect. The underlying circuit parameters here refer to the equivalent resistance, equivalent capacitance, etc. of the breakdown path, whose values are affected by the TDDB damage state. For each performance critical parameter P (e.g. critical path delay τ_d) and each underlying circuit parameter X (e.g. breakdown path resistance R_breakdown) affected by TDDB, the sensitivity of P with respect to X is calculated as The sensitivity can be calculated by parameter sweeping or small signal analysis on the modified circuit netlist. For example, a small perturbation ΔR is applied to the value of the breakdown path resistance R_breakdown in the modified circuit netlist, the circuit simulation is run, and the change in critical path delay Δτ_d is measured, then the sensitivity is approximated as S≈Δτ_d / ΔR. This process is repeated to calculate the sensitivity of all performance critical parameters to all relevant underlying circuit parameters. The results are organized into a matrix, with rows corresponding to performance critical parameters, columns corresponding to underlying circuit parameters, and matrix elements corresponding to the sensitivity values, forming the parameter sensitivity matrix.

[0151] Define failure modes according to parameter sensitivity matrix. This step defines the functional failure modes of the circuit based on the performance parameter response to the underlying changes revealed by the parameter sensitivity matrix. A failure mode is a specific manifestation of the circuit's function deviating from the specification. For example, if a critical path delay is highly sensitive to the decrease of the breakdown path resistance, and the delay increase exceeds the clock period, a "timing violation" failure mode is defined. If the static power consumption is highly sensitive to the increase of the breakdown path conductance, and the power consumption increase exceeds the system heat dissipation capability, a "power consumption overrun" failure mode is defined. By analyzing the parameter sensitivity matrix in combination with the expected function and design specification of the circuit, determine which significant degradation of performance parameters will directly or indirectly lead to specific functional anomalies, and define these anomalies as specific failure modes.

[0152] Determine parameter failure thresholds according to failure modes and parameter sensitivity matrix. This step sets the critical values of the performance key parameters that lead to the occurrence of each defined failure mode, i.e. the parameter failure thresholds. These thresholds are usually derived from the design specification, system requirements, or reliability standards of the circuit. For example, for the "timing violation" failure mode, its parameter failure threshold is the delay value of the relevant critical path reaching or exceeding the working clock period of the circuit. For the "power consumption overrun" failure mode, its parameter failure threshold is the total power consumption of the circuit reaching or exceeding the maximum power consumption allowed by the system. By combining the definition of failure modes and the parameter sensitivity matrix, determine which performance key parameters correspond to each failure mode, and what numerical level these parameters need to reach to trigger the failure mode. These critical numerical values are the parameter failure thresholds.

[0153] Construct failure cause-effect path graph according to parameter failure thresholds. This step constructs a graph model that represents the logical relationship between the underlying causes of performance degradation and the final functional failure. The graph is a directed graph, whose nodes represent events or states (e.g. underlying parameter changes, performance parameter degradation, failure mode occurrence), and edges represent causal relationships. For example, a node represents "breakdown path resistance decrease", which points to the node "critical path delay increase" through an edge, and the node "critical path delay increase" points to the node "timing violation" through another edge. On the nodes corresponding to performance key parameters and failure modes, mark their parameter failure thresholds. This graph clearly depicts the chain reaction of how TDDB effects affect the underlying circuit parameters, then lead to the degradation of performance key parameters, and finally trigger specific failure modes and cause overall functional failure.

[0154] The failure cause and effect graph calculates a node failure probability table. This step quantifies the probability of each event node in the failure cause and effect graph to occur. Using the timing degradation data table and the performance parameter distribution information in the power consumption variation data, in combination with the parameter failure threshold, the probability of a performance critical parameter exceeding its threshold is calculated. For example, if the timing degradation data gives the distribution of the delay of the critical path over time, the probability of the delay value exceeding the clock period at a specific time point can be calculated. The calculation is performed on the failure cause and effect graph using a probability propagation method (e.g., based on Boolean logic or probabilistic inference algorithms). Based on the probability of the underlying events to occur, and the causal relationships between the nodes (defined by the graph structure), the probability of all nodes in the graph to occur is calculated, in particular the nodes representing the failure modes and the final functional failure. The results are aggregated to form the node failure probability table, which lists each event node in the graph and its corresponding probability of occurrence.

[0155] The functional failure risk graph is generated from the node failure probability table. This step translates the probability information in the node failure probability table into an intuitive representation of the functional failure risk. The functional failure risk is typically directly related to the probability of the node representing the final functional failure or the primary failure mode. For example, the probability value of the functional failure node is taken as a quantitative indicator of the overall functional failure risk. This risk value is associated with the affected circuit region or functional module. Using the circuit layout information, the risk information is mapped onto the physical space to generate a heat map or a highlight display, identifying the circuit region with the highest functional failure risk due to TDDB effects. Alternatively, a report or table is generated listing the primary failure modes, their probability of occurrence, potential impact, and the related circuit location. This output is the functional failure risk graph, which provides a comprehensive assessment of the functional reliability status of the circuit under the influence of TDDB, pointing out the most vulnerable parts and failure modes.

[0156] Referring to Fig. 2 Fig. 1 shows a schematic diagram of an integrated circuit board, including a chip 101, a color ring resistor 102, a power resistor 103, and a capacitor 104.

[0157] The chip 101 is an operational amplifier chip that can amplify and operate (e.g., proportional operation, summation operation, etc.) input electrical signals, and is a core device for implementing signal processing and signal amplification functions in analog circuits.

[0158] The color ring resistor 102 and the power resistor 103 function as current limiting and voltage dividing, such as providing appropriate bias voltage for the chip and limiting branch current to make each part of the circuit work in a suitable electrical state. The color ring resistor 102, for example, represents resistance values through different color ring combinations, participates in the construction of resistance networks in signal loops or power supply loops; the power resistor 103 is used in power loops and bears certain power loss to stabilize power-related parameters of the circuit.

[0159] The capacitor 104 functions as filtering (removing noise in the power supply or signal), coupling (transferring alternating current signals and blocking direct current signals), and energy storage (storing a certain amount of energy to stabilize voltage). If it is an electrolytic capacitor, it can filter in the power supply circuit to make the chip power supply more stable, and can couple alternating current signals in the signal circuit to achieve signal transmission.

[0160] The layout includes multiple functional modules and circuit units, which are connected by metal interconnection lines to form a complete circuit structure. The layout clearly divides different areas, including the arrangement of transistor arrays, resistors, capacitors, and the distribution of power supply and ground lines. The overall layout is compact and regular, reflecting the high integration and modularization characteristics in integrated circuit design. It provides a foundation for circuit simulation and reliability analysis.

[0161] Therefore, from any point of view, the embodiments should be regarded as exemplary and non-limiting, the scope of the present application being defined by the appended claims rather than the above description, and it is intended to encompass all variations falling within the meaning and scope of the equivalent elements of the application file. Therefore, the scope of the present application is not limited to the above-mentioned embodiments, but is intended to encompass all variations falling within the meaning and scope of the equivalent elements of the application file.

[0162] The above description is only a specific embodiment of the present application, enabling those skilled in the art to understand or implement the present application. Various modifications to these embodiments will be apparent to those skilled in the art, and the general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present application. Therefore, the present application will not be limited to these embodiments shown herein, but will conform to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A simulation method for integrated circuits, characterized in that, Includes the following steps: Step S1: Obtain and digitize the circuit topology of the integrated circuit layout to obtain the oxide layer region mesh diagram; The stress distribution mapping of the integrated circuit is generated based on the oxide layer region mesh diagram to obtain the region stress state matrix; Step S2: Perform electric field-temperature co-analysis based on the regional stress state matrix to obtain the random fluctuation factor matrix; perform time-series damage accumulation calculation based on the random fluctuation factor matrix to obtain the time-series damage state set; perform damage state feedback adjustment on the time-series damage state set to obtain the time-varying damage evolution map. Step S3: Identify breakdown points and analyze regional adjacency based on the time-varying damage evolution map to obtain a spatiotemporal distribution table of breakdown points and a regional adjacency graph; perform temporal connectivity analysis based on the spatiotemporal distribution table of breakdown points and the regional adjacency graph to obtain a set of connected paths; wherein, the temporal connectivity analysis based on the spatiotemporal distribution table of breakdown points and the regional adjacency graph includes: The time window is divided according to the spatiotemporal distribution table of the breakdown points to obtain the time window breakdown point table; Based on the time window breakdown point table, the region adjacency graph is clustered with breakdown points within the window to obtain a breakdown point clustering table. Perform inter-cluster connectivity determination on the breakdown point cluster table to obtain the cluster connectivity table; Calculate the connectivity threshold standard based on the cluster connectivity table; Initial path identification is performed based on connectivity threshold criteria to obtain a potential path table; Path duration continuity analysis is performed based on the potential path table to obtain path continuity data; Generate a set of connected paths based on path continuity data; construct a breakdown path topology based on the set of connected paths; Step S4: Perform performance simulation on the integrated circuit based on the breakdown path topology to obtain a timing degradation data table; evaluate the performance impact based on the timing degradation data table to obtain the circuit performance degradation characteristic spectrum.

2. The method for simulating integrated circuits according to claim 1, characterized in that, Step S1 includes the following steps: Step S12: Calculate the initial node electrical parameters of each node of the integrated circuit in its initial operating state; Step S13: Generate a regional electric field distribution map and a regional temperature distribution map based on the initial node electrical parameters and the oxide layer region mesh map; Step S14: Parameterize the material properties based on the oxide layer region mesh diagram to obtain the region material property table; Step S15: Based on the regional material properties, regional electric field distribution map, and regional temperature distribution map, perform an initial damage state assessment to obtain an initial damage coefficient map; Step S16: Generate the regional stress state matrix based on the initial damage coefficient map, the regional electric field distribution map, and the regional temperature distribution map.

3. The method for simulating integrated circuits according to claim 1, characterized in that, Step S2, which involves combined electric field-temperature analysis based on the regional stress state matrix, includes: The electric field-temperature synergy coefficient of the oxide layer region is calculated based on the regional stress state matrix to obtain the synergy coefficient. A regional random baseline table is set based on the synergistic effect coefficient; Generate an initial random number sequence based on the regional random baseline table; Determine the spatial adjacency matrix based on the oxide layer region grid diagram; The initial random number sequence is spatially correlated and adjusted based on the spatial adjacency matrix to obtain a spatially adjusted random sequence. Calculate the electric field influence factor table based on the regional stress state matrix; Calculate the temperature influence factor table based on the regional stress state matrix; Tables of spatially adjusted random sequences, electric field influence factors, and temperature influence factors.

4. The method for simulating integrated circuits according to claim 1, characterized in that, Step S2, which involves calculating the time-series damage accumulation based on the stochastic fluctuation factor matrix, includes: Calculate the unit damage increment table based on the random fluctuation factor matrix; Construct a set of nonlinear cumulative functions based on the unit damage increment table; The time-point iterative calculation is performed on the set of nonlinear cumulative functions to obtain the time-series damage state set.

5. The method for simulating integrated circuits according to claim 1, characterized in that, Step S2, which involves adjusting the damage state feedback of the temporal damage state set, includes: A damage-conductivity conversion was performed on the temporal damage state set to obtain a regional conductivity variation table. The electric field redistribution is calculated based on the regional conductivity variation table to obtain an updated electric field distribution map. Local thermal effects are assessed based on the updated electric field distribution map, resulting in an updated temperature distribution map. Calculate the electric field acceleration factor based on the updated electric field distribution map; Calculate the temperature acceleration factor based on the updated temperature distribution map; The adjusted stress state sequence is obtained by superimposing the acceleration effects of the electric field acceleration factor and the temperature acceleration factor on the temporal damage state set. A time-varying damage evolution map is generated based on the adjusted stress state sequence.

6. The method for simulating integrated circuits according to claim 1, characterized in that, Step S3, which involves identifying the breakdown point and resolving the region adjacency based on the time-varying damage evolution map, includes: Obtain material breakdown characteristic parameters, combine them with time-varying damage evolution maps to set damage threshold criteria, and obtain a breakdown threshold standard table; Local breakdown points are identified based on the breakdown threshold standard table to obtain a spatiotemporal distribution table of breakdown points; Based on the spatiotemporal distribution table of breakdown points, the adjacency relationships of the oxide layer region grid are established, resulting in a region adjacency relationship diagram.

7. The method for simulating integrated circuits according to claim 1, characterized in that, Step S3, which involves constructing the breakdown path topology based on the connected path set, includes: The path dominance of the connected path set is evaluated to obtain a path dominance ranking table. Determine the dominant path identification map based on the path dominance ranking table; Calculate the path resistance distribution table based on the dominant path identification diagram and the spatiotemporal distribution table of breakdown points; Calculate capacitance change data based on the dominant path identification map; The breakdown path topology is generated based on the dominant path identification diagram, path resistance distribution table, and capacitance change data.

8. The method for simulating integrated circuits according to claim 1, characterized in that, Step S4 includes the following steps: Step S41: Obtain the original circuit netlist; update the original circuit netlist with circuit equivalence based on the breakdown path topology to obtain the corrected circuit netlist; Step S42: Simulate the original circuit and extract the benchmark performance index table; Step S43: Perform timing degradation simulation based on the baseline performance index table and the modified circuit netlist to obtain the timing degradation data table; Step S44: Calculate the bias point drift table based on the baseline performance index table and the corrected circuit netlist; Step S45: Calculate the power consumption change data of the circuit based on the bias point drift table; Step S46: Perform a functional failure risk assessment based on the power consumption change data and timing degradation data table to obtain a functional failure risk map; Step S47: Generate a circuit performance degradation characteristic spectrum based on the functional failure risk map.

9. A simulation system for integrated circuits, characterized in that, For performing the integrated circuit simulation method as described in claim 1, the integrated circuit simulation system comprises: The stress distribution mapping module is used to acquire and digitize the circuit topology data of the integrated circuit layout to obtain the oxide layer region mesh map; and to generate the region stress state matrix by performing stress distribution mapping on the integrated circuit based on the oxide layer region mesh map. The damage evolution calculation module is used to perform electric field-temperature co-analysis based on the regional stress state matrix to obtain the random fluctuation factor matrix; perform time-series damage accumulation calculation based on the random fluctuation factor matrix to obtain the time-series damage state set; and perform damage state feedback adjustment on the time-series damage state set to obtain the time-varying damage evolution map. The breakdown path analysis module is used to identify breakdown points and resolve regional adjacency based on the time-varying damage evolution map, resulting in a spatiotemporal distribution table of breakdown points and a regional adjacency graph; it performs temporal connectivity analysis based on the spatiotemporal distribution table of breakdown points and the regional adjacency graph to obtain a set of connected paths; and it constructs the breakdown path topology based on the set of connected paths. The circuit performance evaluation module is used to perform performance simulation of the integrated circuit based on the breakdown path topology to obtain a timing degradation data table; and to evaluate the performance impact based on the timing degradation data table to obtain the circuit performance degradation characteristic spectrum.

Images