1 /**************************************************************************
3 * Java Grande Forum Benchmark Suite - Thread Version 1.0 *
7 * Java Grande Benchmarking Project *
11 * Edinburgh Parallel Computing Centre *
13 * email: epcc-javagrande@epcc.ed.ac.uk *
16 * This version copyright (c) The University of Edinburgh, 2001. *
17 * All rights reserved. *
19 **************************************************************************/
20 public class JGFMolDynBench {
38 public int interactions;
39 public int[] interacts;
42 public JGFInstrumentor instr;
44 public JGFMolDynBench(int nthreads) {
45 this.nthreads=nthreads;
48 public void JGFsetsize(int size){
52 public void JGFinitialise(){
54 datasizes = new int[2];
59 PARTSIZE = mm*mm*mm*4;
69 public static void JGFapplication(JGFMolDynBench mold) {
71 mold.epot = new double [mold.nthreads];
72 mold.vir = new double [mold.nthreads];
73 mold.ek = new double [mold.nthreads];
74 mold.interacts = new int [mold.nthreads];
76 int partsize, numthreads;
77 partsize = mold.PARTSIZE;
78 numthreads = mold.nthreads;
81 double sh_force2 [][][];
82 sh_force = new double[3][partsize];
83 sh_force2 = new double[3][numthreads][partsize];
88 thobjects = new mdRunner[numthreads];
89 br= new Barrier(numthreads);
91 int[] mid = new int[2];
92 mid[0] = (128<<24)|(195<<16)|(175<<8)|73;
93 mid[1] = (128<<24)|(195<<16)|(175<<8)|69;
96 for(int i=1;i<numthreads;i++) {
97 thobjects[i] = new mdRunner(i,mold.mm,sh_force,sh_force2,br,mold.nthreads,mold);
99 //System.printString("Starting thread "+ i + "\n");
102 //System.printString("Finished starting rest threads\n");
104 thobjects[0] = new mdRunner(0,mold.mm,sh_force,sh_force2,br,mold.nthreads,mold);
106 //System.printString("Starting thread 0\n");
109 //System.printString("Finishing start\n");
111 for(int i=1;i<numthreads;i++) {
112 //System.printString("Joining thread "+ i + "\n");
116 //System.printString("Finished joining all threads\n");
119 public void JGFvalidate(){
120 double[] refval = new double[2];
121 refval[0] = 1731.4306625334357;
122 refval[1] = 7397.392307839352;
123 double dev = Math.fabs(ek[0] - refval[size]);
125 //System.printString("Validation failed\n");
126 //System.printString("Kinetic Energy = " + (long)ek[0] + " " + (long)dev + " " + size + "\n");
131 class mdRunner extends Thread {
134 int id,i,j,k,lg,mdsize,mm;
135 double l,rcoff,rcoffs,side,sideh,hsq,hsq2,vel,velt;
136 double a,r,sum,tscale,sc,ekin,ts,sp;
140 double vaver,vaverh,rand;
141 double etot,temp,pres,rp;
142 double u1,u2,v1,v2,s, xx, yy, zz;
143 double xvelocity, yvelocity, zvelocity;
145 double [][] sh_force;
146 double [][][] sh_force2;
148 int ijk,npartm,iseed,tint;
160 public mdRunner(int id, int mm, double [][] sh_force, double [][][] sh_force2,Barrier br,
161 int nthreads, JGFMolDynBench mymd) {
164 this.sh_force=sh_force;
165 this.sh_force2=sh_force2;
167 this.nthreads = nthreads;
180 //System.printString("Start run method\n");
182 /* Parameter determination */
190 mdsize = mymd.PARTSIZE;
191 one = new particle[mdsize];
195 side = Math.pow((tmpmdsize/tmpden),0.3333333);
202 npartm = tmpmdsize - 1;
203 rcoffs = rcoff * rcoff;
204 tscale = 16.0 / (1.0 * tmpmdsize - 1.0);
205 vaver = 1.13 * Math.sqrt(tref / 24.0);
208 /* Particle Generation */
215 for (lg=0; lg<=1; lg++) {
216 for (i=0; i<mm; i++) {
217 for (j=0; j<mm; j++) {
218 for (k=0; k<mm; k++) {
219 one[ijk] = new particle((i*a+lg*a*0.5),(j*a+lg*a*0.5),(k*a),
220 xvelocity,yvelocity,zvelocity,sh_force,sh_force2,id,this);
227 for (lg=1; lg<=2; lg++) {
228 for (i=0; i<mm; i++) {
229 for (j=0; j<mm; j++) {
230 for (k=0; k<mm; k++) {
231 one[ijk] = new particle((i*a+(2-lg)*a*0.5),(j*a+(lg-1)*a*0.5),
232 (k*a+a*0.5),xvelocity,yvelocity,zvelocity,sh_force,sh_force2,id,this);
239 /* Initialise velocities */
244 randnum = new random(iseed,v1,v2);
246 for (i=0; i<tmpmdsize; i+=2) {
248 one[i].xvelocity = r*randnum.v1;
249 one[i+1].xvelocity = r*randnum.v2;
252 for (i=0; i<tmpmdsize; i+=2) {
254 one[i].yvelocity = r*randnum.v1;
255 one[i+1].yvelocity = r*randnum.v2;
258 for (i=0; i<tmpmdsize; i+=2) {
260 one[i].zvelocity = r*randnum.v1;
261 one[i+1].zvelocity = r*randnum.v2;
265 /* velocity scaling */
270 for(i=0;i<tmpmdsize;i++) {
271 sp = sp + one[i].xvelocity;
275 for(i=0;i<tmpmdsize;i++) {
276 one[i].xvelocity = one[i].xvelocity - sp;
277 ekin = ekin + one[i].xvelocity*one[i].xvelocity;
281 for(i=0;i<tmpmdsize;i++) {
282 sp = sp + one[i].yvelocity;
286 for(i=0;i<tmpmdsize;i++) {
287 one[i].yvelocity = one[i].yvelocity - sp;
288 ekin = ekin + one[i].yvelocity*one[i].yvelocity;
293 for(i=0;i<tmpmdsize;i++) {
294 sp = sp + one[i].zvelocity;
298 for(i=0;i<tmpmdsize;i++) {
299 one[i].zvelocity = one[i].zvelocity - sp;
300 ekin = ekin + one[i].zvelocity*one[i].zvelocity;
304 sc = h * Math.sqrt(tref/ts);
307 for(i=0;i<tmpmdsize;i++) {
309 one[i].xvelocity = one[i].xvelocity * sc;
310 one[i].yvelocity = one[i].yvelocity * sc;
311 one[i].zvelocity = one[i].zvelocity * sc;
315 /* Synchronise threads and start timer before MD simulation */
317 Barrier.enterBarrier(tmpbr);
318 //System.clearPrefetchCache();
323 //TournamentBarrier.enterBarrier(myid, tmpbr);
324 //if (id == 0) JGFInstrumentor.startTimer("Section3:MolDyn:Run", instr.timers);
325 //Barrier.enterBarrier(tmpbr);
329 for (int move=0;move<movemx;move++) {
330 /* move the particles and update velocities */
332 for (i=0;i<tmpmdsize;i++) {
333 one[i].domove(side,i);
337 //System.printString("Barrier #2\n");
338 Barrier.enterBarrier(tmpbr);
339 //System.clearPrefetchCache();
340 //TournamentBarrier.enterBarrier(myid, tmpbr);
344 for (i=0;i<tmpmdsize;i++) {
345 sh_force[j][i] = 0.0;
352 mymd.interacts[id] = 0;
356 //System.printString("Barrier #3\n");
357 Barrier.enterBarrier(tmpbr);
358 //System.clearPrefetchCache();
359 //TournamentBarrier.enterBarrier(myid, tmpbr);
363 for (i=0+id;i<tmpmdsize;i+=nthreads) {
364 one[i].force(side,rcoff,tmpmdsize,i,xx,yy,zz,mymd);
368 //System.printString("Barrier #4\n");
369 Barrier.enterBarrier(tmpbr);
370 //System.clearPrefetchCache();
371 //TournamentBarrier.enterBarrier(myid, tmpbr);
373 /* update force arrays */
375 for(int k=0;k<3;k++) {
376 for(i=0;i<tmpmdsize;i++) {
377 for(j=0;j<nthreads;j++) {
378 sh_force[k][i] += sh_force2[k][j][i];
385 for(int k=0;k<3;k++) {
386 for(i=0;i<tmpmdsize;i++) {
387 for(j=0;j<nthreads;j++) {
388 sh_force2[k][j][i] = 0.0;
395 for(j=1;j<nthreads;j++) {
396 mymd.epot[0] += mymd.epot[j];
397 mymd.vir[0] += mymd.vir[j];
399 for(j=1;j<nthreads;j++) {
400 mymd.epot[j] = mymd.epot[0];
401 mymd.vir[j] = mymd.vir[0];
403 for(j=0;j<nthreads;j++) {
404 mymd.interactions += mymd.interacts[j];
409 //System.printString("Barrier #5\n");
410 Barrier.enterBarrier(tmpbr);
411 //System.clearPrefetchCache();
412 //TournamentBarrier.enterBarrier(myid, tmpbr);
416 for (i=0;i<tmpmdsize;i++) {
417 sh_force[j][i] = sh_force[j][i] * hsq2;
426 //System.printString("Barrier #6\n");
427 Barrier.enterBarrier(tmpbr);
428 //System.clearPrefetchCache();
429 //TournamentBarrier.enterBarrier(myid, tmpbr);
431 /*scale forces, update velocities */
433 for (i=0;i<tmpmdsize;i++) {
434 sum = sum + one[i].mkekin(hsq2,i);
442 /* average velocity */
444 for (i=0;i<tmpmdsize;i++) {
445 velt = one[i].velavg(vaverh,h);
446 if(velt > vaverh) { count = count + 1.0; }
452 /* temperature scale if required */
454 if((move < istop) && (((move+1) % irep) == 0)) {
455 sc = Math.sqrt(tref / (tscale*ekin));
456 for (i=0;i<tmpmdsize;i++) {
459 ekin = tref / tscale;
462 /* sum to get full potential energy and virial */
464 if(((move+1) % iprint) == 0) {
465 mymd.ek[id] = 24.0*ekin;
466 mymd.epot[id] = 4.0*mymd.epot[id];
467 etot = mymd.ek[id] + mymd.epot[id];
468 temp = tscale * ekin;
469 pres = tmpden * 16.0 * (ekin - mymd.vir[id]) / tmpmdsize;
470 vel = vel / tmpmdsize;
471 rp = (count / tmpmdsize) * 100.0;
473 //System.printString("Barrier #7\n");
474 Barrier.enterBarrier(tmpbr);
475 //System.clearPrefetchCache();
476 //TournamentBarrier.enterBarrier(myid, tmpbr);
479 //System.printString("Barrier #8\n");
480 Barrier.enterBarrier(tmpbr);
481 //System.clearPrefetchCache();
482 //TournamentBarrier.enterBarrier(myid, tmpbr);
483 //if (id == 0) JGFInstrumentor.stopTimer("Section3:MolDyn:Run", instr.timers);
484 //System.printString("End run method\n");
494 public double xcoord, ycoord, zcoord;
495 public double xvelocity,yvelocity,zvelocity;
498 double [][] sh_force;
499 double [][][] sh_force2;
502 public particle(double xcoord, double ycoord, double zcoord, double xvelocity,
503 double yvelocity,double zvelocity,double [][] sh_force,
504 double [][][] sh_force2,int id,mdRunner runner) {
506 this.xcoord = xcoord;
507 this.ycoord = ycoord;
508 this.zcoord = zcoord;
509 this.xvelocity = xvelocity;
510 this.yvelocity = yvelocity;
511 this.zvelocity = zvelocity;
512 this.sh_force = sh_force;
513 this.sh_force2 = sh_force2;
518 public void domove(double side,int part_id) {
520 xcoord = xcoord + xvelocity + sh_force[0][part_id];
521 ycoord = ycoord + yvelocity + sh_force[1][part_id];
522 zcoord = zcoord + zvelocity + sh_force[2][part_id];
524 if(xcoord < 0) { xcoord = xcoord + side; }
525 if(xcoord > side) { xcoord = xcoord - side; }
526 if(ycoord < 0) { ycoord = ycoord + side; }
527 if(ycoord > side) { ycoord = ycoord - side; }
528 if(zcoord < 0) { zcoord = zcoord + side; }
529 if(zcoord > side) { zcoord = zcoord - side; }
531 xvelocity = xvelocity + sh_force[0][part_id];
532 yvelocity = yvelocity + sh_force[1][part_id];
533 zvelocity = zvelocity + sh_force[2][part_id];
537 public void force(double side, double rcoff,int mdsize,int x, double xx, double yy, double zz, JGFMolDynBench mymd) {
543 double rd,rrd,rrd2,rrd3,rrd4,rrd6,rrd7,r148;
544 double forcex,forcey,forcez;
547 rcoffs = rcoff*rcoff;
553 for (int i=x+1;i<mdsize;i++) {
554 xx = this.xcoord - runner.one[i].xcoord;
555 yy = this.ycoord - runner.one[i].ycoord;
556 zz = this.zcoord - runner.one[i].zcoord;
558 if(xx < (-sideh)) { xx = xx + side; }
559 if(xx > (sideh)) { xx = xx - side; }
560 if(yy < (-sideh)) { yy = yy + side; }
561 if(yy > (sideh)) { yy = yy - side; }
562 if(zz < (-sideh)) { zz = zz + side; }
563 if(zz > (sideh)) { zz = zz - side; }
566 rd = xx*xx + yy*yy + zz*zz;
575 mymd.epot[id] = mymd.epot[id] + (rrd6 - rrd3);
576 r148 = rrd7 - 0.5*rrd4;
577 mymd.vir[id] = mymd.vir[id] - rd*r148;
581 sh_force2[0][id][i] = sh_force2[0][id][i] - forcex;
586 sh_force2[1][id][i] = sh_force2[1][id][i] - forcey;
591 sh_force2[2][id][i] = sh_force2[2][id][i] - forcez;
593 mymd.interacts[id]++;
598 sh_force2[0][id][x] = sh_force2[0][id][x] + fxi;
599 sh_force2[1][id][x] = sh_force2[1][id][x] + fyi;
600 sh_force2[2][id][x] = sh_force2[2][id][x] + fzi;
604 public double mkekin(double hsq2,int part_id) {
608 xvelocity = xvelocity + sh_force[0][part_id];
609 yvelocity = yvelocity + sh_force[1][part_id];
610 zvelocity = zvelocity + sh_force[2][part_id];
612 sumt = (xvelocity*xvelocity)+(yvelocity*yvelocity)+(zvelocity*zvelocity);
616 public double velavg(double vaverh,double h) {
621 sq = Math.sqrt(xvelocity*xvelocity + yvelocity*yvelocity +
622 zvelocity*zvelocity);
628 public void dscal(double sc,int incx) {
630 xvelocity = xvelocity * sc;
631 yvelocity = yvelocity * sc;
632 zvelocity = zvelocity * sc;
645 public random(int iseed,double v1,double v2) {
651 public double update() {
654 double scale= 4.656612875e-10;
658 int imod = 2147483647;
660 if (iseed<=0) { iseed = 1; }
663 is1 = (iseed-is2)/32768;
666 is1 = (is1*imult+(iss2-is2)/32768) % (65536);
668 iseed = (is1*32768+is2) % imod;
670 rand = scale * iseed;
676 public double seed() {
690 r = Math.sqrt(-2.0*Math.log(s)/s);