A near-field acceleration method for random batch Ewald algorithm

Abstract

The present invention discloses a near-field acceleration method of a random batch Ewald algorithm, and relates to the technical field of molecular dynamics simulation. In the molecular dynamics simulation Ewald algorithm, for the long-range part, the frequency is randomly sampled in the Fourier space to calculate the corresponding force; for the short-range part, a double-layer neighbor list method is used. The present invention uses a random method to approximate the force. In the numerical integration process of the system simulation, the expectation of the random force fluctuation is zero and the variance is controllable, so that the final distribution is close to the actual distribution, so it has good convergence, and the calculation complexity is linear, which is faster than the existing method; only a limited number of frequency samples need to be communicated between CPU cores, no intensive communication is required, and the parallel efficiency is high; the effect of the force within a small cutoff radius is accurately calculated, and the particles in the spherical shell are randomly selected and the force is calculated, which reduces the CPU memory usage and the amount of calculation.

Description

Technical Field

[0001] The invention relates to the technical field of molecular dynamics simulation, in particular to a near-field acceleration method of a random batch Ewald algorithm. Background Art

[0002] In molecular dynamics simulation, the application scenario is generally an overall electrically neutral charged ion system, plus periodic boundary conditions, and the Ewald algorithm is generally used to deal with this problem with long-range interactions. The Ewald algorithm was invented by German physicist Paul Peter Ewald as a method for calculating the electrostatic interaction potential of a periodic system. The core idea is to split the Coulomb interaction between charges into two parts, long-range and short-range. The long-range part is quickly calculated by Fourier transform; the short-range part and the Lennard-Jones potential are both rapidly decaying terms, usually processed using the truncation method; in addition, there are interactions such as covalent bonds and bond angles between atoms. The Lennard-Jones potential was proposed by mathematician John Lennard-Jones to describe the interaction between inert gas molecules. The specific algorithm is as follows:

[0003] The simulation system contains N ions, and the overall electrical neutrality condition is satisfied:

[0004]

[0005] By adding periodic boundary conditions, the interaction potential energy of the whole system can be obtained:

[0006]

[0007] The first three items are the bond energies between particles, representing the bond length, bond angle, and the energy corresponding to the dihedral angle formed by multiple bond angles. This part of the energy is only related to adjacent particles. The last two items represent the Coulomb interaction and the Leonard Jones potential energy, respectively. 3 represents a three-dimensional integer vector, ∑ n ' indicates that the case where i = j and n = 0 is excluded. For Coulomb interaction, it is inaccurate to use the cutoff directly. The classic Ewald method is to use It is divided into two parts: long-distance and short-distance: and Where erfc(x) is the error complement function, which decays rapidly, so the direct truncation method can be used; erf(x) is the error function, which decays slowly, so The Fourier transform is converted into a series with rapid decay in Fourier space, and then the frequency is truncated for calculation. Therefore, the overall force expression for particle i is:

[0008]

[0009] The second term is the force corresponding to the Leonard Jones potential energy, and the structural factor Im means taking the imaginary part, r ij =r j -r i In the Ewald method, choosing an appropriate α and cutoff length k c and r s , making the computational complexity of long-range and short-range interactions both O(N 3 / 2 ).

[0010] The current method reduces the overall computational complexity by applying fast Fourier transform and interpolation to the last force in Fourier space, and can eventually reach O(NlogN). For short-range forces, the current method is to establish a neighbor list for each particle and calculate the overall force generated by the particles in the neighbor list. In addition, the cell list method is also applied to reduce the overall computational complexity to linear; and the neighbor list is considered to be slightly larger than the actual need, so as to reduce the frequency of updating the neighbor list.

[0011] The existing technology has high computational complexity: For the long-range part of Coulomb interaction, the existing method has high complexity; when dealing with short-range forces, the number of interacting particles within the cutoff radius is still large (hundreds), and the computational complexity is high. When dealing with the long-range part of Coulomb interaction, the existing method uses multiple fast Fourier transforms, which consumes a lot of time for information communication during multi-core parallel simulation, reducing parallel efficiency. Each particle needs to store a large list of neighboring particles, which occupies a high amount of CPU memory. When simulating very large particle systems, it will be limited by memory bottlenecks.

