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 {
31 public DoubleWrapper[] epot;
32 public DoubleWrapper[] vir;
33 public DoubleWrapper[] ek;
38 public int interactions;
39 public IntWrapper[] 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 = global new int[2];
59 PARTSIZE = mm*mm*mm*4;
69 public static void JGFapplication(JGFMolDynBench mold) {
72 int[] mid = new int[8];
73 mid[0] = (128<<24)|(195<<16)|(136<<8)|162; //dw-10
74 mid[1] = (128<<24)|(195<<16)|(136<<8)|163; //dw-11
75 mid[2] = (128<<24)|(195<<16)|(136<<8)|164; //dw-12
76 mid[3] = (128<<24)|(195<<16)|(136<<8)|165; //dw-13
77 mid[4] = (128<<24)|(195<<16)|(136<<8)|166; //dw-14
78 mid[5] = (128<<24)|(195<<16)|(136<<8)|167; //dw-15
79 mid[6] = (128<<24)|(195<<16)|(136<<8)|168; //dw-16
80 mid[7] = (128<<24)|(195<<16)|(136<<8)|169; //dw-17
83 double sh_force2 [][][];
84 int partsize, numthreads;
86 partsize = mold.PARTSIZE;
87 numthreads = mold.nthreads;
88 mybarr = global new BarrierServer(numthreads);
93 sh_force = global new double[3][partsize];
94 sh_force2 = global new double[3][numthreads][partsize];
95 mold.epot = global new DoubleWrapper[numthreads];
96 mold.vir = global new DoubleWrapper[numthreads];
97 mold.ek = global new DoubleWrapper[numthreads];
98 mold.interacts = global new IntWrapper[numthreads];
99 for(int i=0;i<numthreads;i++) {
100 mold.epot[i]=global new DoubleWrapper();
101 mold.vir[i]=global new DoubleWrapper();
102 mold.ek[i]=global new DoubleWrapper();
103 mold.interacts[i]=global new IntWrapper();
108 MDWrap[] thobjects = new MDWrap[numthreads];
111 for(int i=0;i<numthreads;i++) {
112 thobjects[i] = new MDWrap(global new mdRunner(i,mold.mm,sh_force,sh_force2,mold.nthreads,mold));
116 boolean waitfordone=true;
124 for(int i=0;i<numthreads;i++) {
125 thobjects[i].md.start(mid[i]);
128 for(int i=0;i<numthreads;i++) {
129 thobjects[i].md.join();
133 public void JGFvalidate(){
134 double[] refval = new double[2];
135 refval[0] = 1731.4306625334357;
136 refval[1] = 7397.392307839352;
137 double dev = Math.fabs(ek[0].d - refval[size]);
139 //System.printString("Validation failed\n");
140 //System.printString("Kinetic Energy = " + (long)ek[0] + " " + (long)dev + " " + size + "\n");
145 class mdRunner extends Thread {
149 double l,rcoff,rcoffs,side,sideh,hsq,hsq2,vel,velt;
150 double a,r,sum,tscale,sc,ekin,ts,sp;
154 double vaver,vaverh,rand;
155 double etot,temp,pres,rp;
156 double u1, u2, s, xx, yy, zz;
157 double xvelocity, yvelocity, zvelocity;
159 double [][] sh_force;
160 double [][][] sh_force2;
162 int ijk,npartm,iseed,tint;
170 public mdRunner(int id, int mm, double [][] sh_force, double [][][] sh_force2,
171 int nthreads, JGFMolDynBench mymd) {
174 this.sh_force=sh_force;
175 this.sh_force2=sh_force2;
176 this.nthreads = nthreads;
187 public void init(particle[] one, int mdsize) {
189 for (int lg=0; lg<=1; lg++) {
190 for (int i=0; i<mm; i++) {
191 for (int j=0; j<mm; j++) {
192 for (int k=0; k<mm; k++) {
193 one[ijk] = new particle((i*a+lg*a*0.5),(j*a+lg*a*0.5),(k*a),
194 xvelocity,yvelocity,zvelocity,sh_force,sh_force2,id,one);
201 for (int lg=1; lg<=2; lg++) {
202 for (int i=0; i<mm; i++) {
203 for (int j=0; j<mm; j++) {
204 for (int k=0; k<mm; k++) {
205 one[ijk] = new particle((i*a+(2-lg)*a*0.5),(j*a+(lg-1)*a*0.5),
206 (k*a+a*0.5),xvelocity,yvelocity,zvelocity,sh_force,sh_force2,id,one);
213 /* Initialise velocities */
218 random randnum = new random(iseed,v1,v2);
221 for (int i=0; i<mdsize; i+=2) {
223 one[i].xvelocity = r*randnum.v1;
224 one[i+1].xvelocity = r*randnum.v2;
227 for (int i=0; i<mdsize; i+=2) {
229 one[i].yvelocity = r*randnum.v1;
230 one[i+1].yvelocity = r*randnum.v2;
233 for (int i=0; i<mdsize; i+=2) {
235 one[i].zvelocity = r*randnum.v1;
236 one[i+1].zvelocity = r*randnum.v2;
240 /* velocity scaling */
245 for(int i=0;i<mdsize;i++) {
246 sp = sp + one[i].xvelocity;
250 for(int i=0;i<mdsize;i++) {
251 one[i].xvelocity = one[i].xvelocity - sp;
252 ekin = ekin + one[i].xvelocity*one[i].xvelocity;
256 for(int i=0;i<mdsize;i++) {
257 sp = sp + one[i].yvelocity;
261 for(int i=0;i<mdsize;i++) {
262 one[i].yvelocity = one[i].yvelocity - sp;
263 ekin = ekin + one[i].yvelocity*one[i].yvelocity;
268 for(int i=0;i<mdsize;i++) {
269 sp = sp + one[i].zvelocity;
273 for(int i=0;i<mdsize;i++) {
274 one[i].zvelocity = one[i].zvelocity - sp;
275 ekin = ekin + one[i].zvelocity*one[i].zvelocity;
279 sc = h * Math.sqrt(tref/ts);
282 for(int i=0;i<mdsize;i++) {
284 one[i].xvelocity = one[i].xvelocity * sc;
285 one[i].yvelocity = one[i].yvelocity * sc;
286 one[i].zvelocity = one[i].zvelocity * sc;
292 public void doinit(int mdsize) {
293 for(int j=0;j<3;j++) {
294 double[] sh=sh_force[j];
295 for (int i=0;i<mdsize;i++) {
302 public void doinit2(int mdsize) {
303 for(int k=0;k<3;k++) {
304 double[] sh=sh_force[k];
305 double [][] sha=sh_force2[k];
306 for(int j=0;j<nthreads;j++) {
307 double[] sha2=sha[j];
308 for(int i=0;i<mdsize;i++) {
314 for(int k=0;k<3;k++) {
315 double [][] sh1=sh_force2[k];
316 for(int j=0;j<nthreads;j++) {
318 for(int i=0;i<mdsize;i++) {
326 for(int j=1;j<nthreads;j++) {
327 mymd.epot[0].d += mymd.epot[j].d;
328 mymd.vir[0].d += mymd.vir[j].d;
330 for(int j=1;j<nthreads;j++) {
331 mymd.epot[j].d = mymd.epot[0].d;
332 mymd.vir[j].d = mymd.vir[0].d;
334 for(int j=0;j<nthreads;j++) {
335 mymd.interactions += mymd.interacts[j].i;
338 for (int j=0;j<3;j++) {
339 double sh[]=sh_force[j];
340 for (int i=0;i<mdsize;i++) {
341 sh[i] = sh[i] * hsq2;
347 /* Parameter determination */
352 Barrier barr=new Barrier("128.195.175.84");
357 mdsize = mymd.PARTSIZE;
358 one=new particle[mdsize];
361 side = Math.pow((mdsize/tmpden),0.3333333);
369 rcoffs = rcoff * rcoff;
370 tscale = 16.0 / (1.0 * mdsize - 1.0);
371 vaver = 1.13 * Math.sqrt(tref / 24.0);
374 /* Particle Generation */
384 /* Synchronise threads and start timer before MD simulation */
386 Barrier.enterBarrier(barr);
390 for (int move=0;move<movemx;move++) {
392 /* move the particles and update velocities */
393 for (int i=0;i<mdsize;i++) {
394 one[i].domove(side,i);
399 Barrier.enterBarrier(barr);
406 mymd.epot[id].d = 0.0;
407 mymd.vir[id].d = 0.0;
408 mymd.interacts[id].i = 0;
413 Barrier.enterBarrier(barr);
418 for (int i=0+id;i<mdsize;i+=nthreads) {
419 one[i].force(side,rcoff,mdsize,i,xx,yy,zz,mymd);
424 Barrier.enterBarrier(barr);
426 /* update force arrays */
434 Barrier.enterBarrier(barr);
437 /*scale forces, update velocities */
439 for (int i=0;i<mdsize;i++) {
440 sum = sum + one[i].mkekin(hsq2,i);
448 /* average velocity */
450 for (int i=0;i<mdsize;i++) {
451 velt = one[i].velavg(vaverh,h);
452 if(velt > vaverh) { count = count + 1.0; }
458 /* temperature scale if required */
460 if((move < istop) && (((move+1) % irep) == 0)) {
461 sc = Math.sqrt(tref / (tscale*ekin));
462 for (int i=0;i<mdsize;i++) {
465 ekin = tref / tscale;
468 /* sum to get full potential energy and virial */
470 if(((move+1) % iprint) == 0) {
471 mymd.ek[id].d = 24.0*ekin;
472 mymd.epot[id].d = 4.0*mymd.epot[id].d;
473 etot = mymd.ek[id].d + mymd.epot[id].d;
474 temp = tscale * ekin;
475 pres = tmpden * 16.0 * (ekin - mymd.vir[id].d) / mdsize;
477 rp = (count / mdsize) * 100.0;
480 Barrier.enterBarrier(barr);
483 //if (id == 0) JGFInstrumentor.stopTimer("Section3:MolDyn:Run", instr.timers);
490 public double xcoord, ycoord, zcoord;
491 public double xvelocity,yvelocity,zvelocity;
494 global double [][] sh_force;
495 global double [][][] sh_force2;
498 public particle(double xcoord, double ycoord, double zcoord, double xvelocity,
499 double yvelocity,double zvelocity, double [][] sh_force,
500 double [][][] sh_force2,int id, particle[] one) {
502 this.xcoord = xcoord;
503 this.ycoord = ycoord;
504 this.zcoord = zcoord;
505 this.xvelocity = xvelocity;
506 this.yvelocity = yvelocity;
507 this.zvelocity = zvelocity;
508 this.sh_force = sh_force;
509 this.sh_force2 = sh_force2;
514 public void domove(double side,int part_id) {
516 xcoord = xcoord + xvelocity + sh_force[0][part_id];
517 ycoord = ycoord + yvelocity + sh_force[1][part_id];
518 zcoord = zcoord + zvelocity + sh_force[2][part_id];
520 if(xcoord < 0) { xcoord = xcoord + side; }
521 if(xcoord > side) { xcoord = xcoord - side; }
522 if(ycoord < 0) { ycoord = ycoord + side; }
523 if(ycoord > side) { ycoord = ycoord - side; }
524 if(zcoord < 0) { zcoord = zcoord + side; }
525 if(zcoord > side) { zcoord = zcoord - side; }
527 xvelocity = xvelocity + sh_force[0][part_id];
528 yvelocity = yvelocity + sh_force[1][part_id];
529 zvelocity = zvelocity + sh_force[2][part_id];
533 public void force(double side, double rcoff,int mdsize,int x, double xx, double yy, double zz, JGFMolDynBench mymd) {
539 double rd,rrd,rrd2,rrd3,rrd4,rrd6,rrd7,r148;
540 double forcex,forcey,forcez;
543 rcoffs = rcoff*rcoff;
549 for (int i=x+1;i<mdsize;i++) {
550 xx = this.xcoord - one[i].xcoord;
551 yy = this.ycoord - one[i].ycoord;
552 zz = this.zcoord - one[i].zcoord;
554 if(xx < (-sideh)) { xx = xx + side; }
555 if(xx > (sideh)) { xx = xx - side; }
556 if(yy < (-sideh)) { yy = yy + side; }
557 if(yy > (sideh)) { yy = yy - side; }
558 if(zz < (-sideh)) { zz = zz + side; }
559 if(zz > (sideh)) { zz = zz - side; }
562 rd = xx*xx + yy*yy + zz*zz;
571 mymd.epot[id].d += (rrd6 - rrd3);
572 r148 = rrd7 - 0.5*rrd4;
573 mymd.vir[id].d += - rd*r148;
577 sh_force2[0][id][i] = sh_force2[0][id][i] - forcex;
582 sh_force2[1][id][i] = sh_force2[1][id][i] - forcey;
587 sh_force2[2][id][i] = sh_force2[2][id][i] - forcez;
589 mymd.interacts[id].i++;
594 sh_force2[0][id][x] = sh_force2[0][id][x] + fxi;
595 sh_force2[1][id][x] = sh_force2[1][id][x] + fyi;
596 sh_force2[2][id][x] = sh_force2[2][id][x] + fzi;
600 public double mkekin(double hsq2,int part_id) {
604 xvelocity = xvelocity + sh_force[0][part_id];
605 yvelocity = yvelocity + sh_force[1][part_id];
606 zvelocity = zvelocity + sh_force[2][part_id];
608 sumt = (xvelocity*xvelocity)+(yvelocity*yvelocity)+(zvelocity*zvelocity);
612 public double velavg(double vaverh,double h) {
617 sq = Math.sqrt(xvelocity*xvelocity + yvelocity*yvelocity +
618 zvelocity*zvelocity);
624 public void dscal(double sc,int incx) {
625 xvelocity = xvelocity * sc;
626 yvelocity = yvelocity * sc;
627 zvelocity = zvelocity * sc;
636 public random(int iseed,double v1,double v2) {
642 public double update() {
645 double scale= 4.656612875e-10;
649 int imod = 2147483647;
651 if (iseed<=0) { iseed = 1; }
654 is1 = (iseed-is2)/32768;
657 is1 = (is1*imult+(iss2-is2)/32768) % (65536);
659 iseed = (is1*32768+is2) % imod;
661 rand = scale * iseed;
667 public double seed() {
681 r = Math.sqrt(-2.0*Math.log(s)/s);