A molecular dynamics simulation system

Abstract

The present invention relates to molecular dynamics simulation, and particularly to a molecular dynamics simulation system. The CPU host side reads the system topology file and the system parameter file by using the memory mapping file technology, and saves the system topology file in the form of an array structure (Structure of Arrays, SoA). After receiving the data copied from the CPU host side, the DCU / GPU device side performs molecular simulation through a simulation channel, and the simulation channel is a std::vector of an ensemble base class, and this data records the identifier of each segment of ensemble operation. The sequence of execution of each ensemble in one step is the list builder, the velocity controller, the position controller, the particle force controller, the velocity controller, and the temperature / pressure controller. Compared with the prior art, the present invention has the advantage of high compatibility and can stably run on domestic systems, domestic CPUs, and domestic DCUs, while being compatible with non-domestic systems, CPUs, and GPUs.

Description

Technical Field

[0001] The present invention relates to molecular dynamics simulation, and particularly to a molecular dynamics simulation system. Background Art

[0002] Currently, in the design and development process of molecular dynamics software on the market, the differences in domestic chip toolchains have not been fully considered, resulting in many versions being unable to adapt to the specific requirements of domestic compilers and being difficult to directly compile. Such software usually has complex dependencies, involving not only the code of the software itself but also a large number of third-party libraries. Due to these intricate dependencies, it is often necessary to modify the software itself and the third-party code it depends on during the transplantation process, increasing the difficulty and workload of transplantation.

[0003] Compared with the existing framework, the framework architecture proposed by the present invention is more concise and has fewer dependencies, thus greatly reducing the risk brought by the incompatibility of third-party libraries with domestic devices. In addition, the framework of the present invention is specifically designed for the domestic platform and is compatible with the compilation toolchain of domestic chips. In terms of algorithm design, the present invention makes full use of the characteristics of domestic chips (dynamic allocation of device memory adapts to the memory rate being inferior to foreign GPUs, and the neighbor list construction algorithm and calculation of short-range forces for the 64-thread warp of the DCU are adapted), giving play to the advantages of the device. Summary of the Invention

[0004] In order to be able to perform molecular dynamics simulation using domestic chips, the present invention proposes a molecular dynamics simulation system, including a CPU host side and a DCU / GPU device side, wherein:

[0005] The CPU host side reads the system topology file and the system parameter file using the memory mapping file technology, and saves the system topology file in the form of an array structure (Structure of Arrays, SoA).

[0006] After receiving the data copied from the CPU host side, the DCU / GPU device side performs molecular simulation through a simulation channel, and the simulation channel is a std::vector of an ensemble base class, which records the identifier of each segment of ensemble operation; for each ensemble, within one time step, the list builder, velocity controller, position controller, particle force controller, velocity controller, and temperature / pressure controller execute actions in sequence.

[0007] Further, when using the memory - mapped file technology for reading, the types of files read at least include basic - type files, charge - type files, and topology - type files. Among them, the basic - type file at least includes the total number of particles required for simulation, the total number of particle types, the three - dimensional box size information of the simulation system, force - field parameters, particle index arrays, particle type arrays, particle three - dimensional coordinate arrays, and particle initial velocity arrays. The charge - type file adds a particle charge array on the basis of the basic - type file. The topology - type file adds topological information composed of multiple particles on the basis of the charge - type file. Each type of file is read in segments. After file mapping, the ReadHeader function is configured to automatically identify and parse the file - header information and establish an index structure.

[0008] Further, when the CPU host stores data, the data on the CPU host is stored using the dynamic array std::vector of the C++ standard library; the data on the DCU device is stored using thrust::device_vector built into Hygon DTK; similarly, the same data structure is used for AMD GPUs, that is, thrust::device_vector built into Cuda is used for NVIDIA GPUs; for the data transferred from the CPU host to the DCU / GPU device, the thrust::uninitialized_copy method in the MemoryScheduler base class is used for copying.

[0009] Further, the velocity controller, position controller, particle force controller, velocity controller, and temperature / pressure controller are divided onto the threads of the DCU / GPU device according to the number of simulated particles. At the same time, for the operators that need to handle neighbor relationships, a number of threads equal to the number of warps assigned to a single particle are used to handle its neighbors simultaneously.

[0010] Further, the list builder builds a neighbor list for a particle, specifically including:

