public int PARTSIZE;
- public double[] epot;
- public double[] vir;
- public double[] ek;
+ public DoubleWrapper[] epot;
+ public DoubleWrapper[] vir;
+ public DoubleWrapper[] ek;
int size,mm;
int[] datasizes;
public int interactions;
- public int[] interacts;
+ public IntWrapper[] interacts;
public int nthreads;
public JGFInstrumentor instr;
public static void JGFapplication(JGFMolDynBench mold) {
// Create new arrays
- mold.epot = new double [mold.nthreads];
- mold.vir = new double [mold.nthreads];
- mold.ek = new double [mold.nthreads];
- mold.interacts = new int [mold.nthreads];
-
- int partsize, numthreads;
- partsize = mold.PARTSIZE;
- numthreads = mold.nthreads;
+ //BarrierServer mybarr;
+ /*
+ int[] mid = new int[8];
+ mid[0] = (128<<24)|(195<<16)|(175<<8)|84;//dw-10
+ mid[1] = (128<<24)|(195<<16)|(175<<8)|85;//dw-11
+ mid[2] = (128<<24)|(195<<16)|(175<<8)|86;//dw-12
+ mid[3] = (128<<24)|(195<<16)|(175<<8)|87;//dw-13
+ mid[4] = (128<<24)|(195<<16)|(175<<8)|88;//dw-14
+ mid[5] = (128<<24)|(195<<16)|(175<<8)|89;//dw-15
+ mid[6] = (128<<24)|(195<<16)|(175<<8)|90;//dw-16
+ mid[7] = (128<<24)|(195<<16)|(175<<8)|91;//dw-17
+ */
double sh_force [][];
double sh_force2 [][][];
- sh_force = new double[3][partsize];
- sh_force2 = new double[3][numthreads][partsize];
+ int partsize, numthreads;
+ partsize = mold.PARTSIZE;
+ numthreads = mold.nthreads;
+ //mybarr = new BarrierServer(numthreads);
+ //mybarr.start(mid[0]);
+
+ sh_force = new double[3][partsize];
+ sh_force2 = new double[3][numthreads][partsize];
+ mold.epot = new DoubleWrapper[numthreads];
+ mold.vir = new DoubleWrapper[numthreads];
+ mold.ek = new DoubleWrapper[numthreads];
+ mold.interacts = new IntWrapper[numthreads];
+ for(int i=0;i<numthreads;i++) {
+ mold.epot[i]=new DoubleWrapper();
+ mold.vir[i]=new DoubleWrapper();
+ mold.ek[i]=new DoubleWrapper();
+ mold.interacts[i]=new IntWrapper();
+ }
// spawn threads
- mdRunner[] thobjects;
- Barrier br;
- thobjects = new mdRunner[numthreads];
- br= new Barrier(numthreads);
-
- int[] mid = new int[2];
- mid[0] = (128<<24)|(195<<16)|(175<<8)|73;
- mid[1] = (128<<24)|(195<<16)|(175<<8)|69;
- mdRunner tmp;
-
- for(int i=1;i<numthreads;i++) {
- thobjects[i] = new mdRunner(i,mold.mm,sh_force,sh_force2,br,mold.nthreads,mold);
- tmp = thobjects[i];
- //System.printString("Starting thread "+ i + "\n");
- tmp.start();
- }
- //System.printString("Finished starting rest threads\n");
-
- thobjects[0] = new mdRunner(0,mold.mm,sh_force,sh_force2,br,mold.nthreads,mold);
- tmp = thobjects[0];
- //System.printString("Starting thread 0\n");
- tmp.start();
- tmp.join();
- //System.printString("Finishing start\n");
-
- for(int i=1;i<numthreads;i++) {
- //System.printString("Joining thread "+ i + "\n");
- tmp = thobjects[i];
- tmp.join();
- }
- //System.printString("Finished joining all threads\n");
- }
+ MDWrap[] thobjects = new MDWrap[numthreads];
+
+ for(int i=0;i<numthreads;i++) {
+ thobjects[i] = new MDWrap(new mdRunner(i,mold.mm,sh_force,sh_force2,mold.nthreads,mold));
+ }
+
+ /*
+ boolean waitfordone=true;
+ while(waitfordone) {
+ if (mybarr.done)
+ waitfordone=false;
+ }
+ */
+
+ for(int i=0;i<numthreads;i++) {
+ //thobjects[i].md.start(mid[i]);
+ thobjects[i].md.run();
+ }
+ }
public void JGFvalidate(){
double[] refval = new double[2];
refval[0] = 1731.4306625334357;
refval[1] = 7397.392307839352;
- double dev = Math.fabs(ek[0] - refval[size]);
+ double dev = Math.fabs(ek[0].d - refval[size]);
if (dev > 1.0e-10 ){
//System.printString("Validation failed\n");
//System.printString("Kinetic Energy = " + (long)ek[0] + " " + (long)dev + " " + size + "\n");
}
}
-class mdRunner extends Thread {
+class mdRunner {
double count;
- int id,i,j,k,lg,mdsize,mm;
+ int id,i,j,k,lg,mm;
double l,rcoff,rcoffs,side,sideh,hsq,hsq2,vel,velt;
double a,r,sum,tscale,sc,ekin,ts,sp;
double den;
double h;
double vaver,vaverh,rand;
double etot,temp,pres,rp;
- double u1,u2,v1,v2,s, xx, yy, zz;
+ double u1, u2, s, xx, yy, zz;
double xvelocity, yvelocity, zvelocity;
double [][] sh_force;
int istop;
int iprint;
- Barrier br;
- random randnum;
JGFMolDynBench mymd;
int nthreads;
- particle[] one;
-
- public mdRunner(int id, int mm, double [][] sh_force, double [][][] sh_force2,Barrier br,
- int nthreads, JGFMolDynBench mymd) {
+ public mdRunner(int id, int mm, double [][] sh_force, double [][][] sh_force2,
+ int nthreads, JGFMolDynBench mymd) {
this.id=id;
this.mm=mm;
this.sh_force=sh_force;
this.sh_force2=sh_force2;
- this.br=br;
this.nthreads = nthreads;
this.mymd = mymd;
count = 0.0;
iprint = 10;
}
- public void run() {
-
- //System.printString("Start run method\n");
-
- /* Parameter determination */
-
- int tmpmdsize;
- double tmpden;
- int movemx=50;
- Barrier tmpbr;
-
- tmpbr=br;
- mdsize = mymd.PARTSIZE;
- one = new particle[mdsize];
- l = mymd.LENGTH;
- tmpmdsize = mdsize;
- tmpden = den;
- side = Math.pow((tmpmdsize/tmpden),0.3333333);
- rcoff = mm/4.0;
-
- a = side/mm;
- sideh = side*0.5;
- hsq = h*h;
- hsq2 = hsq*0.5;
- npartm = tmpmdsize - 1;
- rcoffs = rcoff * rcoff;
- tscale = 16.0 / (1.0 * tmpmdsize - 1.0);
- vaver = 1.13 * Math.sqrt(tref / 24.0);
- vaverh = vaver * h;
-
- /* Particle Generation */
-
- xvelocity = 0.0;
- yvelocity = 0.0;
- zvelocity = 0.0;
- ijk = 0;
-
- for (lg=0; lg<=1; lg++) {
- for (i=0; i<mm; i++) {
- for (j=0; j<mm; j++) {
- for (k=0; k<mm; k++) {
- one[ijk] = new particle((i*a+lg*a*0.5),(j*a+lg*a*0.5),(k*a),
- xvelocity,yvelocity,zvelocity,sh_force,sh_force2,id,this);
- ijk = ijk + 1;
- }
+ public void init(particle[] one, int mdsize) {
+ int id=this.id;
+ for (int lg=0; lg<=1; lg++) {
+ for (int i=0; i<mm; i++) {
+ for (int j=0; j<mm; j++) {
+ for (int k=0; k<mm; k++) {
+ one[ijk] = new particle((i*a+lg*a*0.5),(j*a+lg*a*0.5),(k*a),
+ xvelocity,yvelocity,zvelocity,sh_force,sh_force2,id,one);
+ ijk = ijk + 1;
}
}
}
+ }
- for (lg=1; lg<=2; lg++) {
- for (i=0; i<mm; i++) {
- for (j=0; j<mm; j++) {
- for (k=0; k<mm; k++) {
- one[ijk] = new particle((i*a+(2-lg)*a*0.5),(j*a+(lg-1)*a*0.5),
- (k*a+a*0.5),xvelocity,yvelocity,zvelocity,sh_force,sh_force2,id,this);
- ijk = ijk + 1;
- }
+ for (int lg=1; lg<=2; lg++) {
+ for (int i=0; i<mm; i++) {
+ for (int j=0; j<mm; j++) {
+ for (int k=0; k<mm; k++) {
+ one[ijk] = new particle((i*a+(2-lg)*a*0.5),(j*a+(lg-1)*a*0.5),
+ (k*a+a*0.5),xvelocity,yvelocity,zvelocity,sh_force,sh_force2,id,one);
+ ijk = ijk + 1;
}
}
}
+ }
- /* Initialise velocities */
+ /* Initialise velocities */
- iseed = 0;
- v1 = 0.0;
- v2 = 0.0;
- randnum = new random(iseed,v1,v2);
+ iseed = 0;
+ double v1 = 0.0;
+ double v2 = 0.0;
+ random randnum = new random(iseed,v1,v2);
- for (i=0; i<tmpmdsize; i+=2) {
- r = randnum.seed();
- one[i].xvelocity = r*randnum.v1;
- one[i+1].xvelocity = r*randnum.v2;
- }
- for (i=0; i<tmpmdsize; i+=2) {
- r = randnum.seed();
- one[i].yvelocity = r*randnum.v1;
- one[i+1].yvelocity = r*randnum.v2;
- }
-
- for (i=0; i<tmpmdsize; i+=2) {
- r = randnum.seed();
- one[i].zvelocity = r*randnum.v1;
- one[i+1].zvelocity = r*randnum.v2;
- }
+ for (int i=0; i<mdsize; i+=2) {
+ r = randnum.seed();
+ one[i].xvelocity = r*randnum.v1;
+ one[i+1].xvelocity = r*randnum.v2;
+ }
+ for (int i=0; i<mdsize; i+=2) {
+ r = randnum.seed();
+ one[i].yvelocity = r*randnum.v1;
+ one[i+1].yvelocity = r*randnum.v2;
+ }
- /* velocity scaling */
+ for (int i=0; i<mdsize; i+=2) {
+ r = randnum.seed();
+ one[i].zvelocity = r*randnum.v1;
+ one[i+1].zvelocity = r*randnum.v2;
+ }
- ekin = 0.0;
- sp = 0.0;
- for(i=0;i<tmpmdsize;i++) {
- sp = sp + one[i].xvelocity;
- }
- sp = sp / tmpmdsize;
+ /* velocity scaling */
- for(i=0;i<tmpmdsize;i++) {
- one[i].xvelocity = one[i].xvelocity - sp;
- ekin = ekin + one[i].xvelocity*one[i].xvelocity;
- }
+ ekin = 0.0;
+ sp = 0.0;
- sp = 0.0;
- for(i=0;i<tmpmdsize;i++) {
- sp = sp + one[i].yvelocity;
- }
- sp = sp / tmpmdsize;
+ for(int i=0;i<mdsize;i++) {
+ sp = sp + one[i].xvelocity;
+ }
+ sp = sp / mdsize;
- for(i=0;i<tmpmdsize;i++) {
- one[i].yvelocity = one[i].yvelocity - sp;
- ekin = ekin + one[i].yvelocity*one[i].yvelocity;
- }
+ for(int i=0;i<mdsize;i++) {
+ one[i].xvelocity = one[i].xvelocity - sp;
+ ekin = ekin + one[i].xvelocity*one[i].xvelocity;
+ }
+ sp = 0.0;
+ for(int i=0;i<mdsize;i++) {
+ sp = sp + one[i].yvelocity;
+ }
+ sp = sp / mdsize;
- sp = 0.0;
- for(i=0;i<tmpmdsize;i++) {
- sp = sp + one[i].zvelocity;
- }
- sp = sp / tmpmdsize;
+ for(int i=0;i<mdsize;i++) {
+ one[i].yvelocity = one[i].yvelocity - sp;
+ ekin = ekin + one[i].yvelocity*one[i].yvelocity;
+ }
- for(i=0;i<tmpmdsize;i++) {
- one[i].zvelocity = one[i].zvelocity - sp;
- ekin = ekin + one[i].zvelocity*one[i].zvelocity;
- }
- ts = tscale * ekin;
- sc = h * Math.sqrt(tref/ts);
+ sp = 0.0;
+ for(int i=0;i<mdsize;i++) {
+ sp = sp + one[i].zvelocity;
+ }
+ sp = sp / mdsize;
+ for(int i=0;i<mdsize;i++) {
+ one[i].zvelocity = one[i].zvelocity - sp;
+ ekin = ekin + one[i].zvelocity*one[i].zvelocity;
+ }
- for(i=0;i<tmpmdsize;i++) {
+ ts = tscale * ekin;
+ sc = h * Math.sqrt(tref/ts);
- one[i].xvelocity = one[i].xvelocity * sc;
- one[i].yvelocity = one[i].yvelocity * sc;
- one[i].zvelocity = one[i].zvelocity * sc;
- }
+ for(int i=0;i<mdsize;i++) {
- /* Synchronise threads and start timer before MD simulation */
+ one[i].xvelocity = one[i].xvelocity * sc;
+ one[i].yvelocity = one[i].yvelocity * sc;
+ one[i].zvelocity = one[i].zvelocity * sc;
- Barrier.enterBarrier(tmpbr);
- //System.clearPrefetchCache();
- //int myid;
- //atomic {
- // myid = id;
- //}
- //TournamentBarrier.enterBarrier(myid, tmpbr);
- //if (id == 0) JGFInstrumentor.startTimer("Section3:MolDyn:Run", instr.timers);
- //Barrier.enterBarrier(tmpbr);
+ }
- /* MD simulation */
+ }
- for (int move=0;move<movemx;move++) {
- /* move the particles and update velocities */
+ public void doinit(int mdsize) {
+ for(int j=0;j<3;j++) {
+ double[] sh=sh_force[j];
+ for (int i=0;i<mdsize;i++) {
+ sh[i] = 0.0;
+ }
+ }
+ }
- for (i=0;i<tmpmdsize;i++) {
- one[i].domove(side,i);
- }
- /* Barrier */
- //System.printString("Barrier #2\n");
- Barrier.enterBarrier(tmpbr);
- //System.clearPrefetchCache();
- //TournamentBarrier.enterBarrier(myid, tmpbr);
-
- if(id==0) {
- for(j=0;j<3;j++) {
- for (i=0;i<tmpmdsize;i++) {
- sh_force[j][i] = 0.0;
- }
- }
+ public void doinit2(int mdsize) {
+ for(int k=0;k<3;k++) {
+ double[] sh=sh_force[k];
+ double [][] sha=sh_force2[k];
+ for(int j=0;j<nthreads;j++) {
+ double[] sha2=sha[j];
+ for(int i=0;i<mdsize;i++) {
+ sh[i] += sha2[i];
}
+ }
+ }
- mymd.epot[id] = 0.0;
- mymd.vir[id] = 0.0;
- mymd.interacts[id] = 0;
-
-
- /* Barrier */
- //System.printString("Barrier #3\n");
- Barrier.enterBarrier(tmpbr);
- //System.clearPrefetchCache();
- //TournamentBarrier.enterBarrier(myid, tmpbr);
+ for(int k=0;k<3;k++) {
+ double [][] sh1=sh_force2[k];
+ for(int j=0;j<nthreads;j++) {
+ double[] sh2=sh1[j];
+ for(int i=0;i<mdsize;i++) {
- /* compute forces */
- for (i=0+id;i<tmpmdsize;i+=nthreads) {
- one[i].force(side,rcoff,tmpmdsize,i,xx,yy,zz,mymd);
+ sh2[i] = 0.0;
}
+ }
+ }
- /* Barrier */
- //System.printString("Barrier #4\n");
- Barrier.enterBarrier(tmpbr);
- //System.clearPrefetchCache();
- //TournamentBarrier.enterBarrier(myid, tmpbr);
+ for(int j=1;j<nthreads;j++) {
+ mymd.epot[0].d += mymd.epot[j].d;
+ mymd.vir[0].d += mymd.vir[j].d;
+ }
+ for(int j=1;j<nthreads;j++) {
+ mymd.epot[j].d = mymd.epot[0].d;
+ mymd.vir[j].d = mymd.vir[0].d;
+ }
+ for(int j=0;j<nthreads;j++) {
+ mymd.interactions += mymd.interacts[j].i;
+ }
- /* update force arrays */
- if(id == 0) {
- for(int k=0;k<3;k++) {
- for(i=0;i<tmpmdsize;i++) {
- for(j=0;j<nthreads;j++) {
- sh_force[k][i] += sh_force2[k][j][i];
- }
- }
- }
- }
+ for (int j=0;j<3;j++) {
+ double sh[]=sh_force[j];
+ for (int i=0;i<mdsize;i++) {
+ sh[i] = sh[i] * hsq2;
+ }
+ }
+ }
- if(id == 0) {
- for(int k=0;k<3;k++) {
- for(i=0;i<tmpmdsize;i++) {
- for(j=0;j<nthreads;j++) {
- sh_force2[k][j][i] = 0.0;
- }
- }
- }
- }
+ public void run() {
+ /* Parameter determination */
- if(id==0) {
- for(j=1;j<nthreads;j++) {
- mymd.epot[0] += mymd.epot[j];
- mymd.vir[0] += mymd.vir[j];
- }
- for(j=1;j<nthreads;j++) {
- mymd.epot[j] = mymd.epot[0];
- mymd.vir[j] = mymd.vir[0];
- }
- for(j=0;j<nthreads;j++) {
- mymd.interactions += mymd.interacts[j];
- }
- }
+ int mdsize;
+ double tmpden;
+ int movemx=50;
+ //Barrier barr=new Barrier("128.195.175.84");
+ particle[] one;
+ int id;
+ id=this.id;
+ mdsize = mymd.PARTSIZE;
+ one=new particle[mdsize];
+ l = mymd.LENGTH;
+ tmpden = den;
+ side = Math.pow((mdsize/tmpden),0.3333333);
+ rcoff = mm/4.0;
+
+ a = side/mm;
+ sideh = side*0.5;
+ hsq = h*h;
+ hsq2 = hsq*0.5;
+ npartm = mdsize - 1;
+ rcoffs = rcoff * rcoff;
+ tscale = 16.0 / (1.0 * mdsize - 1.0);
+ vaver = 1.13 * Math.sqrt(tref / 24.0);
+ vaverh = vaver * h;
+
+ /* Particle Generation */
+
+ xvelocity = 0.0;
+ yvelocity = 0.0;
+ zvelocity = 0.0;
+ ijk = 0;
+ init(one, mdsize);
- /* Barrier */
- //System.printString("Barrier #5\n");
- Barrier.enterBarrier(tmpbr);
- //System.clearPrefetchCache();
- //TournamentBarrier.enterBarrier(myid, tmpbr);
-
- if(id == 0) {
- for (j=0;j<3;j++) {
- for (i=0;i<tmpmdsize;i++) {
- sh_force[j][i] = sh_force[j][i] * hsq2;
- }
- }
- }
- sum = 0.0;
+ /* Synchronise threads and start timer before MD simulation */
+ //Barrier.enterBarrier(barr);
- /* Barrier */
- //System.printString("Barrier #6\n");
- Barrier.enterBarrier(tmpbr);
- //System.clearPrefetchCache();
- //TournamentBarrier.enterBarrier(myid, tmpbr);
+ /* MD simulation */
- /*scale forces, update velocities */
+ for (int move=0;move<movemx;move++) {
+ /* move the particles and update velocities */
+ for (int i=0;i<mdsize;i++) {
+ one[i].domove(side,i);
+ }
- for (i=0;i<tmpmdsize;i++) {
- sum = sum + one[i].mkekin(hsq2,i);
- }
+ /* Barrier */
+ //Barrier.enterBarrier(barr);
- ekin = sum/hsq;
+ if(id==0) {
+ doinit(mdsize);
+ }
- vel = 0.0;
- count = 0.0;
+ mymd.epot[id].d = 0.0;
+ mymd.vir[id].d = 0.0;
+ mymd.interacts[id].i = 0;
- /* average velocity */
- for (i=0;i<tmpmdsize;i++) {
- velt = one[i].velavg(vaverh,h);
- if(velt > vaverh) { count = count + 1.0; }
- vel = vel + velt;
- }
+ /* Barrier */
+ //Barrier.enterBarrier(barr);
- vel = vel / h;
+ /* compute forces */
- /* temperature scale if required */
+ for (int i=0+id;i<mdsize;i+=nthreads) {
+ one[i].force(side,rcoff,mdsize,i,xx,yy,zz,mymd);
+ }
- if((move < istop) && (((move+1) % irep) == 0)) {
- sc = Math.sqrt(tref / (tscale*ekin));
- for (i=0;i<tmpmdsize;i++) {
- one[i].dscal(sc,1);
- }
- ekin = tref / tscale;
- }
+ /* Barrier */
+ //Barrier.enterBarrier(barr);
- /* sum to get full potential energy and virial */
+ /* update force arrays */
+ if(id == 0) {
+ doinit2(mdsize);
+ }
- if(((move+1) % iprint) == 0) {
- mymd.ek[id] = 24.0*ekin;
- mymd.epot[id] = 4.0*mymd.epot[id];
- etot = mymd.ek[id] + mymd.epot[id];
- temp = tscale * ekin;
- pres = tmpden * 16.0 * (ekin - mymd.vir[id]) / tmpmdsize;
- vel = vel / tmpmdsize;
- rp = (count / tmpmdsize) * 100.0;
- }
- //System.printString("Barrier #7\n");
- Barrier.enterBarrier(tmpbr);
- //System.clearPrefetchCache();
- //TournamentBarrier.enterBarrier(myid, tmpbr);
- }
-
- //System.printString("Barrier #8\n");
- Barrier.enterBarrier(tmpbr);
- //System.clearPrefetchCache();
- //TournamentBarrier.enterBarrier(myid, tmpbr);
- //if (id == 0) JGFInstrumentor.stopTimer("Section3:MolDyn:Run", instr.timers);
- //System.printString("End run method\n");
- }
+ /* Barrier */
+ //Barrier.enterBarrier(barr);
-}
+ /*scale forces, update velocities */
+ sum = 0.0;
+ for (int i=0;i<mdsize;i++) {
+ sum = sum + one[i].mkekin(hsq2,i);
+ }
+ ekin = sum/hsq;
+ vel = 0.0;
+ count = 0.0;
+ /* average velocity */
-class particle {
+ for (int i=0;i<mdsize;i++) {
+ velt = one[i].velavg(vaverh,h);
+ if(velt > vaverh) { count = count + 1.0; }
+ vel = vel + velt;
+ }
- public double xcoord, ycoord, zcoord;
- public double xvelocity,yvelocity,zvelocity;
- int part_id;
- int id;
- double [][] sh_force;
- double [][][] sh_force2;
- mdRunner runner;
-
- public particle(double xcoord, double ycoord, double zcoord, double xvelocity,
- double yvelocity,double zvelocity,double [][] sh_force,
- double [][][] sh_force2,int id,mdRunner runner) {
-
- this.xcoord = xcoord;
- this.ycoord = ycoord;
- this.zcoord = zcoord;
- this.xvelocity = xvelocity;
- this.yvelocity = yvelocity;
- this.zvelocity = zvelocity;
- this.sh_force = sh_force;
- this.sh_force2 = sh_force2;
- this.id=id;
- this.runner=runner;
- }
+ vel = vel / h;
- public void domove(double side,int part_id) {
+ /* temperature scale if required */
- xcoord = xcoord + xvelocity + sh_force[0][part_id];
- ycoord = ycoord + yvelocity + sh_force[1][part_id];
- zcoord = zcoord + zvelocity + sh_force[2][part_id];
+ if((move < istop) && (((move+1) % irep) == 0)) {
+ sc = Math.sqrt(tref / (tscale*ekin));
+ for (int i=0;i<mdsize;i++) {
+ one[i].dscal(sc,1);
+ }
+ ekin = tref / tscale;
+ }
- if(xcoord < 0) { xcoord = xcoord + side; }
- if(xcoord > side) { xcoord = xcoord - side; }
- if(ycoord < 0) { ycoord = ycoord + side; }
- if(ycoord > side) { ycoord = ycoord - side; }
- if(zcoord < 0) { zcoord = zcoord + side; }
- if(zcoord > side) { zcoord = zcoord - side; }
+ /* sum to get full potential energy and virial */
- xvelocity = xvelocity + sh_force[0][part_id];
- yvelocity = yvelocity + sh_force[1][part_id];
- zvelocity = zvelocity + sh_force[2][part_id];
+ if(((move+1) % iprint) == 0) {
+ mymd.ek[id].d = 24.0*ekin;
+ mymd.epot[id].d = 4.0*mymd.epot[id].d;
+ etot = mymd.ek[id].d + mymd.epot[id].d;
+ temp = tscale * ekin;
+ pres = tmpden * 16.0 * (ekin - mymd.vir[id].d) / mdsize;
+ vel = vel / mdsize;
+ rp = (count / mdsize) * 100.0;
+ }
+ //Barrier.enterBarrier(barr);
+ //if (id == 0) JGFInstrumentor.stopTimer("Section3:MolDyn:Run", instr.timers);
+ }
}
+}
- public void force(double side, double rcoff,int mdsize,int x, double xx, double yy, double zz, JGFMolDynBench mymd) {
+ class particle {
+
+ public double xcoord, ycoord, zcoord;
+ public double xvelocity,yvelocity,zvelocity;
+ int part_id;
+ int id;
+ double [][] sh_force;
+ double [][][] sh_force2;
+ particle[] one;
+
+ public particle(double xcoord, double ycoord, double zcoord, double xvelocity,
+ double yvelocity,double zvelocity, double [][] sh_force,
+ double [][][] sh_force2,int id, particle[] one) {
+
+ this.xcoord = xcoord;
+ this.ycoord = ycoord;
+ this.zcoord = zcoord;
+ this.xvelocity = xvelocity;
+ this.yvelocity = yvelocity;
+ this.zvelocity = zvelocity;
+ this.sh_force = sh_force;
+ this.sh_force2 = sh_force2;
+ this.id=id;
+ this.one=one;
+ }
- double sideh;
- double rcoffs;
+ public void domove(double side,int part_id) {
- double fxi,fyi,fzi;
- double rd,rrd,rrd2,rrd3,rrd4,rrd6,rrd7,r148;
- double forcex,forcey,forcez;
+ xcoord = xcoord + xvelocity + sh_force[0][part_id];
+ ycoord = ycoord + yvelocity + sh_force[1][part_id];
+ zcoord = zcoord + zvelocity + sh_force[2][part_id];
- sideh = 0.5*side;
- rcoffs = rcoff*rcoff;
+ if(xcoord < 0) { xcoord = xcoord + side; }
+ if(xcoord > side) { xcoord = xcoord - side; }
+ if(ycoord < 0) { ycoord = ycoord + side; }
+ if(ycoord > side) { ycoord = ycoord - side; }
+ if(zcoord < 0) { zcoord = zcoord + side; }
+ if(zcoord > side) { zcoord = zcoord - side; }
- fxi = 0.0;
- fyi = 0.0;
- fzi = 0.0;
+ xvelocity = xvelocity + sh_force[0][part_id];
+ yvelocity = yvelocity + sh_force[1][part_id];
+ zvelocity = zvelocity + sh_force[2][part_id];
- for (int i=x+1;i<mdsize;i++) {
- xx = this.xcoord - runner.one[i].xcoord;
- yy = this.ycoord - runner.one[i].ycoord;
- zz = this.zcoord - runner.one[i].zcoord;
+ }
- if(xx < (-sideh)) { xx = xx + side; }
- if(xx > (sideh)) { xx = xx - side; }
- if(yy < (-sideh)) { yy = yy + side; }
- if(yy > (sideh)) { yy = yy - side; }
- if(zz < (-sideh)) { zz = zz + side; }
- if(zz > (sideh)) { zz = zz - side; }
+ public void force(double side, double rcoff,int mdsize,int x, double xx, double yy, double zz, JGFMolDynBench mymd) {
+ double sideh;
+ double rcoffs;
- rd = xx*xx + yy*yy + zz*zz;
+ double fxi,fyi,fzi;
+ double rd,rrd,rrd2,rrd3,rrd4,rrd6,rrd7,r148;
+ double forcex,forcey,forcez;
+ int id=this.id;
+ sideh = 0.5*side;
+ rcoffs = rcoff*rcoff;
- if(rd <= rcoffs) {
- rrd = 1.0/rd;
- rrd2 = rrd*rrd;
- rrd3 = rrd2*rrd;
- rrd4 = rrd2*rrd2;
- rrd6 = rrd2*rrd4;
- rrd7 = rrd6*rrd;
- mymd.epot[id] = mymd.epot[id] + (rrd6 - rrd3);
- r148 = rrd7 - 0.5*rrd4;
- mymd.vir[id] = mymd.vir[id] - rd*r148;
- forcex = xx * r148;
- fxi = fxi + forcex;
+ fxi = 0.0;
+ fyi = 0.0;
+ fzi = 0.0;
- sh_force2[0][id][i] = sh_force2[0][id][i] - forcex;
+ for (int i=x+1;i<mdsize;i++) {
+ xx = this.xcoord - one[i].xcoord;
+ yy = this.ycoord - one[i].ycoord;
+ zz = this.zcoord - one[i].zcoord;
- forcey = yy * r148;
- fyi = fyi + forcey;
+ if(xx < (-sideh)) { xx = xx + side; }
+ if(xx > (sideh)) { xx = xx - side; }
+ if(yy < (-sideh)) { yy = yy + side; }
+ if(yy > (sideh)) { yy = yy - side; }
+ if(zz < (-sideh)) { zz = zz + side; }
+ if(zz > (sideh)) { zz = zz - side; }
- sh_force2[1][id][i] = sh_force2[1][id][i] - forcey;
- forcez = zz * r148;
- fzi = fzi + forcez;
+ rd = xx*xx + yy*yy + zz*zz;
- sh_force2[2][id][i] = sh_force2[2][id][i] - forcez;
+ if(rd <= rcoffs) {
+ rrd = 1.0/rd;
+ rrd2 = rrd*rrd;
+ rrd3 = rrd2*rrd;
+ rrd4 = rrd2*rrd2;
+ rrd6 = rrd2*rrd4;
+ rrd7 = rrd6*rrd;
+ mymd.epot[id].d += (rrd6 - rrd3);
+ r148 = rrd7 - 0.5*rrd4;
+ mymd.vir[id].d += - rd*r148;
+ forcex = xx * r148;
+ fxi = fxi + forcex;
- mymd.interacts[id]++;
- }
+ sh_force2[0][id][i] = sh_force2[0][id][i] - forcex;
- }
+ forcey = yy * r148;
+ fyi = fyi + forcey;
- sh_force2[0][id][x] = sh_force2[0][id][x] + fxi;
- sh_force2[1][id][x] = sh_force2[1][id][x] + fyi;
- sh_force2[2][id][x] = sh_force2[2][id][x] + fzi;
+ sh_force2[1][id][i] = sh_force2[1][id][i] - forcey;
- }
+ forcez = zz * r148;
+ fzi = fzi + forcez;
- public double mkekin(double hsq2,int part_id) {
+ sh_force2[2][id][i] = sh_force2[2][id][i] - forcez;
- double sumt = 0.0;
+ mymd.interacts[id].i++;
+ }
- xvelocity = xvelocity + sh_force[0][part_id];
- yvelocity = yvelocity + sh_force[1][part_id];
- zvelocity = zvelocity + sh_force[2][part_id];
+ }
- sumt = (xvelocity*xvelocity)+(yvelocity*yvelocity)+(zvelocity*zvelocity);
- return sumt;
- }
+ sh_force2[0][id][x] = sh_force2[0][id][x] + fxi;
+ sh_force2[1][id][x] = sh_force2[1][id][x] + fyi;
+ sh_force2[2][id][x] = sh_force2[2][id][x] + fzi;
- public double velavg(double vaverh,double h) {
+ }
- double velt;
- double sq;
+ public double mkekin(double hsq2,int part_id) {
- sq = Math.sqrt(xvelocity*xvelocity + yvelocity*yvelocity +
- zvelocity*zvelocity);
+ double sumt = 0.0;
- velt = sq;
- return velt;
- }
+ xvelocity = xvelocity + sh_force[0][part_id];
+ yvelocity = yvelocity + sh_force[1][part_id];
+ zvelocity = zvelocity + sh_force[2][part_id];
- public void dscal(double sc,int incx) {
+ sumt = (xvelocity*xvelocity)+(yvelocity*yvelocity)+(zvelocity*zvelocity);
+ return sumt;
+ }
- xvelocity = xvelocity * sc;
- yvelocity = yvelocity * sc;
- zvelocity = zvelocity * sc;
+ public double velavg(double vaverh,double h) {
+ double velt;
+ double sq;
+ sq = Math.sqrt(xvelocity*xvelocity + yvelocity*yvelocity +
+ zvelocity*zvelocity);
- }
+ velt = sq;
+ return velt;
+ }
-}
+ public void dscal(double sc,int incx) {
+ xvelocity = xvelocity * sc;
+ yvelocity = yvelocity * sc;
+ zvelocity = zvelocity * sc;
+ }
+ }
-class random {
+ class random {
- public int iseed;
- public double v1,v2;
+ public int iseed;
+ public double v1,v2;
- public random(int iseed,double v1,double v2) {
- this.iseed = iseed;
- this.v1 = v1;
- this.v2 = v2;
- }
+ public random(int iseed,double v1,double v2) {
+ this.iseed = iseed;
+ this.v1 = v1;
+ this.v2 = v2;
+ }
- public double update() {
+ public double update() {
- double rand;
- double scale= 4.656612875e-10;
+ double rand;
+ double scale= 4.656612875e-10;
- int is1,is2,iss2;
- int imult=16807;
- int imod = 2147483647;
+ int is1,is2,iss2;
+ int imult=16807;
+ int imod = 2147483647;
- if (iseed<=0) { iseed = 1; }
+ if (iseed<=0) { iseed = 1; }
- is2 = iseed % 32768;
- is1 = (iseed-is2)/32768;
- iss2 = is2 * imult;
- is2 = iss2 % 32768;
- is1 = (is1*imult+(iss2-is2)/32768) % (65536);
+ is2 = iseed % 32768;
+ is1 = (iseed-is2)/32768;
+ iss2 = is2 * imult;
+ is2 = iss2 % 32768;
+ is1 = (is1*imult+(iss2-is2)/32768) % (65536);
- iseed = (is1*32768+is2) % imod;
+ iseed = (is1*32768+is2) % imod;
- rand = scale * iseed;
+ rand = scale * iseed;
- return rand;
+ return rand;
- }
+ }
- public double seed() {
+ public double seed() {
- double s,u1,u2,r;
- s = 1.0;
- do {
- u1 = update();
- u2 = update();
+ double s,u1,u2,r;
+ s = 1.0;
+ do {
+ u1 = update();
+ u2 = update();
- v1 = 2.0 * u1 - 1.0;
- v2 = 2.0 * u2 - 1.0;
- s = v1*v1 + v2*v2;
+ v1 = 2.0 * u1 - 1.0;
+ v2 = 2.0 * u2 - 1.0;
+ s = v1*v1 + v2*v2;
- } while (s >= 1.0);
+ } while (s >= 1.0);
- r = Math.sqrt(-2.0*Math.log(s)/s);
+ r = Math.sqrt(-2.0*Math.log(s)/s);
- return r;
+ return r;
+ }
}
-}
projects / IRC.git / commitdiff