A method and system for molecular dynamics modeling of a porous medium

Abstract

This invention discloses a molecular dynamics modeling method and system for porous media, belonging to the field of molecular dynamics modeling technology. The method includes the following steps: identifying porous and non-porous voxel regions based on a three-dimensional matrix; dividing multiple non-porous voxel regions into boundary pools and non-boundary layer anchor pools; dividing the boundary layer anchor pools and non-boundary layer anchor pools into multiple grid cells; filling the anchor points within the grid cells with molecules to be filled based on the target filling quantity and target mask value; randomly rotating the current molecule to be filled; if its distance from already placed atoms is less than the corresponding minimum distance threshold, abandoning the current placement; after all molecules are filled, calculating the absolute spatial coordinates of all atoms and outputting a coordinate file. This invention achieves precise, non-overlapping placement of molecules in complex porous structures, enabling the establishment of complex random porous microscopic molecular dynamics models for porous media.

Description

Technical Field

[0001] This invention relates to the field of molecular dynamics modeling technology, and in particular to a molecular dynamics modeling method and system for porous media. Background Technology

[0002] Porous media are materials composed of a solid framework and numerous densely packed micropores separated by this framework. These micropores may be interconnected or partially connected and partially disconnected. Porous media have extremely wide applications, leveraging their complex pore structure to achieve functions such as material separation, adsorption, and fluid transport, finding applications in many key fields including energy, environment, biology, and chemical engineering. By establishing a quantitative correlation between complex microscopic pore structures and macroscopic physical phenomena, we can provide an important foundation for predicting the distribution and flow behavior of fluids within complex pores such as rocks, soil, and biological tissues, thereby solving a series of major engineering, scientific, and environmental problems. Modeling porous media allows us to analyze the impact of pore structure on material transport and energy conversion at the microscopic scale, revealing microscopic mechanisms that cannot be reached by macroscopic experiments, and identifying macroscopic performance control factors to guide optimization in practical applications. By modeling porous media, we can set boundary conditions to simulate extreme scenarios, predict potential risks, and provide a basis for countermeasures, thus avoiding safety hazards in practical applications.

[0003] With the development of pore structure models and the improvement of computational efficiency, establishing the intrinsic connection and quantitative relationship between microscopic mechanisms and macroscopic characteristics, and then using the characteristics of porous media to predict and analyze macroscopic seepage and mass transport behavior, is an important research direction for the future. Molecular dynamics techniques have been widely used to study the behavior and laws of microscopic fluid flow. Molecular dynamics simulations can track and analyze the motion behavior of atoms at the nanometer and micrometer scales, and are an effective means to study the fluid flow behavior inside the micropores of porous media.

[0004] Current microscopic molecular dynamics modeling methods establish regularly shaped pore structures through molecular arrangement and extract microscopic molecular properties as research factors. However, the constructed microscopic models are mainly regular models, and there is a lack of effective molecular filling methods for irregular models. It has always been impossible to break through complex shapes. This has led to the molecular dynamics research of porous media remaining at the level of simplified models for a long time, and it is impossible to establish molecular models of the microscopic pore structure of porous materials that conform to reality. Summary of the Invention

[0005] Based on the shortcomings of the existing technology, the present invention provides a molecular dynamics modeling method and system for porous media, which solves the problem that the existing molecular dynamics research on porous media has long been stuck at the level of simplified models, lacks effective molecular filling methods for complex models, and is difficult to reflect the complex pore molecular microscopic behavior of real materials.

[0006] The present invention adopts the following technical solution:

[0007] In a first aspect, the present invention provides a molecular dynamics modeling method for porous media, comprising the following steps:

[0008] The three-dimensional voxel data of the pore structure of the target porous medium are obtained based on the four-parameter growth method. The three-dimensional voxel data of the pore structure is transformed into a three-dimensional matrix corresponding to the three-dimensional pore structure space. The values ​​in the three-dimensional matrix are assigned mask values ​​of 0 or 1, where 1 represents the pore voxel region and 0 represents the non-pore voxel region.

[0009] Set the size of the simulation box used for molecular dynamics modeling of the target porous medium; within the simulation box, set the type of molecules to be filled, including the number of molecules, the relative coordinates of atoms within the molecules, and the corresponding target mask value; and set the minimum distance between atoms and the anchor point for a single effective spatial position for molecule placement, and obtain the minimum distance threshold between molecules based on the minimum distance between atoms.

[0010] Based on the identification of porous and non-porous voxel regions using a 3D matrix, multiple non-porous voxel regions are delineated to obtain boundary layer anchor pools and non-boundary layer anchor pools. These pools are then divided into multiple grid cells. The target filling quantity for each grid cell is assigned based on the number of molecules filling. Based on the target filling quantity and target mask value, the anchor points within each grid cell are filled using the molecules to be filled. During the filling process, the molecule to be filled is randomly rotated; if its distance from already placed atoms is less than the corresponding minimum distance threshold, the current placement is abandoned. The boundary layer anchor pool is filled first.

[0011] After all molecules have been filled, the absolute spatial coordinates of all atoms are calculated, and the coordinate file is output to complete the molecular dynamics modeling of the target porous medium.