[0011] Create a neighbor - list object, that is, assign particles to their respective cells by the linked - cell method and assign a cell ID to each particle. This cell ID is the ID of the cell to which the particle is assigned. Then, use thrust::make_tuple to encapsulate all the information of each particle into a tuple. Next, apply thrust::make_zip_iterator to process the encapsulated tuple so that it can be traversed one by one. Then, perform the thrust::stable_sort_by_key operation to sort all the tuples and their cell - ID arrays according to the cell ID.

[0012] Calculate the start and end indices of the particles in each cell. Use the kernel function of DTK / Cuda to calculate the start and end positions of each cell in the particle list. That is, in the present invention, particles are assigned to each cell, and each cell stores the start index and end index of multiple particles. One index position in each cell corresponds to the index information of one particle. The number of indices in each cell is related to the number of particles filled in the cell. The sequence number of the first stored particle in the cell is used as its start index, and the sequence number of the last stored particle is used as its end index;

[0013] For each particle, assign a warp to it. All threads in the warp search for particles with a distance less than the cut-off distance in the configuration file in 27 cells of the current particle as neighbor particles;

[0014] For the number of neighbor particles calculated by each thread in the warp, use the WarpReduce function in the hipcub library of the DTK function interface or the cub library of the Cuda function interface to sum them to obtain the total number of neighbor particles of the current particle, and calculate the maximum number of neighbor particles of the particle;

[0015] Accumulate the maximum number of neighbor particles of all particles through the DeviceReduce::Sum function in the hipcub library of the DTK function interface or the cub library of the Cuda function interface to calculate the total number of particles included in the neighbor list;

[0016] Then use the calculated total number of particles as the input parameter of the resize method of the device_vector in the thrust library to complete the allocation of the memory of the neighbor list;

[0017] Add the neighbors of the particle to its neighbor list in turn.

[0018] Furthermore, the process of calculating the maximum number of neighbor particles of the current particle according to the total number of neighbor particles of the current particle includes:

[0019]

[0020] Among them, represents finding the maximum value, represents the total number of neighbor particles of the current particle; is the storage space reserved for the neighbors of each atom; represents the floor operation; represents the number of threads in a warp; represents the minimum storage space reserved for atom neighbors.

[0021] Further, the process of sequentially adding the neighbors of a particle to its neighbor list includes:

[0022] Assign a warp to each particle. For the particles within the cut-off distance of the current particle, instead of directly recording the number of neighbor particles, use a flag bit is_neighbor to mark whether it is a neighbor of the current particle. If it is a neighbor of the current particle, the value of is_neighbor is 1, otherwise it is 0;

[0023] Perform a prefix sum calculation on is_neighbor of each thread through the WarpScan function in the hipcub library of the DTK function interface or the cub library of the Cuda function interface to obtain an offset, ensuring that each thread fills the index of the neighbor particle into the correct position in the neighbor list;

[0024] Broadcast the value of offset + is_neighbor of the last thread within the warp through the ShuffleIndex function in the hipcub library of the DTK function interface or the cub library of the Cuda function interface. Accumulate this value within each warp loop, and use this value to update the total number of neighbor particles of the corresponding particle. Update the write position for the next update of its warp using the total number of neighbor particles of the particle and its write position. The write position of each particle is initialized to the starting index of the particle;

[0025] During the simulation process, if the number of neighbors of a certain particle exceeds the maximum number of neighbor particles, modify the should_realloc flag through an atomicOr atomic operation to re-trigger the estimation to dynamically adjust the neighbor list.

[0026] Further, when calculating the short-range interaction force of particles, traverse and calculate the short-range interaction force of particles according to the start and end index arrays of neighbors in the neighbor particle list of the particles; when calculating the long-range interaction force of particles, use the thrust::reduce operator of DTK / Cuda for multi-threaded parallel calculation.

[0027] The present invention is designed based on the characteristics of domestic chips. Compared with the prior art, it has the advantage of high compatibility, can stably run on domestic systems, domestic CPUs, and domestic DCUs, and is also compatible with non-domestic systems, CPUs, and GPUs. The present invention uses the MMAP technology, stores data in a structured object array, streamlines the architecture to organize DCU / GPU operators, and designs an efficient neighbor list construction method and a particle force calculation method, giving full play to the advantages of domestic devices and achieving accurate molecular dynamics simulation. Compared with the prior art, the present invention has the following beneficial effects:

[0028] 1. The present invention realizes fast data reading through memory - mapped file technology, reduces the I / O operation time, and improves the data loading efficiency.