[0012] Therefore, technicians in this field are committed to developing a near-field acceleration method (IRBE) of the random batch Ewald algorithm. The computational complexity is linear, which is faster than the existing method; only a limited number of frequency samples need to be communicated between CPU cores, no intensive communication is required, and the parallel efficiency is high; a two-layer neighbor list is constructed, and only the information of the inner neighbor needs to be stored, which reduces the CPU memory usage and the amount of calculation. Summary of the invention

[0013] In view of the above-mentioned defects of the prior art, the technical problem to be solved by the present invention is how to reduce the computational complexity of the Ewald algorithm for molecular dynamics simulation, improve parallel efficiency, and reduce CPU memory usage and computational complexity.

[0014] To achieve the above object, the present invention provides a near-field acceleration method of a random batch Ewald algorithm, and a molecular dynamics simulation of the Ewald algorithm, for the short-range part, a double-layer neighbor list method is used.

[0015] Furthermore, for the long-range part, the frequencies are randomly sampled in Fourier space and the corresponding forces are calculated.

[0016] Furthermore, the frequency is randomly sampled, and the order of magnitude of frequency sample calculations performed during the numerical simulation is hundreds or thousands.

[0017] Furthermore, the frequency is randomly sampled, and global communication between CPU cores is performed only once.

[0018] Furthermore, the double-layer neighbor list method accurately calculates the effect of forces within a small cutoff radius and randomly selects and calculates the forces on particles within the spherical shell.

[0019] Furthermore, the double-layer neighbor list method stores neighbor particle information, including position and speed.

[0020] Furthermore, the double-layer neighbor list only stores information of neighbor particles within a small cutoff radius.

[0021] Furthermore, the double-layer neighbor list randomly selects only a portion of external particles and does not store neighbor information.

[0022] Furthermore, for the short-range interaction, the Coulomb short-range part and the Lennard Jones part are included.

[0023] Furthermore, for short-range interactions, two cutoff radii are selected.

[0024] In a preferred embodiment of the present invention, in the molecular dynamics simulation, the application scenario is an overall electrically neutral charged ion system, with periodic boundary conditions. The prior art Ewald algorithm has a high computational complexity: for the long-range part of the Coulomb interaction, the existing method has a high complexity; when dealing with short-range forces, the number of interacting particles within the cutoff radius is still large (hundreds), and the computational complexity is large. The present invention utilizes random approximate forces. In the numerical integration process of the system simulation, the expectation of the random force fluctuation is zero and the variance is controllable, so that the final distribution is close to the actual distribution, and therefore has good convergence. For the long-range partial force, the frequency is randomly sampled in Fourier space, and the corresponding force is calculated; for the short-range part, the double-layer neighbor list method is considered to accurately calculate the effect of the force within the small cutoff radius, and the particles in the spherical shell are randomly selected and the force is calculated, which can reduce the overall computational complexity. The computational complexity is linear, which is faster than the existing method.

[0025] In terms of multi-core parallelism, the existing technology needs to communicate a lot of information, resulting in a decrease in parallel efficiency: the existing method uses multiple fast Fourier transforms when processing the long-range part of Coulomb interaction, which consumes a lot of time for information communication during multi-core parallel simulation, reducing parallel efficiency. The present invention randomly samples the frequency, and only a limited number (hundreds to thousands) of frequency samples are used for calculation in each numerical simulation process, which only requires one communication between CPU cores, thereby improving parallel efficiency.

[0026] In the prior art, each particle needs to store a large list of neighboring particles, which occupies a high CPU memory. When simulating a very large particle system, it will be limited by the memory bottleneck. The present invention accurately calculates the effect of the force within a small cutoff radius, and randomly selects and calculates the force on particles in the spherical shell. By establishing a double-layer neighbor list, only particles within a small cutoff radius need to be stored. For external particles, only a part needs to be randomly selected, and there is no need to store neighbor information, thereby reducing CPU memory usage.

[0027] Compared with the prior art, the present invention has the following obvious substantial features and significant advantages:

[0028] 1. The force is approximated by random methods. In the numerical integration process of the system simulation, the force errors cancel each other out, making the final distribution close to the actual distribution, so it has good convergence. The computational complexity is linear, which is faster than the existing methods.

[0029] 2. Only a limited number of frequency samples need to be communicated between CPU cores, no intensive communication is required, and the parallel efficiency is high.

