A molecular dynamics integration method based on a pipelined parallel computing architecture

Abstract

The application relates to the technical field of high-performance computing, and discloses a molecular dynamics integration method based on a pipeline parallel computing architecture. First, a simulation system is constructed in a GPU; then, atomic coordinates in the simulation system are corrected; next, the simulation system is subjected to cell expansion, cutting, package division and coding operations; after that, the GPU calls special hardware (FPGA / ASIC) for potential energy surface reasoning to solve the atomic potential energy surface in a pipeline mode; finally, the GPU solves the Newton motion equation in parallel based on the data returned by the special hardware. Iterative calculation is continuously carried out according to the above process until the set simulation step number is reached. The application can solve the coprocessor performance bottleneck problem caused by the performance improvement of the special hardware for potential energy surface reasoning, and has the characteristics of high hardware resource utilization and high simulation efficiency.

Description

Technical Field

[0001] This invention relates to the field of high-performance computing technology, and provides a molecular dynamics integration method based on a pipelined parallel computing architecture, which can significantly accelerate the computational efficiency of molecular dynamics. Background Technology

[0002] Investigating the physical and chemical properties of materials through experiments is not only time-consuming and costly, but also often difficult due to the lack of suitable experimental equipment. Molecular dynamics (MD) is a computer simulation method that uses supercomputers and other high-powered platforms to solve for the interaction potential functions and Newton's equations of motion of tens or even hundreds of millions of atoms. Finally, it analyzes the macroscopic physical and chemical properties of materials based on theoretical methods of statistical physics. Users only need to input the initial configuration of the simulation system and the corresponding simulation parameters to simulate the phase space trajectories of atoms, observe and reveal the formation mechanism of material properties at the microscopic scale, and assist in the research and design of related materials. Currently, molecular dynamics has developed into the third important method for studying the structure and properties of molecular systems, after experimental and theoretical methods, and is widely used in physics, chemistry, materials science, and other fields.

[0003] Based on simulation accuracy and efficiency, molecular dynamics methods can be divided into classical molecular dynamics (CMD) based on empirical force fields, first-principles molecular dynamics (AIMD) based on density functional theory, and machine learning molecular dynamics (MLMD) based on machine learning. CMD directly calculates the interactions between different atoms using empirical potential function formulas. This method is simple and can quickly predict atomic potential energy and force information based on the chemical environment of atoms, resulting in high computational efficiency. However, the formulas used in CMD are mostly obtained by fitting experimental results, leading to significant empirical errors and thus lower simulation accuracy. AIMD performs first-principles calculations based on density functional theory, using quantum mechanics to describe electronic behavior and calculating the forces on atomic nuclei through electronic structure analysis, achieving very high theoretical and computational accuracy. This method can accurately reflect the microscopic evolution of the system within the required accuracy threshold and can further accurately predict the macroscopic properties of the system based on relevant theories and methods of statistical physics. However, this method suffers from significant drawbacks such as high computational cost and low computational efficiency. With the development of artificial intelligence technology, machine learning molecular dynamics (MLMD) effectively overcomes the shortcomings of the above two methods. By training a neural network potential function model using data obtained through first-principles calculations and inferring the potential energy surface of atoms based on this model, MLMD can achieve AIMD-level accuracy and approach CMD in computational efficiency.