[0029] 2. The present invention uses a structured object array to replace the traditional array structure, which can significantly improve the data cache hit rate and memory access efficiency. It is particularly suitable for data - intensive tasks. This optimization method is particularly obvious for improving the performance of domestic DCU devices, can better adapt to their hardware characteristics, and give full play to the computing potential.

[0030] 3. The present invention abstracts the entire simulation process into a simulation pipeline. The simulation pipeline is a combination of ensembles. Each ensemble contains multiple controllers, and each controller contains several DCU / GPU operators, and may also contain CPU operators adapted to the domestic chip instruction set. This architecture is convenient for transplantation and expansion on domestic hardware platforms and reduces complexity at the same time.

[0031] 4. The present invention combines the built - in operators of DTK / Cuda and uses the Linked Cell method to dynamically allocate neighbor list memory, thus significantly improving the neighbor search efficiency. The Linked Cell method divides the space into several cells and only searches for particles within adjacent cells, greatly reducing the amount of calculation. In particular, this method performs well on domestic DCU devices. In the case of limited memory rate, the present invention effectively alleviates the memory bandwidth bottleneck by optimizing the memory allocation strategy and further improves the overall computing performance. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 It is an architecture diagram of a molecular dynamics simulation system of the present invention;

[0033] Figure 2 It is a working flow diagram of a molecular dynamics simulation system of the present invention;

[0034] Figure 3 It is a potential energy comparison diagram between the present invention and lammps software;

[0035] Figure 4 It is an RDF comparison diagram between the present invention and lammps software. DETAILED DESCRIPTION OF THE INVENTION

[0036] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0037] The present invention provides a molecular dynamics simulation system, including a CPU host side and a DCU / GPU device side, where:

[0038] The CPU host side reads the system topology file and the system parameter file by using the memory-mapped file technology, and saves the system topology file as an array in the SoA form;

[0039] After receiving the data copied from the CPU host side, the DCU / GPU device side performs molecular simulation through the computing unit - simulation channel. The simulation channel is a std::vector of an ensemble base class, and this data records the identifier of each segment of ensemble operation; for each ensemble, the list builder, velocity controller, position controller, particle force controller, velocity controller, and temperature / pressure controller execute actions in sequence within one time step.

[0040] Such as Figure 1 , a molecular dynamics simulation system of the present invention mainly includes a CPU host side and a DCU / GPU device side. The present invention maps the system topology file and the system parameter file to the memory of the CPU host side through an optimized memory-mapped file (Memory-Mapped File, MMAP) technology, and then saves the system topology file as an array in the SoA form on the CPU host side; the data saved as an array in the SoA form is copied to the DCU / GPU device side and stored in the storage unit of the DCU / GPU device side, and then the molecular dynamics simulation calculation is performed in the computing - simulation pipeline of the DCU / GPU device side.

[0041] Such as Figure 2 , in this embodiment, a schematic diagram of the working process of a molecular dynamics simulation system is given, which specifically includes the following steps:

[0042] S1. The MMAP reads the topology, coordinate, and force field parameter files.

[0043] Molecular dynamics software usually needs to input a topology file containing particle coordinates, a simulation configuration file, etc. to set the simulation conditions. In order to significantly accelerate the reading process of these files, the present invention adopts an optimized memory-mapped file technology and makes innovative improvements and adaptive modifications in the following aspects:

[0044] (1) A unified and efficient reading mechanism for multiple types of files. The present invention supports reading three types of files: basic type files, charge type files, and topology type files. The basic type files at least include the total number of particles required for simulation, the total number of particle types, the three-dimensional box size information of the simulation system, force field parameters, particle index arrays, particle type arrays, particle three-dimensional coordinate arrays, and particle initial velocity arrays. The charge type files add a particle charge array based on the basic type files. The topology type files add topological information composed of multiple particles based on the charge type files, such as topological information such as bonding lists, angle lists, and dihedral lists composed of particles. By designing a unified memory mapping and parsing mechanism, the system can automatically call the corresponding parsing function according to the file type, realize efficient and flexible reading of multiple types of files, and configure the AllocateDataSpace function. This function dynamically allocates the corresponding data structure according to atom_style (basic, charge, topology) to meet the needs of different file types.

[0045] (2) Segment mapping and parallel processing optimization: In view of the specific structure of molecular dynamics files, the present invention maps the files into logical segments. Each segment corresponds to a different part of the file (such as Masses, Force-field parameters, Atoms, bond-list and angle-list, etc.). Multi-threading is used to parse different data segments in parallel, making full use of the computing power of multi-core processors, significantly improving the data parsing speed and reducing the overall reading time.