[0030] 3. Accurately calculate the force within a small cutoff radius, and randomly select and calculate the force on particles in the spherical shell. Reduce CPU memory usage and calculation amount.

[0031] The concept, specific structure and technical effects of the present invention will be further described below in conjunction with the accompanying drawings to fully understand the purpose, characteristics and effects of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 It is an algorithm flow chart of a preferred embodiment of the present invention;

[0033] Figure 2 is the radial distribution function of oxygen atoms in a water molecule of a preferred embodiment of the present invention;

[0034] Figure 3 is the mean square displacement of the center of mass motion of water molecules in a preferred embodiment of the present invention;

[0035] Figure 4This is the relationship between the computational efficiency of a preferred embodiment of the present invention and the classical algorithm and the number of CPU cores. DETAILED DESCRIPTION

[0036] The following describes several preferred embodiments of the present invention with reference to the drawings in the specification, so that the technical content is clearer and easier to understand. The present invention can be embodied in many different forms of embodiments, and the protection scope of the present invention is not limited to the embodiments mentioned in the text.

[0037] In the drawings, components with the same structure are indicated by the same numerical reference numerals, and components with similar structures or functions are indicated by similar numerical reference numerals. The size and thickness of each component shown in the drawings are arbitrarily shown, and the present invention does not limit the size and thickness of each component. In order to make the illustration clearer, the thickness of the components is appropriately exaggerated in some places in the drawings.

[0038] For the simulated particle system, the most computationally complex part is the calculation of the force on each particle. The overall force expression for particle i is:

[0039]

[0040] The current method can reduce the overall computational complexity to O(NlogN), but for simulating large particle systems, the computational cost is still relatively high. In addition, when multi-core parallel computing is applied, multiple fast Fourier transforms will cause a large amount of communication between CPU cores, reducing computational efficiency.

[0041] To address these issues, stochastic methods are used for both parts. For the long-range part of the Coulomb interaction, we observe that there is a discrete Gaussian distribution term Therefore, we can sample the frequency k in Fourier space according to this distribution, and randomly select one hundred to one thousand frequencies to calculate the force of this part, which can reduce the amount of calculation. On the other hand, the more frequencies are randomly selected, the more accurate the calculation result will be. The force expression corresponding to each frequency is According to the importance sampling principle, assuming that p frequency samples are extracted, the long-range part of the Coulomb effect received by particle i can be calculated as:

[0042]

[0043] It can be calculated that the expectation of the force is equal to the actual force on the particle, thereby reducing the amount of calculation of the long-range force.

[0044] Consider the short-range force, including the Coulomb short-range part and the Lennard-Jones part

[0045]

[0046] The truncation method is generally used, that is, the interaction between two particles that are far apart is not considered, where the truncation radius is r s Before each simulation, it is necessary to create a neighbor list for each particle, determine the numbers of neighbor particles whose distance is less than the cutoff distance, and store the neighbor particle information, including position and velocity. This results in the CPU memory occupied by storage becoming the bottleneck of the simulation for large-scale systems. Therefore, it is proposed to select two stage radii, where r c <r s , thus the particles within the cutoff radius are divided into two parts. The forces generated by the inner particles are calculated explicitly, while the particles in the spherical shell are randomly selected and the random batch principle is used to obtain the corresponding force expression:

[0047]

[0048] Where B(i) represents the particle in the spherical shell, N I is the corresponding number of particles, F cor represents the correction term for the force, which also ensures the conservation of the total momentum of the system.

[0049] Adding up the two forces and the bonding interaction force mentioned above, we get the resultant force received by each particle. Then, according to Newton's second law, the overall system evolves, and the position and velocity of the particles are continuously updated until the simulation stops. Based on the system configuration at each time step and the simulated trajectory of each particle, the physical properties of the simulated system and the dynamic behavior of the particles can be calculated.

[0050] like Figure 1 As shown, the specific calculation steps are as follows:

[0051] 1. Generate the initial state of the system;

[0052] 2. Calculate the long-range and short-range forces on each particle at each time step (as shown in steps 3-5);

[0053] 3. For the long-range part of the Coulomb interaction, randomly select p frequencies in Fourier space and calculate the corresponding particle forces:

[0054]

[0055] 4. For short-range interactions, including the Coulomb short-range part and the Lennard Jones part, two cutoff radii are selected: the large cutoff radius r s Generally, 12 angstroms are selected, and the small cutoff radius r c You can choose 5 angstroms to 7 angstroms;