[0012] Preferably, the boundary layer anchor pool consists of multiple non-porous voxel regions directly adjacent to the porous voxel region, and there is no need to divide the anchor pools for the multiple porous voxel regions.

[0013] Preferably, the step of filling the anchor points within the grid cells with molecules to be filled based on the target fill quantity and target mask value specifically includes the following steps:

[0014] Randomly traverse multiple grid cells;

[0015] For the current grid cell, randomly select anchor points within it to fill it;

[0016] After all grid cells are filled, obtain the number of molecules that have been placed; if it is less than the set number of molecules, then place the remaining molecules.

[0017] Preferably, the step of converting the three-dimensional voxel data of the pore structure into a three-dimensional matrix corresponding to the three-dimensional pore structure space includes the following steps:

[0018] The fread function reads uint8 data that matches the size of the simulation box from the three-dimensional voxel data of the pore structure;

[0019] Store uint8 data as a one-dimensional vector, and then reshape the one-dimensional vector into a three-dimensional matrix using the reshape function.

[0020] Preferably, after all molecules have been filled, the absolute spatial coordinates of all atoms are calculated, specifically including the following steps:

[0021] The extreme values ​​of the relative coordinates of the atoms within each molecule are calculated, and the extreme values ​​of each molecule are averaged to obtain the coordinates of the geometric center of each molecule.

[0022] Subtract the corresponding geometric center coordinates from the relative coordinates of all atoms within each molecule to obtain the coordinate difference between each atom within each molecule;

[0023] The absolute spatial coordinates of all atoms within the newly placed molecule are obtained by adding the differences between the geometric center coordinates of the newly placed molecule and the coordinates of each atom.

[0024] Preferably, the acquisition of three-dimensional voxel data of the pore structure of the target porous medium based on the four-parameter growth method specifically includes the following steps:

[0025] An initial model of the target porous medium is established based on the four-parameter growth method. The initial model is then modified to obtain a numerical model of the pore structure.

[0026] The initial model of the target porous medium established based on the four-parameter growth method includes:

[0027] The growth domain is defined, and the porosity, number of solid particles, and generation probability of solid particles are set. Within the growth domain, the solid particles are randomly distributed according to the number of solid particles and generation probability. The grid within the growth domain is traversed to obtain the three-dimensional spatial coordinates of all solid particles.

[0028] Anisotropic growth parameters for solid particles are set, including the growth probability along the principal axis, the growth probability along the face diagonal, and the growth probability along the volume diagonal.

[0029] Using the three-dimensional spatial coordinates of solid particles as the reference node, random numbers of 0-1 are randomly generated uniformly in the adjacent grid cells in each direction around the reference node. If the random number in the current direction is less than the growth probability set in that direction, the particle state of the adjacent grid cells in the current direction is changed to porous particles, and the coordinates of the newly generated particles are stored.

[0030] By iterating through all solid particles, an initial model is obtained.

[0031] Preferably, the step of correcting the pore size of the initial model to obtain a numerical model of the pore structure specifically includes the following steps:

[0032] Spatial connectivity analysis was performed on multiple pores in the initial model to obtain the analysis results; the location and number of openings of each pore on the material surface were detected, and the volume of each pore was obtained;

[0033] Set a minimum pore volume threshold. If the current pore volume is less than the minimum pore volume threshold, remove it. Remove pores with openings on the material edges.

[0034] Based on the analysis results and the number of surface openings, the pores are divided into multi-opening interconnected pores, single-opening pores, and blind pores. Blind pores are removed to obtain a numerical model of the pore structure.

[0035] Preferably, the initial model of the target porous medium established based on the four-parameter growth method further includes:

[0036] Set edge suppression parameters, including edge suppression distance and edge suppression factor, and constrain the position of the pores through the edge suppression parameters.

[0037] Secondly, the present invention provides a molecular dynamics modeling system for porous media, comprising:

[0038] The acquisition module is used to acquire three-dimensional voxel data of the pore structure of the target porous medium based on the four-parameter growth method; the three-dimensional voxel data of the pore structure is converted into a three-dimensional matrix corresponding to the three-dimensional pore structure space, and the values ​​in the three-dimensional matrix are assigned mask values ​​of 0 or 1, where 1 represents the pore voxel region and 0 represents the non-pore voxel region.

[0039] The configuration module is used to set the size of the simulation box for molecular dynamics modeling of the target porous medium; within the simulation box, the type of molecules to be filled is set, including the number of molecules, the relative coordinates of atoms within the molecules, and the corresponding target mask value; and the minimum distance between atoms and the anchor point for a single effective spatial position for molecule placement are set, and the minimum distance threshold between molecules is obtained based on the minimum distance between atoms.

[0040] The filling module is used to identify porous and non-porous voxel regions based on a 3D matrix, perform boundary delineation on multiple non-porous voxel regions to obtain boundary layer anchor pools and non-boundary layer anchor pools; divide the boundary layer anchor pools and non-boundary layer anchor pools into multiple grid cells, allocate the target filling quantity of multiple grid cells based on the number of molecules to be filled, and fill the anchor points in the grid cells with molecules to be filled based on the target filling quantity and target mask value; during the filling process, the current molecule to be filled is randomly rotated, and if its distance from the already placed atoms is less than the corresponding minimum distance threshold, the current placement is abandoned; the boundary layer anchor pool is filled first.