[0046] (3) Dynamic type identification and pre-parsing index establishment: After the file is mapped, configure the ReadHeader function to automatically identify and parse the file header information and establish the necessary index structure.

[0047] S2. Create a SoA storage data structure and data copy suitable for DCU / GPU devices.

[0048] When the files required for simulation are input into the software, the simulation data will be stored in the form of SoA, which is a memory layout method that stores data of the same type in a centralized manner, which is conducive to improving cache hit rate and accelerating data processing. In this embodiment, the host data is stored using the dynamic array std::vector of the C++ standard library, and the device data is stored using thrust::device_vector built into Haiguang DTK. This data structure is also used for AMD GPUs, and this data structure in Cuda is used for NVIDIA GPUs.

[0049] In molecular dynamics, different types of data need to be transmitted when applying different force fields, and the data transmission mechanisms are also different. Therefore, the present invention constructs a base class MemoryScheduler to adapt to the basic data of different force fields, making the data transmission more reusable. For transmitting host data to the device, the thrust::uninitialized_copy method in this base class is used for copying, and this copy operation does not initialize the values of the allocated memory.

[0050] S3. Create a simulation pipeline that supports multi-segment ensemble simulation.

[0051] After the data storage or copying to the device is completed, a simulation pipeline (SimulatePipeline) will be initialized. The simulation pipeline generates a simulation process according to the configuration file input by the user. Specifically, the simulation pipeline is a std::vector of an ensemble base class, and this vector records the identifiers (flags) of each segment of ensemble operation. Such an operation can sequentially configure different ensembles for molecular dynamics simulation in different stages, and the ensemble contains a series of controllers that implement the integration to solve Newton's equations of motion. In a simulation pipeline, a neighbor list builder, a velocity controller, a position controller, a velocity controller, a particle force controller, and a temperature / pressure controller are sequentially executed within a time step. All the calculations of the controllers in the simulation pipeline of the present invention are implemented on the computing units of the DCU / GPU. The computing units of the DCU / GPU contain several blocks, and each block contains several threads. Each block schedules the number of threads in a warp at a time.

[0052] All the controllers of the present invention are divided onto the threads of the DCU / GPU according to the number of simulated particles. At the same time, for the operators that need to handle neighbor relationships (such as neighbor list construction and calculation of short-range forces), the number of threads in a warp is allocated to each particle to handle its neighbors simultaneously.

[0053] S4. Create a neighbor list that supports adaptation.

[0054] In this embodiment, the process of creating a neighbor list includes:

[0055] I. Construct a neighbor list object, including the following steps:

[0056] 1. First, use the linked cell method to allocate particles to their respective cells and assign a cell ID to each particle;

[0057] 2. Then, use thrust::make_tuple to encapsulate all information such as the position and velocity of each particle into a tuple;

[0058] 3. Subsequently, apply thrust::make_zip_iterator to process these tuples so that they can be traversed one by one;

[0059] 4. Finally, perform the thrust::stable_sort_by_key operation to sort all tuples and their cell ID arrays according to the cell ID.

[0060] Through the above operations, the data is arranged tightly and orderly in memory, making it easier to achieve aligned memory access (continuously accessing n bytes optimized by hardware), and improving the cache hit rate.

[0061] II. Calculate the start and end indices of particles in each cell, specifically including: In the present invention, according to the number of cells, each particle is evenly distributed to each cell, and the number of indices of the particles assigned to each cell is included in each cell. Since the cell ID array to which each particle belongs has been sorted, the start and end positions of each cell in the particle list can be calculated using the kernel function of DTK / Cuda. Specifically, if M particles occupy one cell, every M sorted particles are assigned a cell ID. If a particle has a different cell ID from its previous particle, then this particle is the start index of the current cell. Similarly, if a particle has a different cell ID from its subsequent particle, then this particle is the end index of the current cell. The start position and end position of each cell in the particle sequence can be obtained through the kernel function of DTK / Cuda, thereby significantly reducing the number of searches and improving the construction efficiency.

[0062] III. Estimate the space occupancy of the neighbor list, specifically including the following steps:

[0063] The present invention designs a neighbor list builder base class, which can be derived into a full neighbor list builder and a half neighbor list builder. In the full neighbor list, for any two adjacent particles, not only the connection relationship from particle A to particle B is recorded, but also the connection relationship from particle B to particle A is recorded; while in the half neighbor list, only the connection relationship from particle A to particle B is recorded, and the connection relationship from particle B to particle A is not recorded.