[0056] 5. Use the small cutoff radius to create a cell and neighbor list. For particles j whose spacing is less than the small cutoff radius, we only need to consider the particles in the 27 adjacent cells. For particles in the spherical shell, we randomly select them and approximately calculate the corresponding forces:

[0057]

[0058]

[0059] 6. Combining the above long-range and short-range parts, as well as the bonding interaction force, the total force on each particle can be calculated;

[0060] 7. According to Newton's second law, time evolution is performed. In addition, for special systems, it is sometimes necessary to add some control items to keep the system temperature or pressure constant.

[0061] Finally, the system configuration at each time step is obtained and the corresponding physical quantities and dynamic properties can be obtained by calculation.

[0062] The computer code of the present invention has been successfully integrated into the all-atom molecular dynamics simulation software LAMMPS, and has been compared with the PPPM algorithm + direct truncation in the all-atom molecular dynamics simulation of a pure water system with 300 million atoms for accuracy and computational efficiency. Comparing the structural information of water molecules, radial distribution function, such as Figure 2 As shown; and dynamic information, mean square displacement, as Figure 3 As shown in Figure 2, the results obtained by this method are completely consistent with those obtained by the traditional PPPM algorithm + direct truncation. At the same time, the computational efficiency of this method (the time required for a dynamic simulation to run one step) is about 10 times higher than that of the PPPM + direct truncation method on 50,000 CPU cores. Figure 4 As shown. At the same time, it can be seen from Figure C that when the number of CPU cores is relatively small, the computational efficiency of this algorithm is similar to that of the classic PPPM algorithm + direct truncation; when the number of CPU cores is more, the computational efficiency of the PPPM + direct truncation method reaches a bottleneck, and the time required for each step of the calculation has not been significantly improved. Only using random batches to accelerate the long-range part will speed up the calculation to a certain extent and improve the parallel efficiency. At this time, the computational efficiency of this method can be further improved, mainly reflected in the acceleration of the short-range part. This means that this method is more suitable for simulating all-atom systems containing a large number of atoms on larger-scale multi-core CPU supercomputers.

[0063] The preferred specific embodiments of the present invention are described in detail above. It should be understood that ordinary technicians in the field can make many modifications and changes based on the concept of the present invention without creative work. Therefore, all technical solutions that can be obtained by technicians in the technical field based on the concept of the present invention through logical analysis, reasoning or limited experiments on the basis of the prior art should be within the scope of protection determined by the claims.

Claims

1. A near-field acceleration method for random batch Ewald algorithm, It is characterized in that Molecular dynamics simulation Ewald algorithm; for the long-range part of the force, the frequency is randomly sampled in Fourier space to calculate the corresponding force; for the short-range part, a double-layer neighbor list method is used to accurately calculate the force within a small cutoff radius, and the particles in the spherical shell are randomly selected and the force is calculated to reduce the computational complexity; In the long-range part, in terms of multi-core parallelism, the frequency is randomly sampled, and only a limited number of frequency samples are used for calculation in each numerical simulation process to improve parallel efficiency; For the short-range part, the truncation method is used. For particles that are far away, the interaction is not considered. The truncation radius is r s ; Create a neighbor list for each particle, determine the number of neighbor particles whose distance is less than the cutoff radius, and store the neighbor particle information, including position and speed; select two stage radii, where r c <r s , the particles within the cutoff radius are divided into two parts; the forces generated by the inner particles are calculated explicitly, and the particles in the spherical shell are randomly selected using the random batch principle.

2. The near-field acceleration method of the random batch Ewald algorithm as claimed in claim 1, It is characterized in that The frequency is randomly sampled, and the order of magnitude of the frequency sample calculations performed during the numerical simulation is hundreds or thousands.

3. The near-field acceleration method of the random batch Ewald algorithm as claimed in claim 1, It is characterized in that The frequency is randomly sampled, and global communication between CPU cores is performed only once.

4. The near-field acceleration method of the random batch Ewald algorithm as claimed in claim 1, It is characterized in that The double-layer neighbor list only randomly selects a part of the external particles and does not store the neighbor information.

5. The near-field acceleration method of the random batch Ewald algorithm as claimed in claim 1, It is characterized in that For short-range interactions, both the Coulomb short-range part and the Leonard Jones part are included.

Images