gfs@decwrl.dec.com@x.UUCP (Greg F. Shay) (04/05/88)
Submitted-By: "Greg F. Shay" <gfs@decwrl.dec.com@x.UUCP> Archive-Name: molecule comp.sources.misc: Volume 2, Issue 92 Submitted-By: "Greg F. Shay" <gfs@decwrl.dec.com@x.UUCP> Archive-Name: molecule [Some people don't read Administrivia, it seems: it's not shar'ed. ++bsa] For a request on comp.graphics: 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 ! /*--------Program molecule.c-------------------------*/ #include <math.h> main() { float dx,dy,d; float forcefunction(); float fx,fy,F; int aindex[100]; int a,b,c; int start,end; int maxindex,itcount,itmax; short i,j,k,l; short idat[2],maxin,itinmax; float mass,dt,thresh; float atompos[100][2]; float atomv[100][2]; float atomf[100][2]; /* Read in initial positions of atoms*/ scanf("%f",&thresh); scanf("%f",&mass); scanf("%f",&dt); scanf("%d\n",&itmax); scanf("%d\n",&maxindex); maxindex = maxindex-1; maxin = (short)maxindex; itinmax = (short)itmax; /*printf("%d \n %d \n",maxindex,itmax);*/ write(1,&maxin,2); write(1,&itinmax,2); for(a=0;a <= maxindex;a++) { for(i=0;i<2;i++) { scanf ("%f ",&atompos[a][i]); } for(i=0;i<2;i++) { scanf ("%f ",&atomv[a][i]); } scanf ("\n"); } /* Initialize index array */ for (a=0;a <= maxindex;a++) aindex[a]=a; /* Now start iteration loop */ for (itcount =0;itcount<itmax;itcount++) { bsort(maxindex,aindex,atompos); for(b=0;b<=maxindex;b++) { idat[0]=(short)(10*atompos[aindex[b]][0]); idat[1]=(short)(10*atompos[aindex[b]][1]); write(1,idat,4); /*printf("%d %d\n",(int)(10*atompos[aindex[b]][0]),(int)(10*atompos[aindex[b]][1])); */ } start = 0; end = 0; for(a=0;a<=maxindex;a++) { /* Do every atom */ /* Adjust start range pointer */ while ( (atompos[aindex[a]][1]-atompos[aindex[start]][1]) > thresh ) start++; /* Adjust end range pointer */ while ( ((atompos[aindex[end]][1]-atompos[aindex[a]][1]) <= thresh) && (end < maxindex) ) end++; /* Evaluate forces on this atom 'a' from every atom between*/ /* start and end on the index list */ atomf[aindex[a]][0]=0.0; atomf[aindex[a]][1]=0.0; for(b=start;b<=end;b++) if (b != a) { dx=atompos[aindex[b]][0]-atompos[aindex[a]][0]; if (dx <= thresh) { dy=atompos[aindex[b]][1]-atompos[aindex[a]][1]; d=sqrt(dx*dx+dy*dy); F=forcefunction(d); fx = F*(dx/d); fy = F*(dy/d); atomf[aindex[a]][0]= fx+atomf[aindex[a]][0]; atomf[aindex[a]][1]= fy+atomf[aindex[a]][1]; } } } /* Do all atoms */ /* Now resolve forces */ for(a=0;a<=maxindex;a++) for(c=0;c<2;c++) { atomv[a][c]= dt*atomf[a][c]/mass+atomv[a][c]; atompos[a][c]=(dt*atomv[a][c])+atompos[a][c]; } /* Now loop back for more iterations */ } /* End of iteration loop */ } /* End of Main */ bsort(limit,index,data) int limit; int index[]; float data[][2]; { short flag,j; int temp; do { flag = 0; for (j=0;j<=limit-1;j++) if (data[index[j]][1] > data[index[j+1]][1]) { flag = 1; temp = index[j+1]; index[j+1]=index[j]; index[j]=temp; } } while (flag == 1); } /* end of bubble sort routine*/ float forcefunction(d) float d; { /* return force depending on distance d */ if (d > 1.0) return (0.); d= (1.0-d); return(60.0*d-160.0*d*d*d); } /*--------End Program molecule.c-------------------------*/ /*-------- Program Display.c-------------------------------*/ main() { int i,j,k,x,y,sx,sy; int fb,maxnum,iterno; short idat[2],maxin,iterin; system("cleargp"); mginit (0,0); mgimodfunc (3,0,3); mgicm(2,4096*15); mgihue(2); fb = 1; read(0,&maxin,2); read(0,&iterin,2); maxnum = (int)maxin; iterno=(int)iterin; /*scanf("%d %d\n",&maxnum,&iterno);*/ for(k=0;k<iterno;k++) { mgifb(fb,3-fb); fb = 3-fb; mgiclearpln(2,-1,0); for(i=0;i<=maxnum;i++) { /*scanf("%d %d\n",&x,&y);*/ read(0,idat,4); mgibox(idat[0],idat[1],idat[0]+2,idat[1]+2); } } mgideagp (); } /*----------_End program display.c --------------------------*/ /* Sample input files */ 1.0 1.0 0.07 1000 10 16 16 0 0 10 10 0 0 14 14 1 1 24 24 0 0 23.5 23.5 0 0 18 18 -.25 .2 20 20 -.1 -.1 32 32 -.4 -.4 28 28 .4 -.3 40 40 -2 -2 /* Sample input files */ 1.0 .2 0.01 800 10 31 31.2 0 0 31 31 0 0 19 19 9 9. 19.5 19.5 9 9. 19.0 19.5 9 9. 31.2 31 0 0 31.2 31.5 0 0 31.4 32 0 0 30.5 31 0 0 30.4 31.5 0 0 /*------------------------------------------*/