[0064] The system of the present invention uses the neighbor list builder to perform a complete neighbor estimation without storing the actual neighbor list. Specifically, for each particle, warpSize (the size of a thread block) threads are assigned to it, and the warpSize threads are looped to simultaneously search for neighbor particles with a distance less than the cut-off distance in the configuration file in the surrounding 27 cells (one layer in each of the xyz directions), and the number of neighbor particles maintained by the current thread is increased.

[0065] For each thread in the warp, the number of neighbor particles calculated is summed up using WarpReduce of hipcub or cub to obtain the total number of neighbor particles (atom_neighbor_num) of the current particle, and then calculate:

[0066]

[0067] Among them, represents the maximum value, represents the total number of neighbor particles of the current particle; is the storage space reserved for the neighbors of each atom; represents the floor operation; represents the number of threads in a warp; represents the minimum storage space reserved for atom neighbors.

[0068] In the present invention, warpSize is a built-in value of the device; BUFFER_ZONE is a configurable constant that reserves a certain amount of space for the number of neighbors of each atom to accommodate possible newly added neighbor atoms during the simulation. The empirical value of BUFFER_ZONE is often taken as 1.2. If the value of BUFFER_ZONE is too small, it will be more likely to trigger memory reallocation, and if the value is too large, it will cause waste of memory space. At the same time, for atoms with the number of neighbors less than the minimum threshold MIN_NBNUM, the number of their neighbors is set to MIN_NBNUM. This threshold is designed as a multiple of the device warp size to ensure memory alignment and allow sparser regions in the simulation domain to receive more atoms during the simulation, making the estimated space more flexible. The introduction of the MIN_NBNUM threshold can also reduce the number of times of dynamically allocating memory space.

[0069] In this embodiment, the maximum number of neighbors of particles is accumulated through DeviceReduce::Sum of hipcub or cub to calculate the total number of particles contained in the neighbor list (denoted as nb_list_num), and then this value is used as the input parameter of the resize method of device_vector in the thrust library to complete the allocation of neighbor list memory.

[0070] Then, use DeviceScan::ExclusiveSum of hipcub or cub to perform a scan operation on the maximum number of neighbors of the particles, so as to calculate the starting index position of each particle's neighbors in the neighbor list. Here, the scan operation refers to accumulating the maximum number of neighbors of the previous particle to generate an accumulated sequence, and each element of this sequence represents the starting index of the current particle's neighbors in the list. At the same time, for each particle, use the starting index plus its number of neighbor particles (atom_neighbor_num) to obtain its ending index in the neighbor list, and then set the flag should_realloc to false to indicate that there is no need to re-estimate and adjust the memory to prepare for subsequent operations.

[0071] IV. Fill the neighbor list. This step is similar to the estimation process. Similarly, each particle is assigned warpSize (thread block size) threads. However, for particles within the truncation distance, instead of directly recording the number of neighbor particles, the flag is_neighbor is set to 1 to mark them, and for non-neighbor particles outside the truncation distance, the flag is the default value 0. Specifically, it includes the following steps:

[0072] First, perform a prefix sum calculation on is_neighbor of each thread through WarpScan of hipcub or cub. Specifically, use the starting index start of the current particle in the neighbor list plus the offset as the write position of the neighbor. The offset is the prefix sum of the is_neighbor flag bits, that is, judge whether to write to the neighbor list according to the flag bit. If it is necessary, the value of the flag bit is set to 1, otherwise it is set to 0. For example, whether a certain particle in the i-th thread cell of the current particle's thread block is a neighbor of the current particle. If it is not within the intercept range of the current particle, it is not a neighbor of the current particle, and the value of the is_neighbor flag bit is 0 at this time. If it is within the intercept range of the current particle, it is a neighbor of the current particle, and the value of the is_neighbor flag bit is 1 at this time, and the neighbor particle needs to be written into the neighbor list of the current particle. The write position is star+offset(i), where star represents the starting index of the neighbor list of the current particle, and offset(i) represents the offset value of the current particle, and its value is the value obtained by performing a prefix sum calculation on is_neighbor of the i-th thread, ensuring that each thread fills the index of the neighbor particle into the correct position in the neighbor list. As a specific implementation method, if a thread block includes 4 threads, and the sequence of is_neighbor values of each thread is represented as [1, 0, 1, 1], then the sequence of prefix sums is represented as [0, 1, 1, 2], so there is:

[0073] Thread 0: is_neighbor = 1, offset = 0, the writing position of this neighboring particle is start + 0;