[0041] The output module is used to calculate the absolute spatial coordinates of all atoms after all molecules have been filled, output a coordinate file, and complete the molecular dynamics modeling of the target porous medium.

[0042] Compared with the prior art, the above-mentioned at least one technical solution adopted by the present invention can achieve the following beneficial effects:

[0043] This invention uses a four-parameter growth method to obtain three-dimensional voxel data of pore structures with high similarity to actual porous media. It then fills the porous and non-porous voxel regions of the three-dimensional voxel data with molecules through priority pore boundary filling and mesh generation. Simultaneously, a minimum interatomic distance is set during the filling process to reflect the actual molecular behavior. This invention allows for flexible control of molecular distribution through parameters such as pore boundary layer ratio and mesh dimension, adapting to pore structures of different irregular shapes. It enables precise, non-overlapping placement of molecules in complex pore structures, facilitating the establishment of microscopic molecular dynamics models of complex random pores in porous media. Attached Figure Description

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

[0045] Figure 1 This is a flowchart of a molecular dynamics modeling method for porous media according to the present invention;

[0046] Figure 2 This is a flowchart illustrating the process of establishing the pore structure based on the four-parameter growth method of the present invention.

[0047] Figure 3 This is a flowchart of the molecular filling process of the present invention. Detailed Implementation

[0048] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0049] The four-parameter growth method is a probabilistic simulation approach that randomly distributes solid particles in the computational domain, adjusting porosity and pore shape within the computational region by setting different growth probabilities. Current research using the four-parameter growth method can establish three-dimensional random pore structures, suitable for macroscopic simulation of porous media, but it cannot accurately reflect the actual properties of materials. Many researchers have used molecular dynamics software to establish regularly shaped material pore models, extracting microscopic molecular properties as research factors; however, the constructed microscopic models are predominantly regular, failing to overcome the challenges of complex, irregular shapes. Therefore, microscopic molecular dynamics modeling of complex pores remains a technical challenge.

[0050] This invention aims to propose a molecular dynamics modeling method that conforms to the complex pore structure of actual porous media. A four-parameter growth method is used to establish random pores, and then the pore structure is modified to more closely approximate the real internal structure of porous media. Subsequently, a constructed molecular traversal filling strategy is used for molecular filling.

[0051] This invention provides a four-parameter molecular dynamics modeling method for growing porous media for microscale simulation, and its technical solution is as follows: Figure 1As shown. The overall core method can be divided into two parts: the establishment of a numerical model of random pores and molecular ergonomic filling. A numerical model of random pores can be established using a four-parameter growth method, and various parameters such as porosity, anisotropic growth parameters, solid-phase particle growth probability, and edge suppression growth parameters can be set to achieve preliminary control over pore morphology. After particle growth reaches the set porosity, a series of pore structure corrections are performed, such as edge pore cleaning, removal of small-volume pores, and removal of blind pores, to make the pore structure closer to that of actual porous materials. Furthermore, the bwlabeln function connected region labeling algorithm can accurately identify and analyze the connectivity of the pore structure. Through pore connectivity analysis, pores with different connectivity can be classified, such as multi-opening connected pores, single-opening surface pores, and closed blind pores, and labeled, providing effective methodological guidance for studying different types of complex pores. The core strategy of molecular ergonomic filling is a combination of priority pore boundary filling, mesh generation, and random supplementation methods, with collision detection during particle placement to avoid particle overlap. The molecular traversal filling strategy, through various algorithms, can achieve the filling of three-dimensional random porous systems with molecules. While ensuring a certain degree of uniformity in molecular distribution, it simulates the randomness of molecular existence in macroscopic materials, closely approximating the morphology of actual porous media. Specifically, it includes the following steps:

[0052] S1: Establish a numerical model of random pore size in porous media.

[0053] The process of establishing pore structures based on the four-parameter growth method is as follows: Figure 2 As shown, this includes setting initial growth parameters, distinguishing pore types, and correcting pore structure. The computational flow of the pore generation program based on the four-parameter growth method is as follows:

[0054] S11: First, set the porosity within the predefined growth domain. Then, randomly arrange the solid particles according to the initially set number of solid particles and the probability of solid particle generation. After that, fully traverse the grid system within the growth domain to identify the solid particles and record their three-dimensional spatial coordinates.

[0055] S12: Set the anisotropic growth parameters of the solid particles. The anisotropic growth parameters include the growth probability in the principal axis direction, the growth probability in the face diagonal direction, and the growth probability in the volume diagonal direction.

[0056] Based on the particle growth algorithm model, 26 directions (6 principal axes, 12 face diagonals, and 8 volume diagonals) are used. Anisotropic growth parameters are set, including the growth probability along the principal axes, face diagonals, and volume diagonals. Increasing the growth probability in a certain direction makes pore growth more significant in that direction, while suppressing growth in other directions, achieving a "preferential growth in a specified direction" effect. If the growth probabilities in all directions are equal, a "uniform anisotropic growth" effect is achieved. By setting anisotropic growth parameters, the pore morphology can be controlled.

