gfs@abvax.UUCP (Greg F. Shay) (04/05/88)
(I posted the source to comp.sources.misc. If it does not get through, ask me for a posting here or I can email to you. ...Greg) With regard to the posting about a practical partical system, here is a simple 2-d system I wrote for experimentation purposes a while ago that illustrates a few ideas in speeding up the computations and beating the O(n*n) computation complexity. The system is a collection of simple body atoms that exhibit a force function on each other characteristic of a actual atoms: I use a piece of a cubic function to model this. Note beyond a distance D there is no interaction. At closer than D, atoms first attract (like surface tension phenomena), and then when too close, there is repulsion because the atoms cannot occupy the same space. | + | + (repulsion) | + Force ^ | + | + | + | + D -----+-------------------++++----- Distance between atoms -> + + + + + + (attraction) + The program is rough, but it does work. The first program computes the positions of the atoms for every time frame and outputs them into a file. The second reads the file and animates the data on a graphics screen. If the number of atoms is small, you can pipe the two together and compute in real time (I have even run the programs on different machines and ethernetted the intermediate data live. If the computation gets too slow for good animation precompute and save the intermediate binary file. The display program is presently written for a Masscomp system. The basic calls are: mgibox which is used draw a solid box of 2x2 pixels. Mgifb is flipping frame buffers at the next vertical retrace, and mgiclear is used to erase the back (non visible) frame buffer. Adaptations can be made to other graphics output hardware. If anyone writes a display routine for a SUN, let me have a copy please, I have not the time to do it myself. The program computes forces, divides by mass to get acceleration, integrates to get velocity, again to get distance. The input file specifies the distance threshold D, mass, dt (delta time), number of iterations, and number of atoms, in that order. See the the program for more info. If the time is too long or mass too small (forces are fixed), integration errors will get too large and the system will diverge (REALLY diverge, try it, set dt = 0.5, BOOM!). A very interesting input atom set is a square lattice of atoms spaced about 0.25 or so apart (5 x 5 or so), the lattice will oscillate symmetrically before errors creep in and make the mass amorphous. You can adjust the display program to magnify just the area where you put your lattice to show the details better. A VERY interesting modification I made once was to put in a 'cooling' function, where the reacting force was set to say 0.9700 of the sum of forces. This causes the kinetic energy to gradually bleed off. A cluster of atoms will eventually cool into a regular lattice, a crystal! Face centered cubic no less! Just like we learned in chemistry 101! Now, the comments on computational speedup. All the atoms are kept sorted by there y coordinate. An index array is used to make swapping of floats unneccessary. Making use of the distance D, for a given atom, all atoms ahead of it on the list by more than a distance D and all the atoms after it on the list by a distance D or more do not have to be considered. Only the atoms in the slab from (y-D,y+D) are considered. Next, the x distance is tested to eliminate more. Finally, dx^2+dy^2, the true distance is calculated only if inside the bounding box. After the forces are summed for all interacting atoms, the velocity and positions are updated. A bubble sort is used because it has the property that if only a few items moved relative to each other, it is fairly quick. I know a better sort would be appropriate, but I find bubble's easy to write. This speedup by trivial elimination is quite visible when running a realtime (computation on the fly) animation. When there are many atoms clustered together in a multi-body interaction, the program slows down. When they are separated, the program speeds up. I would be curious to know how much of the time is spent in the sorting and how much in the floating computations. I use a 20mhz 68020+68881 system and it can perform a 6x6 (36 atom) lattice in real time over the ethernet quite well at 15-20fps or so. Well, I hope someone has fun with this. Drop me a line if anyone gets it running. Greg Shay mandrill ! pyramid !... abvax!gfs decvax !
seibel@cgl.ucsf.edu (George Seibel%Kollman) (04/07/88)
Regarding the generation of neighbor lists in particle simulations, and avoiding O(N**2) complexity, this is quite the big deal among computational physicists and chemists. Here are a couple of refs: Grid-cell technique and some good info on periodic boundary conditions: W.F. van Gunsteren, et. al., "On Searching Neighbors in Computer Simulations of Macromolecular Systems", J. Comp. Chem. Vol 5, 272-279 (1984). A very interesting O(N) neighbor list method: Jay Boris, "A Vectorized 'Near Neighbors' Algorithm of Order N Using a Monotonic Logical Grid", J. Comp. Phys, vol ? 1-20, (1986). If anyone implements the Boris algorithm, let me know. George Seibel, UCSF seibel@cgl.ucsf.edu