[0074] Thread 1: is_neighbor = 0, this neighboring particle is not written;

[0075] Thread 2: is_neighbor = 1, offset = 1, the writing position of this neighboring particle is start + 1;

[0076] Thread 3: is_neighbor = 1, offset = 2, the writing position of this neighboring particle is start + 2;

[0077] The total number of neighboring particles is: the offset value of the last thread + the is_neighbor value of the last thread = 2 + 1 = 3.

[0078] Next, broadcast the offset + is_neighbor value of the last thread within the warp through hipcub or cub's ShuffleIndex, and accumulate this value within each warp loop to update the current atomic total number of neighboring particles. Use the total number of neighboring particles of the particle and its writing position to update the writing position for the next update of its warp. The writing position of each particle is initialized to the starting index of the particle;

[0079] During the simulation process, if the number of neighbors of a certain particle exceeds the maximum number of neighboring particles, modify the should_realloc flag through the atomicOr atomic operation to re-trigger the estimation to dynamically adjust the neighbor list.

[0080] As an alternative implementation, in this embodiment, the user can configure a safety distance in the configuration file. During the process of updating the neighbor list of the particles, extend the original cut-off distance based on the safety distance to obtain a new neighbor cut-off distance. Generally, when using the original cut-off distance to obtain neighbor information, it is necessary to update the neighbor information at each step, and the overall time consumption will increase. When using the safety distance to obtain the new neighbor cut-off distance, it is only necessary to update the neighbor information every 20 steps or every 100 steps (those skilled in the art can adjust this number of steps according to the size of the simulation system), and the overall time consumption is greatly reduced. In addition, this safety distance can be large or small. A larger safety distance will reduce the update frequency of the neighbor list and the time consumption, but at the same time, more neighbor information will be stored, resulting in an increase in memory occupancy. A smaller safety distance reduces the memory occupied by storing neighbor information, but at the same time, the update frequency will increase. Generally speaking, configuring the safety distance to 2 Å is a good choice.

[0081] S5. Use Transform to sum the particle interaction forces.

[0082] In molecular dynamics simulations, the calculation of the short-range interaction forces of particles often involves the neighbor information of the particles. Therefore, use the neighbor start and end index arrays of the particles prepared in step S4 to traverse and obtain the neighbor particles j of each central particle i, and calculate the relative coordinate r of the central particle i and the neighbor particle j ij (r ij = r j -r i ) and then calculate the short-range interaction forces of the particles. Such a data structure can greatly reduce the time consumption required for calculating the forces; while in molecular dynamics simulations, the long-range interaction forces of particles often involve the local summation information of the particles, such as the local summation of the particle structure factor under each wave vector in k-space. The traditional method is to traverse and sum through a for loop under the total number of wave vectors. When the number of particles and the simulation space are large, the total number of wave vectors will be relatively large. Therefore, the traditional method will be more time-consuming. To reduce the time consumption required for the long-range interaction forces, use the thrust::reduce operator of DTK / Cuda for multi-threaded parallel summation of the structure factor. To sum the contributions of the total interaction forces of the particles in the molecular dynamics simulation force field, the present invention uses the thrust::transform operator of DTK / Cuda, which supports different types of data and operation functions. The transform uses a kernel function to assign a unique thread index to each particle at the bottom layer, and then parallelly sums the total interaction forces corresponding to each particle. If the force field developer develops a new particle interaction force, only need to pass its force array into the thrust::transform operator, and it can be quickly summed, which is efficient and convenient. In addition, use the AtomicAdd operator of DTK / Cuda to ensure the total energy contributed by all particles is summed under thread safety.

[0083] S6. Update the particle velocity and position.

[0084] Use the total force of each particle obtained in step S5 to update the position and velocity of each particle. This step also uses the kernel function of DTK / Cuda to parallelly process the position and velocity of each particle.

[0085] S7. Determine whether the iteration is over. If not, return to step S4. Otherwise, transfer the DCU / GPU data back to the host for post-processing of the data.

[0086] In this embodiment, a water molecule system (H2O) is selected as an example of this method to illustrate the present invention, which specifically includes the following steps:

[0087] 101. Prepare the system topology file and the system parameter file.

[0088] First, prepare the H2O molecular topology file. This file is of the topology type and contains the following information:

[0089] (1) 2000 H2O molecules, a total of 6000 particles; the dimensions of the three-dimensional orthogonal box of the system (the side length of x is 30 Å, the side length of y is 50 Å, and the side length of z is 38 Å);