[0057] S13: Using the stored solid particle coordinates as the reference node, for each of the 26 adjacent mesh elements in different directions, a random number (0-1, uniformly distributed) is generated using the rand function. This number is then compared with the corresponding growth probability set for that direction. If the value is less than the growth probability, the particle state is changed to a porous particle, and the newly generated particle coordinates are stored. Finally, the initial model is obtained by continuously traversing the mesh. For example, the coordinates of a solid particle element are (x... a ,y a ,z a ), where (x a+1 ,y a ,z a ), (x a-1 ,y a ,z a ), (x a ,y a+1 ,z a ), (x a ,y a-1 ,z a ), (x a ,y a ,z a+1 ) and (x a ,y a ,z a-1 Particles along the main axis have growth probabilities P1-P6, (x a+1 ,y a, z a+1 ), (x a-1 ,y a, z a+1 ), (x a ,y a+1 ,z a+1 ), (x a ,y a-1 ,z a+1 ), (x a+1 ,y a ,z a-1 ), (x a-1 ,y a ,z a-1 ), (xa ,y a+1 ,z a-1 ), (x a ,y a-1 ,z a-1 ), (x a+1 ,y a+1 ,z a ), (x a+1 ,y a-1 ,z a ), (x a-1 ,y a+1 ,z a ) and (x a-1 ,y a-1 ,z a For particles located along the face diagonal, the growth probability is P7-P. 18 , (x a+1 ,y a+1 ,z a+1 ), (x a+1 ,y a-1 ,z a+1 ), (x a-1 ,y a+1 ,z a+1 ), (x a-1 ,y a-1 ,z a+1 ), (x a+1 ,y a+1 ,z a-1 ), (x a+1 ,y a-1 ,z a-1 ), (x a-1 ,y a+1 ,z a-1 ) and (x a-1 ,y a-1 ,z a-1 Let ) be a growing particle along the body diagonal, with a growth probability of P. 19 -P 26 .

[0058] Users can set the total porosity to create porous media models with different porosities, and set the number of solid particles to control the dispersion of pores; the more particles, the more dispersed the pore distribution. Users can also set the probability of solid particle formation to control the connectivity of the pores; a higher probability results in better pore connectivity.

[0059] S14: By setting edge suppression growth parameters, including edge suppression distance and edge suppression factor, the pores are kept as far away from the edges as possible, which is more in line with the pore characteristics of actual porous media.

[0060] S15: After the particles reach the initially set porosity, the random pore model generated by the four-parameter growth method can be more consistent with the actual pore structure characteristics of porous media by modifying the pore structure. The modification includes the identification of blind pores, the removal of extremely small volume pores, and the removal of edge-opening pores. Spatial connectivity analysis of multiple pores in the initial model is performed based on the bwlabeln function. Subsequently, the location and number of openings of each pore on the material surface are detected by traversing six surfaces, and the volume of each pore is calculated. Setting a minimum pore volume threshold can remove pores with extremely small volume sizes, improving the realism of the pore structure. Based on connectivity and the number of surface openings, pores can be divided into three categories: multi-opening connected pores, single-opening pores, and blind pores, and are labeled accordingly, providing methodological guidance for studying pores with different connectivity. Pores in porous media materials can be roughly divided into multi-opening through pores with good connectivity and single-opening pores with poor connectivity based on their connectivity. Blind pores (pores without openings) enclosed within the material are meaningless for the function of porous materials, so they are removed after being identified by connectivity analysis to avoid interference from unrealistic pores in the study of porous material performance. Pores with openings at the material edges are detected, but since these types of pores are extremely rare in porous media, they can be removed.

[0061] S16: After a series of pore structure corrections, the porosity is recalculated, including the porosity of various pore types after connectivity analysis. Since the corrections are only minor adjustments to the pore structure and do not significantly change the overall morphology, the error between the corrected porosity and the initially set total porosity is extremely small, within a reasonable and feasible range. Finally, a three-dimensional numerical model of random pores is obtained.

[0062] S17: Save the random pore 3D numerical model as a raw binary file, i.e., pore structure 3D voxel data.

[0063] In summary, the numerical model of random pores established based on the four-parameter growth method has a high degree of similarity with actual porous media materials, and the output three-dimensional voxel data provides a reliable foundation for subsequent molecular filling.

[0064] S2: Molecular traversal filling.

[0065] The molecular filling strategy is the focus of this invention, and the molecular filling process is as follows: Figure 3 As shown, the core filling strategy combines priority pore boundary filling, mesh generation, and random supplementation. It achieves precise, non-overlapping molecular placement within the porous structure through four-parameter grown 3D voxel data reading, parameter setting, step-by-step molecular placement, and collision detection. Specifically, it includes the following steps:

[0066] S21: Read the three-dimensional voxel data of the pore structure.

[0067] First, the raw binary voxel data file of the pore structure output by the four-parameter growth method is read, laying the foundation for the subsequent reshaping of the three-dimensional numerical model of the pore structure, thereby providing the target area for molecular placement.

[0068] S22: Set the simulation box size `boxSize`. `boxSize` defines the number of voxels in the simulation space along the X, Y, and Z directions. `boxSize` must match the size of the growth domain simulated by four-parameter growth. To ensure the efficiency and accuracy of molecular filling, the pore structure volume output by four-parameter growth needs to be proportionally enlarged. The pore morphology remains unchanged after proportional enlargement, making the overall positional distribution of the filled molecules more consistent with the actual target expectation. Proportional enlargement of the pore volume is achieved by setting the actual voxel size `voxelSize`. `voxelSize` represents the side length of each smallest unit in the three-dimensional voxel data in space.

[0069] S23: Then set the radius of each atom and the total number of each type of molecule. Considering the interatomic interactions involved in molecular dynamics, in order to avoid excessive energy caused by atoms being too close together, it is not enough to simply ensure that the molecules do not overlap when placing them. It is necessary to additionally set the minimum distance between atoms to conform to the actual situation of molecules.

[0070] The principle behind minimum interatomic distance detection is achieved through the `noOverlapDistanceBased` function. This function uses the `atomLocationGrid` parameter, a spatial grid index, to detect the distance between a new atom and existing atoms. It then matches the distance to a corresponding minimum distance threshold based on the molecule type. If the distance is less than the threshold, the atom is considered to be overlapping, and the current placement is abandoned. The function does not iterate through all existing atoms; it only detects atoms within a certain grid range adjacent to the new atom, significantly improving detection efficiency.

[0071] S24: Define the molecule type using the `chainTypes` parameter, storing the molecule type name, number of molecules, relative atomic coordinates, atomic labels, and target voxel mask value (0 or 1). The difference between filling random pores with molecules and single atoms lies in the fact that the former has a specific internal atomic distribution. The internal structure of each molecule is different, such as the number of atoms and their relative positions. Therefore, the number of atoms and detailed relative position coordinates inside the molecule need to be defined through the `chainTypes` parameter. These atoms are treated as a whole, representing a molecule, for filling. This is one of the main core innovations of this invention.

[0072] S25: Initialize data structures to store global variables containing information about all placed atoms. Initialize an empty matrix `allAtoms` to store the 3D physical coordinates of all atoms. Initialize an empty cell array `allLabels` to store the labels of all atoms. Initialize an empty cell array `allTypes` to store the molecular type to which all atoms belong. Initialize a 3D cell array `atomLocationGrid`, where each cell stores the index of the atom falling within that grid for overlap detection.

[0073] S26: Subsequently, the pore voxel data size is verified and reshaped. The 8-bit integer uint8 data matching the simulated box size boxSize is read from the pore voxel data file using the fread function and stored as a one-dimensional vector volume. It is verified that the amount of data read is equal to the product of the three dimensions of boxSize. The one-dimensional vector volume is then reshaped into a three-dimensional matrix (with the same dimensions as boxSize) using the reshape function, so that it corresponds to the three-dimensional pore structure space.

[0074] S27: The principle of molecule placement is to ensure uniform distribution and accurate regional distribution. To ensure the accuracy of each molecule placement region, a corresponding mask value is assigned to each molecule. The four-parameter growth method assigns a mask value of 1 to the generated porous voxel region and a mask value of 0 to the matrix material voxel region outside the porous region. In this molecule filling stage, each type of molecule needs to be assigned a mask value corresponding to the voxel region constructed by the four-parameter growth method to ensure the accuracy of the molecule filling region. The mask coordinate index construction step creates a physical coordinate index for different voxel regions, realizing an accurate mapping from "voxel region" to "spatial coordinates," providing a basis for subsequent molecule anchor point selection. Due to the complexity of the atomic structure inside the molecule, the shape of the molecular strip is extremely irregular compared to a single atomic sphere. This leads to a significant difference between the pore morphology formed after the actual placement of molecules in region 0 and the pores established by the four-parameter growth method, thus increasing the porosity and causing errors. Therefore, when placing molecules in non-porous region 0, prioritizing the placement of pore boundaries can ensure that the original pore orbital morphology remains unchanged as much as possible, improving the rationality of the overall pore structure.

[0075] S28: Based on the different voxel mask values, anchor pools are divided into 0-anchor pools and 1-anchor pools. A 0-anchor pool is defined as a "voxel with a mask value of 0," and a 1-anchor pool is defined as a "voxel with a mask value of 1." The 0-anchor pool type classification step is specifically designed for molecules filling non-porous regions 0. By identifying "boundary voxels" and "non-boundary voxels," two types of 0-anchor pools are classified to achieve priority filling of molecules in pore boundary regions. The boundary layer 0-anchor pool is defined as "0-region voxels directly adjacent to the pore region." The coordinates of the boundary layer 0-anchor pool are obtained by performing an "expansion search" of voxels of a certain size in the pore region using the `imdilate` function. The non-boundary layer 0-anchor pool is defined as "the remaining voxels in the 0 region excluding the boundary layer," serving as the next filling region for this type of molecule. For molecules in region 1, during the filling process, the 1-anchor pool does not need further classification; all molecules in region 1 are directly placed in the 1-anchor pool. An anchor point is a single, valid spatial location selected from all voxel data for molecule placement. It is a precise point within the mask value region corresponding to the molecule category. The anchor point is the geometric center of the molecule, and all atoms within a molecule rotate and are placed around this anchor point. The extreme values ​​of the X / Y / Z relative position coordinates of all atoms within the initially set molecule are calculated, and the average of these extreme values ​​yields the geometric center coordinates of the molecule. The geometric center coordinates are then subtracted from the relative coordinates of all atoms within the molecule to obtain the coordinate differences for each atom. When placing subsequent molecules of the same type, the anchor point coordinates of the newly placed molecule are added to the calculated coordinate differences for each atom to obtain the absolute spatial coordinates of all atoms within the newly placed molecule.

[0076] Within a spatial region, a collection of multiple anchor points constitutes an anchor point pool, which is generated in the `buildIndexMap` function. During molecule placement, the system selects anchor points from the pool as candidate placement locations. Not all anchor points become placement locations for the molecule; only anchor points that pass the molecular placement rules (minimum interatomic distance control, random rotation) are actually used to place the molecule.

[0077] S29: The core strategy for molecule placement is a combination of grid partitioning and random replenishment. The core principle of grid partitioning is to divide the target space (anchor pool) for molecule filling into multiple small grid units. Within each unit, the number of molecules to be filled is evenly distributed, and anchor points are randomly selected. This achieves a uniform distribution of molecules throughout the overall space. Simultaneously, "local search within units" reduces the complexity of anchor point selection, balancing filling efficiency and distribution uniformity. The rule for filling non-porous region 0-type molecules during the grid partitioning molecular placement stage is to prioritize placement from the pore boundary layer anchor pool, followed by placement from the non-boundary layer anchor pool, ensuring a sufficient number of 0-type molecules are placed at the boundary layer between the pores and the matrix. Furthermore, the number of 0-type molecules to be placed in the pore boundary layer can be set according to the porosity. The core objective of grid partitioning is to solve the problem of uniformity and efficiency. During molecule filling, directly randomly selecting anchor points from the entire anchor pool may result in localized molecule density and gaps in some areas. Traversing each anchor point individually would lead to inefficiency due to the large number of anchor points. The grid-based approach addresses two key issues simultaneously through spatial partitioning, task allocation, and local filling: it assigns a fixed proportion of molecular fill to each grid cell to ensure uniformity, and it selects anchor points only within the current grid cell, eliminating the need to traverse the entire anchor point pool, thus reducing computation and improving placement efficiency. The `gridDiv` parameter represents the number of grid divisions in the X, Y, and Z directions, dividing the target space into multiple equal small cells and calculating the physical dimensions of each small grid cell. The number of grid cells is controlled by the `gridDiv` parameter; a higher number of cells results in a more uniform molecular distribution, but slightly increases computation. The `gridDiv` parameter typically needs to be adjusted based on the size of the anchor point pool.

[0078] S30: To further ensure uniform filling, a random grid cell traversal method is adopted. This method does not traverse grid cells in a fixed order, but generates a random cell index order and processes each cell in a random order. If traversed in a fixed order, the cells processed earlier may be filled with more molecules due to sufficient anchor points, while the cells processed later may be filled with empty cells due to insufficient anchor points. Therefore, random traversal can avoid this order bias and give each cell an equal chance to be filled. For each randomly selected grid cell, the anchor points falling within the cell are first selected, and then the target filling quantity for the cell is allocated according to "number of molecules filled / total number of cells", completing the local filling within the cell. The custom function fillChainsInGrid is called to place molecules in the grid. Due to the irregularity of molecular morphology, the molecules to be placed need to be rotated to create an "interlacing" effect in the spatial arrangement of molecules, improving the success rate of placement while closely reflecting the randomness of the molecular morphology. Molecular rotation is achieved by generating a random rotation matrix using the randRot function, which generates a three-dimensional random rotation matrix to simulate the random orientation of molecules. The principle of the randRot function is to randomly generate three Euler angles ( θ, φ, ψ To construct rotation matrices (Rz1, Ry, Rz2) around the Z-axis, Y-axis, and Z-axis respectively, the final rotation matrix is ​​obtained by strictly multiplying them in the order Rz2×Ry×Rz1. Regardless of whether it's the grid method or the supplementary method, after the molecule undergoes rotation attempts, it must sequentially pass boundary checks, region mask value checks, and atomic spacing checks before it can be successfully placed.

[0079] S31: After the mesh filling stage, calculate the number of placed molecules (placedCount). If it is less than the set number of molecules (chainCount), then place the remaining molecules. This is because the initial mesh generation method cannot guarantee that all molecules can be successfully placed. Firstly, the irregular shape of the pore structure causes some meshes to cover voxels with different mask values, resulting in insufficient anchor pool space for that molecule mask type within that mesh. Secondly, if an anchor point in a mesh still overlaps with existing atoms after rotation, that anchor point will be abandoned, causing the mesh to fail to reach the "target placement number." In the supplementation stage, the rule for filling non-pore region 0 molecules remains the same: prioritize supplementation from the pore boundary layer anchor pool, then from the non-boundary layer anchor pool. Molecules needing to fill pore region 1 are directly filled using type 1 meshes. For each molecule to be supplemented, the placement rule is the same as in the mesh generation stage: randomly traverse the anchor pool and anchor points, attempt several random rotations for each anchor point, and then successfully place it after checking the atomic spacing and region mask value.

[0080] S32: After all molecules are successfully placed, the relative position coordinates of the atoms are converted into absolute spatial coordinates. The absolute spatial coordinates of all atoms are counted, and the coordinate XYZ file is output. Valid data is extracted from the global variables allAtoms (absolute coordinate matrix), allLabels (atom label array), and allTypes (molecule type). Each atom and molecule is traversed, and a text line is generated for each molecule in the format of "atom label + X coordinate + Y coordinate + Z coordinate". The corresponding molecule type can be written according to the pore type, providing conditions for subsequent research on different types of pore molecular dynamics.

[0081] The core difference between the particle traversal filling strategy proposed in this invention and traditional methods lies in the rational application of multiple filling algorithms and rules to achieve high-efficiency and high-accuracy filling of random porous structures, making the molecular dynamics model closer to the actual material properties. This method employs a composite strategy of "layer-first + two-stage filling," which allows for flexible control of molecular distribution through parameters such as the proportion of the pore boundary layer and mesh dimension, adapting to different pore structures. Traditional methods lack the controllability of molecular distribution. These features make this method superior to traditional methods in terms of the rationality and practicality of molecular distribution, filling efficiency, and stability, meeting the demand for precise filling of large-scale molecules in complex porous structures.

[0082] Based on the same concept, the present invention also provides a molecular dynamics modeling system for porous media, including an acquisition module, a setting module, a filling module and an output module.

[0083] The acquisition module is used to acquire three-dimensional voxel data of the pore structure of the target porous medium based on the four-parameter growth method; the three-dimensional voxel data of the pore structure is converted into a three-dimensional matrix corresponding to the three-dimensional pore structure space, and the values ​​in the three-dimensional matrix are assigned mask values ​​of 0 or 1, where 1 represents the pore voxel region and 0 represents the non-pore voxel region.

[0084] The setting module is used to set the size of the simulation box for molecular dynamics modeling of the target porous medium; within the simulation box, the type of molecules to be filled is set, including the number of molecules, the relative coordinates of atoms within the molecules, and the corresponding target mask value; and the minimum distance between atoms and the anchor point for a single effective spatial position for molecule placement are set, and the minimum distance threshold between molecules is obtained based on the minimum distance between atoms.

[0085] The filling module is used to identify porous and non-porous voxel regions based on a 3D matrix. It performs boundary delineation on multiple non-porous voxel regions to obtain boundary layer anchor pools and non-boundary layer anchor pools. The boundary layer anchor pools and non-boundary layer anchor pools are then divided into multiple grid cells. The target filling quantity for each grid cell is allocated based on the number of molecules to be filled. Based on the target filling quantity and target mask value, the anchor points within the grid cells are filled using the molecules to be filled. During the filling process, the molecule to be filled is randomly rotated. If its distance from already placed atoms is less than the corresponding minimum distance threshold, the current placement is abandoned. The boundary layer anchor pool is filled first.

[0086] The output module is used to calculate the absolute spatial coordinates of all atoms after all molecules have been filled, output a coordinate file, and complete the molecular dynamics modeling of the target porous medium.

[0087] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention.

[0088] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A molecular dynamics modeling method for porous media, characterized in that, Includes the following steps: Three-dimensional voxel data of the pore structure of the target porous medium were obtained based on the four-parameter growth method. The three-dimensional voxel data of the pore structure is transformed into a three-dimensional matrix corresponding to the three-dimensional pore structure space. The values ​​in the three-dimensional matrix are assigned mask values ​​of 0 or 1, where 1 represents the pore voxel region and 0 represents the non-pore voxel region. Set the size of the simulation box used for molecular dynamics modeling of the target porous medium; within the simulation box, set the type of molecules to be filled, including the number of molecules, the relative coordinates of atoms within the molecules, and the corresponding target mask value; It also sets the minimum interatomic distance and the anchor point for a single effective spatial position for molecule placement, and obtains the minimum distance threshold between molecules based on the minimum interatomic distance; Based on the identification of porous voxel regions and non-porous voxel regions using a three-dimensional matrix, boundary delineation is performed on multiple non-porous voxel regions to obtain boundary layer anchor pools and non-boundary layer anchor pools. The boundary layer anchor pool and non-boundary layer anchor pool are divided into multiple grid cells. The target filling quantity of multiple grid cells is allocated based on the number of molecules filled. Based on the target filling quantity and the target mask value, the anchors in the grid cells are filled by the molecules to be filled. During the filling process, the molecules to be filled are randomly rotated. If the distance between the molecule to be filled and the already placed atom is less than the corresponding minimum distance threshold, the current placement is abandoned. The boundary layer anchor pool is filled first. After all molecules have been filled, the absolute spatial coordinates of all atoms are calculated, and the coordinate file is output to complete the molecular dynamics modeling of the target porous medium.

2. The molecular dynamics modeling method for porous media as described in claim 1, characterized in that, The boundary layer anchor pool consists of multiple non-porous voxel regions directly adjacent to the porous voxel region. There is no need to divide the anchor pool for multiple porous voxel regions.

3. The molecular dynamics modeling method for porous media as described in claim 1, characterized in that, The process of filling anchor points within grid cells using molecules to be filled, based on the target fill quantity and target mask value, specifically includes the following steps: Randomly traverse multiple grid cells; For the current grid cell, randomly select anchor points within it to fill it; After all grid cells are filled, obtain the number of molecules that have been placed; if it is less than the set number of molecules, then place the remaining molecules.

4. The molecular dynamics modeling method for porous media as described in claim 1, characterized in that, The process of converting the three-dimensional voxel data of the pore structure into a three-dimensional matrix corresponding to the three-dimensional pore structure space includes the following steps: The fread function reads uint8 data that matches the size of the simulation box from the three-dimensional voxel data of the pore structure; Store uint8 data as a one-dimensional vector, and then reshape the one-dimensional vector into a three-dimensional matrix using the reshape function.

5. The molecular dynamics modeling method for porous media as described in claim 1, characterized in that, After all molecules have been filled, the absolute spatial coordinates of all atoms are calculated, which includes the following steps: The extreme values ​​of the relative coordinates of the atoms within each molecule are calculated, and the extreme values ​​of each molecule are averaged to obtain the coordinates of the geometric center of each molecule. Subtract the corresponding geometric center coordinates from the relative coordinates of all atoms within each molecule to obtain the coordinate difference between each atom within each molecule; The absolute spatial coordinates of all atoms within the newly placed molecule are obtained by adding the differences between the geometric center coordinates of the newly placed molecule and the coordinates of each atom.

6. The molecular dynamics modeling method for porous media as described in claim 1, characterized in that, The method for obtaining three-dimensional voxel data of the pore structure of the target porous medium based on the four-parameter growth method specifically includes the following steps: An initial model of the target porous medium is established based on the four-parameter growth method. The initial model is then modified to obtain a numerical model of the pore structure. The initial model of the target porous medium established based on the four-parameter growth method includes: The growth domain is defined, and the porosity, number of solid particles, and probability of solid particle generation are set. Within the growth domain, the solid particles are randomly distributed according to the number of solid particles and the probability of solid particle generation. The grid within the growth domain is traversed to obtain the three-dimensional spatial coordinates of all solid particles. Anisotropic growth parameters for solid particles are set, including the growth probability along the principal axis, the growth probability along the face diagonal, and the growth probability along the volume diagonal. Using the three-dimensional spatial coordinates of solid particles as the reference node, random numbers of 0-1 are randomly generated uniformly in the adjacent grid cells in each direction around the reference node. If the random number in the current direction is less than the growth probability set in that direction, the particle state of the adjacent grid cells in the current direction is changed to porous particles, and the coordinates of the newly generated particles are stored. By iterating through all solid particles, an initial model is obtained.

7. The molecular dynamics modeling method for porous media as described in claim 6, characterized in that, The process of correcting the pore structure of the initial model to obtain a numerical model of the pore structure includes the following steps: Spatial connectivity analysis was performed on multiple pores in the initial model to obtain the analysis results; the location and number of openings of each pore on the material surface were detected, and the volume of each pore was obtained; Set a minimum pore volume threshold. If the current pore volume is less than the minimum pore volume threshold, remove it. Remove pores with openings on the material edges. Based on the analysis results and the number of surface openings, the pores are divided into multi-opening interconnected pores, single-opening pores, and blind pores. Blind pores are then removed, resulting in a numerical model of the pore structure.

8. The molecular dynamics modeling method for porous media as described in claim 6, characterized in that, The initial model for establishing the target porous medium based on the four-parameter growth method also includes: Set edge suppression parameters, including edge suppression distance and edge suppression factor, and constrain the position of the pores through the edge suppression parameters.

9. A molecular dynamics modeling system for porous media, characterized in that, include: The acquisition module is used to acquire three-dimensional voxel data of the pore structure of the target porous medium based on the four-parameter growth method. The three-dimensional voxel data of the pore structure is transformed into a three-dimensional matrix corresponding to the three-dimensional pore structure space. The values ​​in the three-dimensional matrix are assigned mask values ​​of 0 or 1, where 1 represents the pore voxel region and 0 represents the non-pore voxel region. The configuration module is used to set the size of the simulation box for molecular dynamics modeling of the target porous medium; within the simulation box, the type of molecules to be filled is set, including the number of molecules, the relative coordinates of atoms within the molecules, and the corresponding target mask value; It also sets the minimum interatomic distance and the anchor point for a single effective spatial position for molecule placement, and obtains the minimum distance threshold between molecules based on the minimum interatomic distance; The filling module is used to identify porous voxel regions and non-porous voxel regions based on a three-dimensional matrix, perform boundary division on multiple non-porous voxel regions, and obtain boundary layer anchor pools and non-boundary layer anchor pools. The boundary layer anchor pool and non-boundary layer anchor pool are divided into multiple grid cells. The target filling quantity of multiple grid cells is allocated based on the number of molecules filled. Based on the target filling quantity and the target mask value, the anchors in the grid cells are filled by the molecules to be filled. During the filling process, the molecules to be filled are randomly rotated. If the distance between the molecule to be filled and the already placed atom is less than the corresponding minimum distance threshold, the current placement is abandoned. The boundary layer anchor pool is filled first. The output module is used to calculate the absolute spatial coordinates of all atoms after all molecules have been filled, output a coordinate file, and complete the molecular dynamics modeling of the target porous medium.

Images