[0004] Patent application (application number 202210436539.6) proposes a molecular dynamics computation chip design method and system based on neural networks. It uses non-von Neumann architecture hardware, such as a Field Programmable Gate Array (FPGA) or Application Specific Integrated Circuit (ASIC), as the lower-level machine to handle potential energy surface inference calculations, while a Central Processing Unit (CPU) acts as the upper-level machine to handle integration calculations (i.e., solving Newton's equations of motion). A pipelined computing architecture is employed to improve hardware utilization, ultimately enhancing the efficiency of potential energy surface solution and reducing power consumption. With further improvements in the computing performance of the lower-level machine, such as replacing the FPGA with an ASIC with higher clock frequencies and richer storage and computing resources, the inference efficiency of the lower-level machine can be improved by 2-3 orders of magnitude compared to the FPGA. However, because the CPU's integration calculation speed cannot match the speed of the lower-level machine's potential energy surface solution, the lower-level machine's computing resources become idle. In this situation, the computing efficiency of the upper-level machine becomes the bottleneck for improving the overall system's computing power. Summary of the Invention

[0005] To overcome the shortcomings of existing technologies and further improve the computational efficiency of MLMD while achieving AIMD-level computational accuracy, this invention provides a molecular dynamics integration method based on a pipelined parallel computing architecture. This method first uses a Graphical Processing Unit (GPU) to replace the CPU as the host machine to perform molecular dynamics integration operations in parallel, forming a pipelined computing architecture with the slave machine. Then, the simulation system is expanded, segmented, and divided into sub-modules to match the potential energy surface inference speed of the slave machine. The slave machine includes, but is not limited to, dedicated hardware for potential energy surface inference such as FPGA and ASIC. The steps of this method are as follows:

[0006] Step S1: Input the simulation parameters of the simulation system and construct the initial configuration of the system. Based on the element types, crystal structure, and other information of the simulation system, construct the initial coordinates and initial velocity information of all atoms in the simulation system in parallel in the global memory of the GPU, and establish the simulation box. All atoms constructed at this time are called the central atom.

[0007] Step S2: Correct the atomic coordinates in the simulation system to satisfy the periodic boundary condition (PBC). According to formulas (1)-(3), atoms that have moved outside the simulation box during the iteration process are moved into the simulation system from the opposite direction.

[0008]

[0009]

[0010]

[0011] Where, r i =(x i y i , z i Let ) be the Cartesian coordinates of atom i, and L = (L x L y L z Let r' be the length of the simulated box in the x, y, and z directions. i =(x′) i y′ i , z′ i ) represents the new coordinates of atom i after applying PBC.

[0012] Step S3: Perform cell expansion on the simulation system. Periodically replicate the simulation system in six directions in space, expanding the cell around the simulation system. cs The atomic information within the cell is stored in global memory. The atoms newly added during cell expansion are called shell atoms.

[0013] Step S4: Divide the simulation system into blocks. First, calculate the number of blocks in the simulation system according to formulas (4)-(7).

[0014]

[0015]

[0016]

[0017] Nb=Nb x ·Nb y ·Nb z (7)

[0018] Where Nb is the total number of blocks in the simulation system, (Nb x Nb y Nb z The numbers represent the number of blocks in the simulation system in the x, y, and z directions, respectively. `ceil()` performs a rounding up operation.

[0019] Then, the block number of each atom is calculated according to formulas (8)-(11), and these numbers are stored in the global memory. Finally, all atom indices are first arranged in ascending order according to the block number of the corresponding atom; if the block numbers of the atoms corresponding to the atom indices are the same, the atom indices are further arranged in ascending order.

[0020]

[0021]

[0022]

[0023] b i =bx i +by i ·Nb x +bz i ·Nb x ·Nb y (11)

[0024] Among them, b i For block number, (bx i by i bz i ) are the block indices for the x, y, and z directions, respectively.

[0025] Step S5 involves subdividing the simulation system. First, based on the atomic spatial arrangement characteristics of the simulation system, the upper limit of the number of atoms per unit volume is estimated, denoted as ρ. m The calculation formula is as follows:

[0026]

[0027] Where, n c0 v represents the number of atoms in the unit cell of the simulation system, and v represents the volume of the unit cell.

[0028] Let N be the maximum number of atoms that the lower-level machine can calculate at any given time. cm The number of blocks sent to the lower-level machine is calculated according to formulas (13)-(15).

[0029]

[0030]

[0031]

[0032] Among them, (nb x nb y nb z These represent the number of blocks sent to the lower-level machine in the x, y, and z directions, respectively.

[0033] Then, based on the block index, the corresponding atom index is retrieved for each block. From the atom index, the atom coordinates and atom type are extracted. The data is then encoded according to the communication protocol and stored in the buffer. At this point, each atom data sent to the buffer constitutes a small packet. In the z-direction, the atom information of each small packet is stored in the buffer sequentially according to the above process until all small packets in the z-direction have been sent. At this point, all these small packets constitute a medium packet. New medium packets are constructed and processed according to the same process. Each time the data of a medium packet is encoded and stored, the corresponding data is sent to the lower-level machine for calculation, while simultaneously preparing the data for the next medium packet.

[0034] Step S6: The lower-level machine is called to perform inference calculations of the potential energy surface according to the pipeline mode, and the GPU decodes the data returned by the lower-level machine.

[0035] In step S7, the GPU solves Newton's equations of motion in parallel based on the atomic potential energy and force data returned by the lower-level machine, and updates the atomic coordinates and velocities. After updating the atomic potential energy and forces of the simulation system, the GPU starts all computing cores in parallel to update the coordinates and velocities of all atoms according to formulas (16)-(18).

[0036]

[0037]

[0038]

[0039] Where, r i (t), v i (t), F i (t) represents the coordinates, velocity, and force acting on atom i at time t. Δt is the set time step. m i This represents the mass of atom i.

[0040] Step S8: Perform statistical analysis on the simulation data. The GPU transfers data from global memory to the CPU, which then analyzes the simulation data and calculates the properties of the material system based on theoretical methods of statistical physics.

[0041] Compared with existing technologies, this invention uses a GPU with a many-core architecture as the phase space trajectory of the host computer parallel computing simulation system, and calls the lower computer to perform potential energy surface inference according to a pipelined hardware scheduling strategy. This can solve the performance bottleneck problem of the host computer caused by the improvement of the lower computer performance, further improve the utilization rate of hardware resources of the host and lower computers, and improve the simulation efficiency of molecular dynamics. Attached Figure Description

[0042] Figure 1 This is a flowchart of the operation of the present invention;

[0043] Figure 2 This is a schematic diagram of PBC and cell expansion operations in an embodiment of the present invention;

[0044] Figure 3 This is a schematic diagram of the cutting operation in an embodiment of the present invention;

[0045] Figure 4 This is a schematic diagram of the subcontracting operation according to an embodiment of the present invention;

[0046] Figure 5 This is a schematic diagram of the pipeline computing framework used in the embodiments of the present invention. Detailed Implementation

[0047] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0048] This invention relates to a molecular dynamics integration method based on a pipelined parallel computing architecture. For example... Figure 1 As shown, before the loop begins, an initial simulation system is built in the GPU according to the user-defined simulation parameters, and the data is stored in global memory. Then, the iterative solution process begins. First, the atomic coordinates in the simulation system are corrected to satisfy the periodic boundary conditions. Then, the simulation system is expanded, and the system is divided into blocks according to the size of the simulation box. Further, the atomic data is packaged into small packets based on the block index and encoded, and the lower-level machine is scheduled to perform potential energy surface inference in a pipelined manner. Finally, the coordinates and velocities of the atoms are updated based on the forces transmitted back from the lower-level machine. After completing the set number of iterations, the macroscopic properties of the simulation system are analyzed and calculated based on theoretical methods of statistical physics.

[0049] Compared to previous technologies, this invention, while maintaining AIMD-level computational accuracy, can further improve the computational efficiency of MD. The implementation steps are described in detail using the Ar atom system as an example.

[0050] Step S1: Input the simulation parameters of the simulation system and construct the initial configuration of the system in the GPU. Store the relative atomic mass of Ar atoms (39.948) in the GPU's global memory. Set the size of the simulation system to 32 units in each of the x, y, and z directions, and construct the atomic coordinates based on the characteristics of a face-centered cubic unit cell, generating a total of 131,072 atoms. Set the initial temperature to 300K, and call the GPU kernel function to construct the initial velocities of all atoms according to a normal distribution. Simultaneously, construct the simulation box based on the above atomic system, with dimensions L. x L y L z All Set the cell radius Find adjacent cutoff radii

[0051] Step S2: Correct the atomic coordinates in the simulation system to satisfy the periodic boundary condition (PBC). Figure 2 The effect of periodic boundary conditions is illustrated in a two-dimensional plane. Ar atoms numbered 1-4 are explicitly marked in the original simulation box and its eight surrounding periodically expanding boxes. When Ar atom number 1 moves out from the left boundary of the central box, after correcting the atom coordinates, its mirror image, Ar atom number 1 to the right, moves into the central box from the right boundary. For the three-dimensional structure of this embodiment, the specific calculation formula is as follows:

[0052]

[0053]

[0054]

[0055] Where, r i =(x i y i , z i ) represents the Cartesian coordinates of Ar atom i in the embodiment, L = (L x L y L z Let r' be the length of the simulated box in the x, y, and z directions. i =(x′) i y′ i , z′ i ) represents the new coordinates of Ar atom i after applying PBC.

[0056] Step S3: Perform cell expansion operation on the simulation system. Figure 2 This illustrates the cell expansion scenario in a two-dimensional plane. After expanding the simulation box around itself, the cell will expand beyond the simulation box by r. cs Atoms within the range are stored in the GPU's global memory, along with their atomic coordinates and type. The atoms generated by this cell expansion are called shell atoms. Figure 2 The display shows one shell atom numbered 1, one numbered 2, three numbered 3, and one numbered 4 located outside the central simulation box in the gray shaded area.

[0057] Step S4: Divide the simulation system into meshes and cut them into blocks. Figure 3The meshing method is illustrated using a two-dimensional example. The solid frame represents the simulated box, and the area outside the solid frame represents the outer shell atoms generated by cell expansion. Each small square represents a block. For the three-dimensional structure of this embodiment, the above two-dimensional mesh can be analogously transformed into a three-dimensional mesh. First, the number of blocks in the simulation system is calculated according to formulas (4)-(7): Nb x =Nb y =Nb z =34, Nb=39304.

[0058]

[0059]

[0060]

[0061] Nb=Nb x ·Nb y ·Nb z (7)

[0062] Then, the block number of each atom is calculated according to formulas (8)-(11), and these numbers are stored in global memory. Finally, all atom indices are arranged in ascending order according to the block number of the corresponding atom; if the block numbers of the atoms corresponding to the atom indices are the same, the atom indices are further arranged in ascending order.

[0063]

[0064]

[0065]

[0066] b i =bx i +by i ·Nb x +bz i ·Nb x ·Nb y (11)

[0067] Step S5: Subdivide the simulation system. First, based on the simulation system and simulation parameters, set the number of atoms in each block to not exceed ρ. m =16. Here, the lower-level machine N is set. cm =16384. From formulas (13)-(14), it can be seen that only nb needs to be satisfied. x ≤18, nb y ≤18. Here we take nb. x =18, nb y =18.

[0068] Then, based on the block index, the corresponding atom index of each block is retrieved. Atom coordinates and atom types are extracted from the atom index, and the data is encoded according to the communication protocol and stored in the buffer. In the z-direction, the atom information of each small packet is stored in the buffer sequentially according to the above process until all small packets in the z-direction are sent, forming a medium packet. New medium packets are then constructed and processed using the same process. Each time the data of a medium packet is encoded and stored, the corresponding data is sent to the lower-level machine for calculation, while simultaneously preparing the data for the next medium packet. Figure 4 The four medium packages in the embodiment are shown, including their package numbers and block numbers.

[0069]

[0070]

[0071]

[0072] Step S6: The lower-level machine is called to perform inference calculations of the potential energy surface according to the pipeline mode, and the GPU decodes the data returned by the lower-level machine. Figure 5 A schematic diagram illustrating the forces acting on Ar atoms in a pipelined update embodiment is shown. The lower-level machine is configured with four parallel units capable of simultaneously and independently performing potential energy surface inference calculations on the Ar atom system.

[0073] In step S7, the GPU solves Newton's equations of motion in parallel based on the forces transmitted back from the lower-level machine, and updates the atomic coordinates and velocities. After updating the atomic potential energy and forces of the entire Ar atom system, the GPU calls all computing cores in parallel to update the coordinates and velocities of all atoms according to formulas (16)-(18).

[0074]

[0075]

[0076]

[0077] Where, r i (t), v i (t), F i (t) represents the coordinates, velocity, and force acting on atom i at time t. Δt is the set time step. m i This represents the mass of atom i.

[0078] Step S8: Perform statistical analysis on the simulation data. The GPU transfers data from global memory to the CPU, which then analyzes the simulation data and calculates the properties of the material system based on theoretical methods of statistical physics.

[0079] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit the scope of protection of the present invention. Referring to the description of these embodiments, those skilled in the art should be able to understand and make relevant modifications or substitutions to the technical solutions of the present invention without departing from the spirit and scope of the present invention.

Claims

1. A molecular dynamics integration method based on a pipelined parallel computing architecture, characterized in that, Includes the following steps: Step S1: Input the simulation parameters of the simulation system and build the initial configuration of the simulation system in the GPU; Step S2: Correct the atomic coordinates in the simulation system so that the simulation system satisfies the periodic boundary conditions; Step S3: Perform cell expansion operation on the simulation system; Step S4: Perform a block segmentation operation on the simulation system; Step S5, perform sub-package operation on the simulation system: First, based on the atomic spatial arrangement characteristics of the simulation system, estimate the upper limit of the number of atoms per unit volume according to formula (1), denoted as , （1） in, To simulate the number of atoms in the unit cell of the system, The volume is the unit cell volume; Record the maximum number of atoms that the lower-level machine can calculate at any given time. The number of blocks sent to the lower-level machine is calculated according to formulas (2)-(4); （2） （3） （4） in, These represent the number of blocks sent to the lower-level machine in the x, y, and z directions, respectively. Then, the atomic index of each block is retrieved based on the block index, and the atomic coordinates and atomic type information are retrieved from the atomic index. Finally, the data is encoded according to the communication protocol and stored in the buffer. At this point, each atomic data sent to the buffer constitutes a small packet. In the z-direction, the atomic information of each small packet is stored in the buffer in sequence according to the above process until all small packets in the z-direction are sent. At this point, all these small packets constitute a medium packet. Following the same process, new medium packets are constructed and processed. Whenever the data of a medium packet is encoded and stored, the data of the medium packet is sent to the lower-level machine for calculation, while preparing the data for the next medium packet. Step S6: The GPU calls the lower-level machine to perform inference calculations on the potential energy surface; Step S7: The GPU solves Newton's equations of motion in parallel, generates simulation data of atomic coordinates and velocities for the next time step, and stores them in the GPU's global memory. Step S8: Perform statistical analysis on the simulation data located in the GPU global memory, and calculate the macroscopic physical properties of the simulation system based on the theoretical methods of statistical physics.

2. The molecular dynamics integration method based on a pipelined parallel computing architecture as described in claim 1, characterized in that: In step S1, the GPU constructs initial atomic coordinates and establishes a simulation box in parallel based on the element types and crystal structure information of the simulation system, and stores the generated initial data in the GPU's global memory.

3. The molecular dynamics integration method based on a pipelined parallel computing architecture as described in claim 1, characterized in that: In step S2, the GPU corrects the atomic data in the global memory in parallel according to formulas (5)-(7), and moves the atoms that have moved outside the simulation box during the iteration process back into the simulation system from the opposite direction, so that all atoms in the simulation system can satisfy the periodic boundary conditions. （5） （6） （7） in, Let i be the Cartesian coordinates of atom i. The length of the simulated box in the x, y, and z directions; The new coordinates of atom i after applying PBC.

4. The molecular dynamics integration method based on a pipelined parallel computing architecture as described in claim 1, characterized in that: The cell expansion operation in step S3 is as follows: the periodically expanded cell is located outside the simulation box, at a distance from the boundary of the simulation box. The atomic information within is stored in global memory as an extension of the simulation system; wherein... The cell expansion length is user-defined and is measured in angstroms.

5. The molecular dynamics integration method based on a pipelined parallel computing architecture as described in claim 1, characterized in that: The block-cutting operation in step S4 is as follows: First, calculate the number of blocks in the simulation system according to formulas (8)-(11). （8） （9） （10） （11） in, This represents the total number of blocks in the entire system. These represent the number of blocks in the entire simulation system in the x, y, and z directions, respectively, and ceil() is the rounding up operation; Then, calculate the block number of each atom according to formulas (12)-(15) and store these numbers in the global memory; finally, arrange all atom indices in ascending order according to the block number of the corresponding atom; if the block numbers of the atoms are the same, further arrange the atom indices in ascending order. （12） （13） （14） （15） in, The block number for atom i. This is the block index of the atom in the x, y, and z directions.

6. The molecular dynamics integration method based on a pipelined parallel computing architecture as described in claim 1, characterized in that: In step S6, the GPU calls the lower-level machine in pipeline mode to perform potential energy surface inference calculation and decodes the data returned by the lower-level machine; the lower-level machine includes, but is not limited to, FPGA and ASIC dedicated hardware for potential energy surface inference.

7. The molecular dynamics integration method based on a pipelined parallel computing architecture as described in claim 1, characterized in that: In step S7, the GPU solves Newton's equations of motion in parallel based on the atomic potential energy and force data returned by the lower-level machine, and updates the atomic coordinates and velocities. After updating the atomic potential energy and forces of the simulation system, the GPU starts all computing cores in parallel to update the coordinates and velocities of all atoms according to formulas (16)-(18). （16） （17） （18） in, , , Let be the coordinates, velocity, and force acting on atom i at time t. For the set time step, This represents the mass of atom i.

8. The molecular dynamics integration method based on a pipelined parallel computing architecture as described in claim 1, characterized in that: In step S8, the GPU transmits the simulation data located in global memory to the CPU, which then analyzes the simulation data and calculates the properties of the material system based on the theoretical methods of statistical physics.

Images