[0090] (2) The force field parameters of water molecules; (3) The index arrays of oxygen atoms and hydrogen atoms in water molecules, as well as the corresponding coordinate arrays, initial velocity arrays, and charge arrays; the bonding list and the angle list between oxygen atoms and hydrogen atoms.

[0091] Use the MMAP method in step S1 to read the H2O molecular topology file.

[0092] Secondly, prepare the configuration file for running the water molecule system. This file contains the following information:

[0093] (1) The storage path of the H2O molecular topology file;

[0094] (2) The storage path of the force field type and its force field parameter or potential function parameter file;

[0095] (3) The simulation ensemble;

[0096] (4) The total number of running steps.

[0097] Also use the MMAP method to read this configuration file and store it in memory at once for subsequent simulations.

[0098] 102. Save the system topology file in the form of a SoA array, and use the MemoryScheduler module to copy the topology data and force field parameters of water molecules into the thrust::device_vector array type recognizable by the DCU / GPU device.

[0099] 103. According to the simulation ensemble information of the water molecules read, use step S3 to initialize a SimulatePipeline "simulation pipeline", and the data of the water molecule system starts to flow into the simulation pipeline.

[0100] 104. Use step S4 to construct a neighbor index array for each particle of the water molecule system, and calculate the relative coordinate r of each particle ij .

[0101] 105. The relative coordinate r of the particle ijAnd the incoming steps such as the charge array are input into the particle force interaction controller in step S5. In this embodiment, the CVFF force field is selected. The calculation of this force field includes contributions such as particle short-range force, particle long-range force, and particle bond angle force. The kernel function is used for efficient parallel computing inside various forces, and the atomicAdd operator is used to obtain the energy of the water molecule system. Finally, the thrust::transform operator is used to sum these forces to obtain the total force of the particles, and the position and velocity of the particles are updated iteratively through the total force of the particles.

[0102] The time step of this embodiment is set to 0.5 fs, the cut-off distance of the neighbor list is set to 12 Å, the ensemble is set to the NVT ensemble, the temperature is maintained at 298 K, and a total of 30,000 steps are iteratively run. After the operation is completed, the particle coordinate array (trajectory data) and the velocity array on the DCU / GPU device are transmitted back to the host memory for subsequent calls to the data post-processing controller.

[0103] To prove the effectiveness of the present invention, under the same force field and initial water molecule model settings, the data post-processing results obtained by the present invention are compared with the molecular dynamics software LAMMPS for accuracy. The comparison of the energy (total potential energy) is as Figure 3 shown. The total potential energies of water molecules obtained by the present invention and the LAMMPS software are -2.460 kcal / mol and -2.466 kcal / mol respectively, and the relative error is very small, which is 0.2%. The comparison results of the radial distribution function RDF obtained by the present invention and the LAMMPS software are as Figure 4 shown. The first peak of RDF(O-O) is both located at 2.75 Å, and the second peak of RDF(O-O) is also both located at 4.45 Å. Figure 3 and Figure 4 The simulation shows that the framework structure provided by the present invention has good accuracy when applied to molecular dynamics simulation.

[0104] Although the embodiments of the present invention have been shown and described, it will be understood by those of ordinary skill in the art that various changes, modifications, substitutions and variations can be made to these embodiments without departing from the principles and spirit of the present invention. The scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A molecular dynamics simulation system, characterized in that, It includes the CPU host side and the DCU / GPU device side, including: The CPU host uses memory mapping file technology to read the system topology file and system parameter file, and saves the system topology file in the form of an array of SoA; After receiving the data copied from the CPU host, the DCU / GPU device performs molecular simulation through the simulation channel. The simulation channel is a std::vector of the ensemble base class. The data records the identifier of each ensemble operation. In each ensemble, the list builder, velocity controller, position controller, particle force controller, velocity controller, and temperature / pressure controller perform actions in sequence within a time step. The process of adding the neighbors of the particle to its neighbor list in sequence includes: Create a neighbor list object, that is, assign particles to their respective cells through the linked cell method and assign a cell ID to each particle, then use the thrust::make_tuple operator to encapsulate all the information of each particle into a tuple, then apply the thrust::make_zip_iterator operator to process the encapsulated tuple so that it can be traversed one by one, and then execute the thrust::stable_sort_by_key operator operation to sort all tuples and their cell ID arrays according to the cell ID; Calculate the starting and ending indexes of the particles in each cell, and use the DTK / Cuda kernel function to calculate the starting and ending positions of each cell in the particle list; For each particle, a thread warp is assigned to it, and all threads in the thread warp search for particles whose distance is less than the cutoff distance in the configuration file in the 27 cells of the current particle as neighbor particles; The number of neighbor particles calculated by each thread in the thread bundle is summed using the WarpReduce function of the hipcub function library in the DTK function interface or the cub function library in the Cuda function interface to obtain the total number of neighbor particles of the current particle, and the maximum number of neighbors of each particle is calculated; The DeviceReduce::Sum function of the hipcub function library in the DTK function interface or the cub function library in the Cuda function interface is used to accumulate the maximum number of neighbors of all particles and calculate the total number of particles included in the neighbor list; Then use the calculated total number of particles as the input parameter of the resize method of the device_vector function in the thrust function library to complete the allocation of neighbor list memory; Add the neighbors of the particle to its neighbor list one by one.

