Robust/src/Benchmarks/Scheduling/JGFMolDyn/MD.java

   1 /** Banboo Version */
   2
   3 /**************************************************************************
   4  *                                                                         *
   5  *             Java Grande Forum Benchmark Suite - Version 2.0             *
   6  *                                                                         *
   7  *                            produced by                                  *
   8  *                                                                         *
   9  *                  Java Grande Benchmarking Project                       *
  10  *                                                                         *
  11  *                                at                                       *
  12  *                                                                         *
  13  *                Edinburgh Parallel Computing Centre                      *
  14  *                                                                         *
  15  *                email: epcc-javagrande@epcc.ed.ac.uk                     *
  16  *                                                                         *
  17  *                  Original version of this code by                       *
  18  *                         Dieter Heermann                                 *
  19  *                       converted to Java by                              *
  20  *                Lorna Smith  (l.smith@epcc.ed.ac.uk)                     *
  21  *                   (see copyright notice below)                          *
  22  *                                                                         *
  23  *      This version copyright (c) The University of Edinburgh, 2001.      *
  24  *                         All rights reserved.                            *
  25  *                                                                         *
  26  **************************************************************************/
  27
  28 public class MD {
  29     flag initialise;
  30     flag move;
  31     flag fire;
  32     flag update;
  33     flag scale;
  34     flag validate;
  35     flag finish;
  36
  37     public int mm;
  38     public int group;
  39     public int mdsize;
  40     public int move;
  41     public int movemx;
  42
  43     public float side;
  44     float hsq, hsq2, den, h, vaver,vaverh;
  45     float tref, tscale;
  46     int irep;
  47     int istop;
  48     int iprint;
  49
  50     public Particle[] one;
  51     public float[][] sh_force;
  52     //public float[] epot;
  53     //public float[] vir;
  54     //public float[] ek;
  55     public float epot;
  56     public float vir;
  57     public float ek;
  58     float ekin;
  59
  60     public int counter;
  61
  62     public void MD(int d, int g) {
  63         this.mm = d;
  64         this.group = g;
  65         this.mdsize = this.mm * this.mm * this.mm * 4;
  66         this.one = new Particle [this.mdsize];
  67         this.sh_force = new float[3][this.mdsize];
  68         this.move = 0;
  69         this.movemx = 2;
  70         this.den = (float)0.83134;
  71         this.h = (float)0.064;
  72         this.side = Math.powf((float)(this.mdsize/this.den),(float)0.3333333);
  73         this.hsq = this.h * this.h;
  74         this.hsq2 = this.hsq * (float)0.5;
  75         this.vaver = (float)1.13 * Math.sqrtf((float)this.tref / (float)24.0);
  76         this.vaverh = this.vaver * this.h;
  77         this.irep = 10;
  78         this.istop = 19;
  79         //this.epot = (float)0.0;
  80         //this.vir = (float)0.0;
  81         this.tref = (float)0.722;
  82         this.tscale = (float)16.0 / ((float)1.0 * (float)this.mdsize - (float)1.0);
  83         this.iprint = 10;
  84         this.counter = 0;
  85
  86         this.epot = (float)0.0;//new float[this.group + 1];
  87         this.vir = (float)0.0;//new float[this.group + 1];
  88         this.ek = (float)0.0;//new float[this.group + 1];
  89     }
  90
  91     public void init() {
  92         /*for(int i = 0; i < this.group + 1; i++) {
  93             this.epot[i] = (float)0.0;
  94             this.vir[i] = (float)0.0;
  95         }*/
  96         this.epot = (float)0.0;
  97         this.vir = (float)0.0;
  98         for(int j = 0; j < 3; j++) {
  99             for (int i = 0; i < this.mdsize; i++) {
 100                 this.sh_force[j][i] = (float)0.0;
 101             }
 102         }
 103         this.counter = 0;
 104     }
 105
 106     public void initialise() {
 107
 108         /* Particle Generation */
 109         float xvelocity, yvelocity, zvelocity;
 110         int ijk, lg, i, j, k;
 111         float a = this.side/this.mm;
 112         xvelocity = (float)0.0;
 113         yvelocity = (float)0.0;
 114         zvelocity = (float)0.0;
 115         //System.printString("here 1\n");
 116         //System.printI(0xa0);
 117         ijk = 0;
 118         for (lg=0; lg<=1; lg++) {
 119             for (i=0; i<mm; i++) {
 120                 for (j=0; j<mm; j++) {
 121                     for (k=0; k<mm; k++) {
 122                         one[ijk] = new Particle((float)(i*a+lg*a*(float)0.5),(float)(j*a+lg*a*(float)0.5),(k*a),
 123                                 xvelocity,yvelocity,zvelocity,this);
 124                         ijk = ijk + 1;
 125                     }
 126                 }
 127             }
 128         }
 129         for (lg=1; lg<=2; lg++) {
 130             for (i=0; i<mm; i++) {
 131                 for (j=0; j<mm; j++) {
 132                     for (k=0; k<mm; k++) {
 133                         one[ijk] = new Particle((float)(i*a+(2-lg)*a*(float)0.5),(float)(j*a+(lg-1)*a*(float)0.5),
 134                                 (float)(k*a+a*(float)0.5),xvelocity,yvelocity,zvelocity,this);
 135                         ijk = ijk + 1;
 136                     }
 137                 }
 138             }
 139         }
 140         //System.printString("here 2\n");
 141         //System.printI(0xa1);
 142         /* Initialise velocities */
 143         int iseed;
 144         float v1,v2, r;
 145         iseed = 0;
 146         v1 = (float)0.0;
 147         v2 = (float)0.0;
 148
 149         MyRandom random = new MyRandom(iseed,v1,v2);
 150         //System.printString("here 3\n");
 151         //System.printI(0xa2);
 152         for (i=0; i<this.mdsize; i+=2) {
 153             r  = random.seed();
 154             one[i].xvelocity = r*random.v1;
 155             one[i+1].xvelocity  = r*random.v2;
 156         }
 157         //System.printString("here 4\n");
 158         //System.printI(0xa3);
 159         for (i=0; i<this.mdsize; i+=2) {
 160             r  = random.seed();
 161             one[i].yvelocity = r*random.v1;
 162             one[i+1].yvelocity  = r*random.v2;
 163         }
 164         //System.printString("here 5\n");
 165         //System.printI(0xa4);
 166         for (i=0; i<this.mdsize; i+=2) {
 167             r  = random.seed();
 168             one[i].zvelocity = r*random.v1;
 169             one[i+1].zvelocity  = r*random.v2;
 170         }
 171
 172         /*for(i = 0; i < this.mdsize; i++) {
 173                 System.printString("xvel: " + (int)(one[i].xvelocity*100000) + "; yvel: " + (int)(one[i].yvelocity*100000) + "; zvel: " + (int)(one[i].zvelocity*100000) + "\n");
 174         }*/
 175
 176         //System.printString("here 6\n");
 177         //System.printI(0xa5);
 178         /* velocity scaling */
 179         float sp, ts, sc;
 180         ekin = (float)0.0;
 181         sp = (float)0.0;
 182         //System.printString("here 7\n");
 183         //System.printI(0xa6);
 184         for(i=0;i<this.mdsize;i++) {
 185             sp = sp + one[i].xvelocity;
 186         }
 187         sp = sp / this.mdsize;
 188         //System.printString("here 8\n");
 189         //System.printI(0xa7);
 190         for(i=0;i<this.mdsize;i++) {
 191             one[i].xvelocity = one[i].xvelocity - sp;
 192             ekin = ekin + one[i].xvelocity*one[i].xvelocity;
 193         }
 194         //System.printString("here 9\n");
 195         //System.printI(0xa8);
 196         sp = (float)0.0;
 197         for(i=0;i<this.mdsize;i++) {
 198             sp = sp + one[i].yvelocity;
 199         }
 200         sp = sp / this.mdsize;
 201         //System.printString("here 10\n");
 202         //System.printI(0xa9);
 203         for(i=0;i<this.mdsize;i++) {
 204             one[i].yvelocity = one[i].yvelocity - sp;
 205             ekin = ekin + one[i].yvelocity*one[i].yvelocity;
 206         }
 207         //System.printString("here 11\n");
 208         //System.printI(0xa10);
 209         sp = (float)0.0;
 210         for(i=0;i<this.mdsize;i++) {
 211             sp = sp + one[i].zvelocity;
 212         }
 213         sp = sp / this.mdsize;
 214         //System.printString("here 12\n");
 215         //System.printI(0xa11);
 216         for(i=0;i<this.mdsize;i++) {
 217             one[i].zvelocity = one[i].zvelocity - sp;
 218             ekin = ekin + one[i].zvelocity*one[i].zvelocity;
 219         }
 220         //System.printString("here 13\n");
 221         //System.printI(0xa12);
 222         ts = this.tscale * ekin;
 223         sc = h * Math.sqrtf(this.tref/ts);
 224
 225         for(i=0;i<this.mdsize;i++) {
 226             one[i].xvelocity = one[i].xvelocity * sc;
 227             one[i].yvelocity = one[i].yvelocity * sc;
 228             one[i].zvelocity = one[i].zvelocity * sc;
 229         }
 230         //System.printString("here 14\n");
 231         //System.printI(0xa13);
 232     }
 233
 234
 235     public void domove(){
 236         for (int i=0;i<this.mdsize;i++) {
 237             one[i].domove(this.side,i);
 238         }
 239
 240         /*for(int j=0;j<3;j++) {
 241             for (int i=0;i<mdsize;i++) {
 242                 sh_force[j][i] = (float)0.0;
 243             }
 244         }*/
 245
 246         this.move++;
 247     }
 248
 249     public boolean finish() {
 250         if(this.move == this.movemx) {
 251             return true;
 252         }
 253         return false;
 254     }
 255
 256     public void update(MDRunner runner) {
 257         float sum, vel,velt, count, sc;
 258         float etot,temp,pres,rp;
 259
 260         /* update force arrays */
 261         for(int k=0;k<3;k++) {
 262             for(int i=0;i<this.mdsize;i++) {
 263                 sh_force[k][i] += runner.sh_force2[k][i];
 264             }
 265         }
 266
 267         //runner.init();
 268
 269         //this.epot[runner.id + 1] = runner.epot;
 270         //this.vir[runner.id + 1] = runner.vir;
 271         this.epot += runner.epot;
 272         this.vir += runner.vir;
 273     }
 274
 275     public void sum() {
 276         float sum = (float)0.0;
 277
 278         for (int j=0;j<3;j++) {
 279             for (int i=0;i<this.mdsize;i++) {
 280                 sh_force[j][i] = sh_force[j][i] * hsq2;
 281             }
 282         }
 283
 284         sum = (float)0.0;
 285         for (int i=0;i<this.mdsize;i++) {
 286             sum = sum + this.one[i].mkekin(hsq2,i);
 287         }
 288         ekin = (float)(sum/hsq);
 289     }
 290
 291     public void scale() {
 292         float sum, vel,velt, count, sc;
 293         float etot,temp,pres,rp;
 294
 295         //runner.epot = this.epot[0];
 296         //runner.vir = this.vir[0];
 297
 298         vel = (float)0.0;
 299         count = (float)0.0;
 300
 301         /* average velocity */
 302         for (int i=0;i<this.mdsize;i++) {
 303             velt = this.one[i].velavg(vaverh,h);
 304             if(velt > vaverh) {
 305                 count = count + (float)1.0;
 306             }
 307             vel = vel + velt;
 308         }
 309         vel = (float)(vel / h);
 310
 311         /* temperature scale if required */
 312         if((this.move < this.istop) && (((this.move+1) % this.irep) == 0)) {
 313             sc = Math.sqrtf(this.tref / (this.tscale*ekin));
 314             for (int i=0;i<mdsize;i++) {
 315                 one[i].dscal(sc,1);
 316             }
 317             ekin = (float)(this.tref / this.tscale);
 318         }
 319
 320         /* sum to get full potential energy and virial */
 321         if(((this.move+1) % this.iprint) == 0) {
 322             /*this.ek[runner.id+1] = (float)24.0*ekin;
 323             this.epot[runner.id+1] = (float)4.0*this.epot[runner.id+1];
 324             etot = this.ek[runner.id+1] + this.epot[runner.id+1];
 325             temp = this.tscale * ekin;
 326             pres = den * (float)16.0 * (ekin - this.vir[runner.id+1]) / mdsize;
 327             vel = vel / this.mdsize;
 328             rp = (count / this.mdsize) * (float)100.0;
 329
 330             if(this.counter == this.group) {*/
 331                 this.ek = (float)24.0*ekin;
 332                 //this.epot[0] = (float)4.0*this.epot[0];
 333                 //etot = this.ek[0] + this.epot[0];
 334                 //temp = this.tscale * ekin;
 335                 //pres = den * (float)16.0 * (ekin - this.vir[0]) / mdsize;
 336                 //vel = vel / this.mdsize;
 337                 //rp = (count / this.mdsize) * (float)100.0;
 338                 //System.printString("ek: " + (int)(this.ek*1000000) + " (" + this.move + ")\n");
 339             //}
 340         }
 341     }
 342
 343     public void validate() {
 344         float refval = (float)1731.4306625334357;
 345         float dev = Math.abs(this.ek - refval);
 346         if (dev > 1.0e-10 ){
 347             //System.printString("Validation failed\n");
 348             //System.printString("Kinetic Energy = " + (int)(this.ek*1000000) + ";  " + (int)(refval*1000000) + "\n");
 349             //System.printI(0xdddf);
 350             //System.printI((int)(this.ek[0]*1000000));
 351             //System.printI((int)(refval*1000000));
 352             //System.printI((int)(ek*10000));
 353             //System.printI((int)(refval*10000));
 354         }
 355     }
 356 }