2. The molecular dynamics simulation system according to claim 1, wherein When reading using the memory-mapped file technology, the types of files read at least include basic type files, charge type files, and topology type files. Among them, the basic type files at least include the total number of particles required for the simulation, the total number of particle types, the three-dimensional box size information of the simulation system, force field parameters, particle index arrays, particle type arrays, particle three-dimensional coordinate arrays, and particle initial velocity arrays. The charge type files add a particle charge array on the basis of the basic type files. The topology type files add topological information composed of multiple particles on the basis of the charge type files. Each type of file is read in segments. After the file is mapped, the ReadHeader function is configured to automatically identify and parse the file header information and establish an index structure.

3. A molecular dynamics simulation system according to claim 1, characterized in that, When storing data on the CPU host side, the data on the CPU host side is stored using the dynamic array std::vector of the C++ standard library; the data on the DCU device side is stored using the thrust::device_vector built into the Haiguang DTK; the same data structure is used for the AMD GPU, and the data structure in Cuda is used for the NVIDIA GPU; for the data transferred from the CPU host side to the DCU / GPU device side, the thrust::uninitialized_copy method in the MemoryScheduler base class is used for copying.

4. A molecular dynamics simulation system according to claim 1, characterized in that, The velocity controller, position controller, particle force controller, velocity controller, and temperature / pressure controller are divided among the threads on the DCU / GPU device side according to the number of simulated particles. At the same time, for the operators that need to handle neighbor relationships, a number of threads equal to the number of warps assigned to a single particle are used to handle its neighbors simultaneously.

5. A molecular dynamics simulation system according to claim 1, characterized in that, The calculation of the maximum number of neighbors of a particle includes: ； Among them, represents finding the maximum value, represents the total number of neighboring particles of the current particle; is the storage space reserved for the neighbors of each atom; represents the floor operation; represents the number of threads in a warp; represents the minimum storage space reserved for atomic neighbors.

6. A molecular dynamics simulation system according to claim 1, characterized in that, The process of adding the neighbors of a particle to its neighbor list in sequence includes: Assign a warp to each particle. For the particles within the cut-off distance of the current particle, the number of neighbor particles is not directly recorded. Instead, the flag bit is_neighbor is used to mark whether it is a neighbor of the current particle. If it is a neighbor of the current particle, the value of is_neighbor is 1, otherwise it is 0; The prefix sum calculation of is_neighbor for each thread is performed through the WarpScan function in the hipcub library of the DTK function interface or the cub library of the Cuda function interface to obtain the offset, ensuring that each thread fills the index of the neighbor particle into the correct position in the neighbor list; Broadcast the value of offset + is_neighbor of the last thread within a warp through the ShuffleIndex function of the hipcub library in the DTK function interface or the cub library in the Cuda function interface. Accumulate this value within each warp loop, and use this value to update the total number of neighbor particles of the corresponding particle. Update the write position for the next update of its warp using the total number of neighbor particles of the particle and its write position. Initialize the write position of each particle to the starting index of the particle; During the simulation, if the number of neighbors of a certain particle exceeds the maximum number of neighbor particles, modify the should_realloc flag through an atomicOr atomic operation and re-trigger the estimation to dynamically adjust the neighbor list.

7. A molecular dynamics simulation system according to claim 1, characterized in that, When calculating the tomographic interaction force of particles, traverse and calculate the short-range interaction force of particles according to the start and end index arrays of neighbors in the neighbor particle list of the particles; when calculating the long-range interaction force of particles, use the thrust::reduce operator of the DTK function interface or the Cuda function interface to obtain it in multi-threaded parallel.

Images