summaryrefslogtreecommitdiff
path: root/src/spheres_poly
diff options
context:
space:
mode:
Diffstat (limited to 'src/spheres_poly')
-rw-r--r--src/spheres_poly/box.cpp1157
-rw-r--r--src/spheres_poly/box.h169
-rw-r--r--src/spheres_poly/event.cpp65
-rw-r--r--src/spheres_poly/event.h58
-rw-r--r--src/spheres_poly/grid_field.h148
-rw-r--r--src/spheres_poly/heap.cpp186
-rw-r--r--src/spheres_poly/heap.h42
-rw-r--r--src/spheres_poly/neighbor.cpp40
-rw-r--r--src/spheres_poly/sphere.cpp66
-rw-r--r--src/spheres_poly/sphere.h44
-rw-r--r--src/spheres_poly/spheres.cpp64
-rw-r--r--src/spheres_poly/spheres.h4
-rw-r--r--src/spheres_poly/vector.h471
13 files changed, 2514 insertions, 0 deletions
diff --git a/src/spheres_poly/box.cpp b/src/spheres_poly/box.cpp
new file mode 100644
index 0000000..716090d
--- /dev/null
+++ b/src/spheres_poly/box.cpp
@@ -0,0 +1,1157 @@
+/*
+ Updated July 24, 2009 to include hardwall boundary condition option, as
+ well as polydisperse spheres.
+
+ Packing of hard spheres via molecular dynamics
+ Developed by Monica Skoge, 2006, Princeton University
+ Contact: Monica Skoge (mskoge.princeton.edu) with questions
+ This code may be used, modified and distributed freely.
+ Please cite:
+
+ "Packing Hyperspheres in High-Dimensional Euclidean Spaces"
+ M. Skoge, A. Donev, F. H. Stillinger and S. Torquato, 2006
+
+ if you use these codes.
+*/
+
+#include "box.h"
+#include <iostream>
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+#include <iomanip>
+
+//==============================================================
+//==============================================================
+// Class Box: Fills box with hardspheres to given packing fraction
+// and evolves spheres using molecular dynamics!
+//==============================================================
+//==============================================================
+
+
+//==============================================================
+// Constructor
+//==============================================================
+box::box(int N_i, double r_i, double growthrate_i, double maxpf_i,
+ double bidispersityratio_i, double bidispersityfraction_i,
+ double massratio_i, int hardwallBC_i):
+ r(r_i),
+ N(N_i),
+ growthrate(growthrate_i),
+ h(N_i+1),
+ maxpf(maxpf_i),
+ bidispersityratio(bidispersityratio_i),
+ bidispersityfraction(bidispersityfraction_i),
+ massratio(massratio_i),
+ hardwallBC(hardwallBC_i)
+{
+ ngrids = Optimalngrids2(r);
+ cells.set_size(ngrids);
+
+ s = new sphere[N];
+ binlist = new int[N];
+ x = new vector<DIM>[N];
+ h.s = s;
+
+ gtime = 0.;
+ rtime = 0.;
+ ncollisions = 0;
+ ntransfers = 0;
+ nchecks = 0;
+ neventstot = 0;
+ ncycles = 0;
+ xmomentum = 0.;
+ pressure = 0.;
+ collisionrate = 0.;
+
+ cells.set_size(ngrids);
+ cells.initialize(-1); // initialize cells to -1
+ srandom(::time(0)); // initialize the random number generator
+ for (int i=0; i<N; i++) // initialize binlist to -1
+ binlist[i] = -1;
+
+ time(&start);
+}
+
+
+//==============================================================
+// Destructor
+//==============================================================
+box::~box()
+{
+ delete[] s;
+ delete[] binlist;
+ delete[] x;
+}
+
+
+//==============================================================
+// ReadFile
+//==============================================================
+void box::ReadPositions(const char* filename)
+{
+ // open file to read in arrays
+ std::ifstream infile(filename);
+
+ infile.ignore(256, '\n'); // ignore the dim line
+ infile.ignore(256, '\n'); // ignore the #sphere 1 line
+ infile.ignore(256, '\n'); // ignore the #sphere line
+ infile.ignore(256, '\n'); // ignore the diameter line
+ infile.ignore(1000, '\n'); // ignore the 100 010 001 line
+ infile.ignore(256, '\n'); // ignore the T T T line
+
+ for (int i=0; i<N; i++)
+ {
+ infile >> s[i].r; // read in radius
+ infile >> s[i].gr; // read in growth rate
+ infile >> s[i].m; // read in mass
+ for (int k=0; k<DIM; k++)
+ infile >> s[i].x[k]; // read in position
+ }
+ infile.close();
+}
+
+
+//==============================================================
+// Recreates all N spheres at random positions
+//==============================================================
+void box::RecreateSpheres(const char* filename, double temp)
+{
+ ReadPositions(filename); // reads in positions of spheres
+ VelocityGiver(temp); // gives spheres initial velocities
+ AssignCells(); // assigns spheres to cells
+ SetInitialEvents();
+}
+
+
+//==============================================================
+// Creates all N spheres at random positions
+//==============================================================
+void box::CreateSpheres(double temp)
+{
+ int Ncurrent = 0;
+ for(int i=0; i<N; i++)
+ {
+ CreateSphere(Ncurrent);
+ Ncurrent++;
+ }
+ if (Ncurrent != N)
+ std::cout << "problem! only made " << Ncurrent << " out of " << N << " desired spheres" << std::endl;
+
+ VelocityGiver(temp);
+ SetInitialEvents();
+}
+
+
+//==============================================================
+// Creates a sphere of random radius r at a random unoccupied position
+//==============================================================
+void box::CreateSphere(int Ncurrent)
+{
+ int keeper; // boolean variable: 1 means ok, 0 means sphere already there
+ int counter = 0; // counts how many times sphere already exists
+ vector<DIM> xrand; // random new position vector
+ double d = 0.;
+ double radius;
+ double growth_rate;
+ double mass;
+ int species;
+
+ if (Ncurrent < bidispersityfraction*N)
+ {
+ radius = r;
+ growth_rate = growthrate;
+ mass = 1.;
+ species = 1;
+ }
+ else
+ {
+ radius = r*bidispersityratio;
+ growth_rate = growthrate*bidispersityratio;
+ mass = massratio;
+ species = 2;
+ }
+
+ while (counter<1000)
+ {
+ keeper = 1;
+
+ for(int k=0; k<DIM; k++)
+ xrand[k] = ((double)random()/(double)RAND_MAX)*SIZE;
+
+ for (int i=0; i<Ncurrent; i++) // check if overlapping other spheres
+ {
+ d=0.;
+ if (hardwallBC)
+ {
+ for (int k=0; k<DIM; k++)
+ {
+ d += (xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[k]);
+ }
+ }
+ else
+ {
+ for (int k=0; k<DIM; k++)
+ {
+ if ((xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[k]) > SIZE*SIZE/4.)
+ {
+ if (xrand[k] > SIZE/2.) // look at right image
+ d += (xrand[k] - (s[i].x[k]+SIZE))*
+ (xrand[k] - (s[i].x[k]+SIZE));
+ else // look at left image
+ d += (xrand[k] - (s[i].x[k]-SIZE))*
+ (xrand[k] - (s[i].x[k]-SIZE));
+ }
+ else
+ d += (xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[k]);
+ }
+ }
+ if (d <= (radius + s[i].r)*(radius + s[i].r)) // overlapping!
+ {
+ keeper = 0;
+ counter++;
+ break;
+ }
+ }
+
+ if (hardwallBC)
+ {
+ for (int k=0; k<DIM; k++) // check if overlapping wall
+ {
+ if ((xrand[k] <= radius) || (SIZE - xrand[k] <= radius)) // touching wall
+ {
+ keeper = 0;
+ counter++;
+ break;
+ }
+ }
+ }
+
+ if (keeper == 1)
+ break;
+
+ if (counter >= 1000)
+ {
+ std::cout << "counter >= 1000" << std::endl;
+ exit(-1);
+ }
+ }
+
+ // now convert xrand into index vector for cells
+ vector<DIM,int> cell;
+ cell = vector<DIM>::integer(xrand*((double)(ngrids))/SIZE);
+
+ s[Ncurrent] = sphere(Ncurrent, xrand, cell, gtime, radius, growth_rate,
+ mass, species);
+
+ //first check to see if entry at cell
+ if (cells.get(cell) == -1) //if yes, add Ncurrent to cells gridfield
+ cells.get(cell) = Ncurrent;
+
+ else // if no, add i to right place in binlist
+ {
+ int iterater = cells.get(cell); // now iterate through to end and add Ncurrent
+ int pointer = iterater;
+ while (iterater != -1)
+ {
+ pointer = iterater;
+ iterater = binlist[iterater];
+ }
+ binlist[pointer] = Ncurrent;
+ }
+
+}
+
+
+//==============================================================
+// Assign cells to spheres read in from existing configuration
+//==============================================================
+void box::AssignCells()
+{
+ for (int i=0; i<N; i++)
+ {
+ // now convert x into index vector for cells
+ vector<DIM,int> cell;
+ cell = vector<DIM>::integer(s[i].x*((double)(ngrids))/SIZE);
+ s[i].cell = cell;
+
+ //first check to see if entry at cell
+ if (cells.get(cell) == -1) //if yes, add Ncurrent to cells gridfield
+ cells.get(cell) = i;
+
+ else // if no, add i to right place in binlist
+ {
+ int iterater = cells.get(cell); // now iterate through to end and add Ncurrent
+ int pointer = iterater;
+ while (iterater != -1)
+ {
+ pointer = iterater;
+ iterater = binlist[iterater];
+ }
+ binlist[pointer] = i;
+ }
+ }
+}
+
+
+//==============================================================
+// Velocity Giver, assigns initial velocities from Max/Boltz dist.
+//==============================================================
+void box::VelocityGiver(double T)
+{
+ for (int i=0; i<N; i++)
+ {
+ for (int k=0; k<DIM; k++)
+ {
+ if (T==0.)
+ s[i].v[k] = 0.;
+ else
+ s[i].v[k] = Velocity(T);
+ }
+ }
+}
+
+
+//==============================================================
+// Velocity, gives a single velocity from Max/Boltz dist.
+//==============================================================
+double box::Velocity(double T)
+{
+ double rand; // random number between -0.5 and 0.5
+ double sigmasquared = T; // Assumes M = mass of sphere = 1
+ double sigma = sqrt(sigmasquared); // variance of Gaussian
+ double stepsize = 1000.; // stepsize for discretization of integral
+ double vel = 0.0; // velocity
+ double dv=sigma/stepsize;
+ double p=0.0;
+
+ rand = (double)random() / (double)RAND_MAX - 0.5;
+ if(rand < 0)
+ {
+ rand = -rand;
+ dv = -dv;
+ }
+
+ while(fabs(p) < rand) // integrate until the integral equals rand
+ {
+ p += dv * 0.39894228 * exp(-vel*vel/(2.*sigmasquared))/sigma;
+ vel += dv;
+ }
+ return vel;
+}
+
+
+//==============================================================
+// Finds next events for all spheres..do this once at beginning
+//==============================================================
+void box::SetInitialEvents()
+{
+ for (int i=0; i<N; i++) // set all events to checks
+ {
+ event e(gtime, i, INF);
+ s[i].nextevent = e;
+ h.insert(i);
+ }
+}
+
+
+//==============================================================
+// Finds next event for sphere i
+//==============================================================
+event box::FindNextEvent(int i)
+{
+ double outgrowtime;
+ outgrowtime = (.5/ngrids - (s[i].r+s[i].gr*gtime))/s[i].gr + gtime;
+
+ event t = FindNextTransfer(i);
+ event c = FindNextCollision(i);
+
+ if ((outgrowtime < c.time)&&(outgrowtime < t.time)&&(ngrids>1))
+ {
+ event o = event(outgrowtime,i,INF-1);
+ return o;
+ }
+
+ if ((c.time < t.time)&&(c.j == INF)) // next event is check at DBL infinity
+ return c;
+ else if (c.time < t.time) // next event is collision!
+ {
+ CollisionChecker(c);
+ return c;
+ }
+ else // next event is transfer!
+ return t;
+}
+
+
+//==============================================================
+// Checks events of predicted collision partner to keep collisions
+// symmetric
+//==============================================================
+void box::CollisionChecker(event c)
+{
+ int i = c.i;
+ int j = c.j;
+ event cj(c.time,j,i,c.v*(-1));
+
+ // j should have NO event before collision with i!
+ if (!(c.time < s[j].nextevent.time))
+ std::cout << i << " " << j << " error collchecker, s[j].nextevent.time= " << s[j].nextevent.time << " " << s[j].nextevent.j << ", c.time= " << c.time << std::endl;
+
+ int k = s[j].nextevent.j;
+ if ((k < N) && (k!=i)) // j's next event was collision so give k a check
+ s[k].nextevent.j = INF;
+
+ // give collision cj to j
+ s[j].nextevent = cj;
+ h.upheap(h.index[j]);
+}
+
+
+//==============================================================
+// Find next transfer for sphere i
+//==============================================================
+event box::FindNextTransfer(int i)
+{
+ double ttime = dblINF;
+ int wallindex = INF; // -(k+1) left wall, (k+1) right wall
+
+ vector<DIM> xi = s[i].x + s[i].v*(gtime - s[i].lutime);
+ vector<DIM> vi = s[i].v;
+
+ for (int k=0; k<DIM; k++)
+ {
+ double newtime;
+ if (vi[k]==0.)
+ newtime= dblINF;
+ else if (vi[k]>0) // will hit right wall, need to check if last wall
+ {
+ if ((hardwallBC)&&(s[i].cell[k] == ngrids - 1))
+ newtime = ((double)(s[i].cell[k]+1)*SIZE/((double)(ngrids))
+ - (xi[k]+s[i].r+s[i].gr*gtime))/(vi[k]+s[i].gr);
+ else
+ newtime = ((double)(s[i].cell[k]+1)*SIZE/((double)(ngrids))
+ - xi[k])/(vi[k]);
+
+ if (newtime<ttime)
+ {
+ wallindex = k+1;
+ ttime = newtime;
+ }
+ }
+ else if (vi[k]<0) // will hit left wall
+ {
+ if ((hardwallBC)&&(s[i].cell[k] == 0))
+ newtime = ((double)(s[i].cell[k])*SIZE/((double)(ngrids))
+ - (xi[k]-(s[i].r+s[i].gr*gtime)))/(vi[k]-s[i].gr);
+ else
+ newtime = ((double)(s[i].cell[k])*SIZE/((double)(ngrids))
+ - xi[k])/(vi[k]);
+
+ if (newtime<ttime)
+ {
+ wallindex = -(k+1);
+ ttime = newtime;
+ }
+ }
+ }
+
+
+ if (ttime < 0)
+ ttime = 0;
+ // make the event and return it
+ event e = event(ttime+gtime,i,wallindex+DIM+N+1);
+ return e;
+}
+
+
+//==============================================================
+// Check all nearest neighbor cells for collision partners
+//==============================================================
+void box::ForAllNeighbors(int i, vector<DIM,int> vl, vector<DIM,int> vr,
+ neighbor& operation)
+{
+ vector<DIM,int> cell = s[i].cell;
+
+ // now iterate through nearest neighbors
+ vector<DIM, int> offset; // nonnegative neighbor offset
+ vector<DIM, int> pboffset; // nearest image offset
+
+ vector<DIM,int> grid;
+
+ int ii ;
+
+ grid=vl;
+ while(1)
+ {
+ //if (vr[0] > 1)
+ //std::cout << grid << "..." << cell+grid << "\n";
+ for(int k=0; k<DIM; k++)
+ {
+ offset[k]=grid[k]+ngrids; // do this so no negatives
+ if (cell[k]+grid[k]<0) //out of bounds to left
+ pboffset[k] = -1;
+ else if (cell[k]+grid[k]>=ngrids) // out of bounds to right
+ pboffset[k] = 1;
+ else
+ pboffset[k] = 0;
+ }
+ int j = cells.get((cell+offset)%ngrids);
+ while(j!=-1)
+ {
+ operation.Operation(j,pboffset);
+ j = binlist[j];
+ }
+
+ // A. Donev:
+ // This code makes this loop dimension-independent
+ // It is basically a flattened-out loop nest of depth DIM
+ for(ii=0;ii<DIM;ii++)
+ {
+ grid[ii] += 1;
+ if(grid[ii]<=vr[ii]) break;
+ grid[ii]=vl[ii];
+ }
+ if(ii>=DIM) break;
+ }
+}
+
+
+//==============================================================
+// PredictCollision
+//==============================================================
+void box::PredictCollision(int i, int j, vector<DIM, int> pboffset,
+ double& ctime, int& cpartner,
+ vector<DIM, int>& cpartnerpboffset)
+{
+ double ctimej;
+
+ if (i!=j)
+ {
+ ctimej = CalculateCollision(i,j,pboffset.Double())+gtime;
+
+ if (ctimej < gtime)
+ std::cout << "error in find collision ctimej < 0" << std::endl;
+
+ if ((ctimej < ctime)&&(ctimej < s[j].nextevent.time))
+ {
+ ctime = ctimej;
+ cpartner = j;
+ cpartnerpboffset = pboffset;
+ }
+ }
+}
+
+
+//==============================================================
+// Find next collision
+//==============================================================
+event box::FindNextCollision(int i)
+{
+ collision cc(i, this);
+
+ vector<DIM, int> vl, vr;
+
+ for (int k=0; k<DIM; k++) // check all nearest neighbors
+ {
+ vl[k] = -1;
+ vr[k] = 1;
+ }
+
+ ForAllNeighbors(i,vl,vr,cc);
+
+ event e;
+ if (cc.cpartner == i) // found no collisions in neighboring cells
+ {
+ if (cc.ctime != dblINF)
+ std::cout << "ctime != dblINF" << std::endl;
+ e = event(dblINF,i,INF); // give check at double INF
+ }
+ else
+ e = event(cc.ctime,i,cc.cpartner,cc.cpartnerpboffset);
+
+ return e;
+}
+
+
+//==============================================================
+// Calculates collision time between i and image of j using quadratic formula
+//==============================================================
+double box::CalculateCollision(int i, int j, vector<DIM> pboffset)
+{
+ if ((hardwallBC)&&(pboffset.norm_squared() > 1E-12))
+ {
+ return dblINF;
+ }
+ else
+ {
+ // calculate updated position and velocity of i and j
+ vector<DIM> xi = s[i].x + s[i].v*(gtime - s[i].lutime);
+ vector<DIM> vi = s[i].v;
+ vector<DIM> xj = s[j].x + pboffset*SIZE + s[j].v*(gtime - s[j].lutime);
+ vector<DIM> vj = s[j].v;
+
+ double r_now_i = s[i].r + gtime*s[i].gr;
+ double r_now_j = s[j].r + gtime*s[j].gr;
+
+ double A,B,C;
+ A = vector<DIM>::norm_squared(vi - vj) - (s[i].gr+s[j].gr)*(s[i].gr+s[j].gr);
+ B = vector<DIM>::dot(xi - xj, vi - vj) - (r_now_i+r_now_j)*(s[i].gr+s[j].gr);
+ C = vector<DIM>::norm_squared(xi - xj) - (r_now_i+r_now_j)*(r_now_i+r_now_j);
+
+ if (C < -1E-12*(r_now_i+r_now_j))
+ {
+ std::cout << "error, " << i << " and " << j <<
+ " are overlapping at time "<< gtime << std::cout;
+ std::cout << "A, B, C = " << A << " " << " " << B <<
+ " " << " " << C << std::endl;
+ if (CheckSphereDiameters()>0)
+ std::cout << "a sphere has grown greater than unit cell" <<
+ std::endl;
+ else
+ std::cout << "unknown error" << std::cout;
+ exit(-1);
+ }
+
+ return QuadraticFormula(A, B, C);
+ }
+}
+
+
+//==============================================================
+// Quadratic Formula ax^2 + bx + c = 0
+//==============================================================
+ double box::QuadraticFormula(double a, double b, double c)
+{
+ double x = dblINF;
+ double xpos;
+ double xneg;
+ double det = b*b - a*c;
+
+ if (c <= 0.)
+ {
+ if(b < 0.) // spheres already overlapping and approaching
+ {
+ //std::cout << "spheres overlapping and approaching" << std::endl;
+ //std::cout << "# events= " << neventstot << std::endl;
+ x = 0.;
+ }
+ }
+ else if (det > -10.*DBL_EPSILON)
+ {
+ if (det < 0.) // determinant can be very small for double roots
+ det = 0.;
+ if (b < 0.)
+ x = c/(-b + sqrt(det));
+ else if ((a < 0.)&&(b > 0.))
+ x = -(b + sqrt(det))/a;
+ else
+ x = dblINF;
+ }
+ return x;
+}
+
+
+//==============================================================
+// Returns first event
+//==============================================================
+void box::ProcessEvent()
+{
+ neventstot++;
+ // Extract first event from heap
+ int i = h.extractmax();
+ event e = s[i].nextevent; // current event
+ event f; // replacement event
+
+ if ((e.j>=0)&&(e.j<N)) // collision!
+ {
+ ncollisions++;
+ //std::cout << "collision between " << e.i << " and " << e.j << " at time " << e.time << std::endl;
+ Collision(e);
+ f = FindNextEvent(i);
+ s[i].nextevent = f;
+ h.downheap(1);
+ if (f.time < e.time)
+ {
+ std::cout << "error, replacing event with < time" << std::endl;
+ exit(-1);
+ }
+
+ // make sure collision was symmetric and give j a check
+ if ((s[e.j].nextevent.j != i)||(s[e.j].nextevent.time != gtime))
+ {
+ std::cout << "error collisions not symmetric" << std::endl;
+ std::cout << "collision between " << e.i << " and " << e.j << " at time " << e.time << std::endl;
+ std::cout << "but " << e.j << " thinks it has " << s[e.j].nextevent.j<< " " << s[e.j].nextevent.time << std::endl;
+ exit(-1);
+ }
+ else // give j a check
+ s[e.j].nextevent.j = INF;
+ }
+ else if (e.j==INF) // check!
+ {
+ nchecks++;
+ //std::cout << "check for " << e.i << " at time " << e.time << std::endl;
+ f = FindNextEvent(i);
+ s[i].nextevent = f;
+ h.downheap(1);
+ }
+ else if (e.j==INF-1) // sphere outgrowing unit cell, decrease ngrids!
+ {
+ gtime = e.time;
+ Synchronize(false);
+ ngrids = ngrids - 1;
+// std::cout << "need to reduce ngrids to " << ngrids << std::endl;
+ ChangeNgrids(ngrids);
+ h.downheap(1);
+ }
+ else // transfer!
+ {
+ ntransfers++;
+ //std::cout << "transfer for " << e.i << " at time " << e.time << std::endl;
+ Transfer(e);
+ f = FindNextEvent(i);
+ s[i].nextevent = f;
+ h.downheap(1);
+ //r = FindNextEvent(i, e.j-N-DIM-1);
+ //if (f.time <= e.time)
+ if (f.time < e.time)
+ {
+ std::cout << "error after transfer, replacing new event with < time" << " " << std::endl;
+ std::cout << std::setprecision(16) << "e.time= " << e.time << ", f.time= " << f.time << ", f.i= " << f.i << ", f.j= " << f.j << "e.i= " << e.i << ", e.j= " << e.j << std::endl;
+ std::cout << std::setprecision(16) << "difference= " << e.time - f.time << std::endl;
+ exit(-1);
+ }
+ }
+}
+
+
+//==============================================================
+// Processes a collision
+//=============================================================
+void box::Collision(event e)
+{
+ double ctime = e.time;
+ int i = e.i;
+ int j = e.j;
+ vector<DIM,int> v = e.v; // virtual image
+ gtime = ctime;
+
+ // Update positions and cells of i and j to ctime
+ s[i].x += s[i].v*(gtime-s[i].lutime);
+ s[j].x += s[j].v*(gtime-s[j].lutime);
+
+ // Check to see if a diameter apart
+ double r_sum = s[i].r + s[j].r + (s[i].gr+s[j].gr)*gtime;
+ double distance = vector<DIM>::norm_squared(s[i].x - s[j].x- v.Double()*SIZE) - r_sum*r_sum;
+ if (distance*distance > 10.*DBL_EPSILON)
+ std::cout << "overlap " << distance << std::endl;
+
+ s[i].lutime = gtime;
+ s[j].lutime = gtime;
+
+ vector<DIM,double> vipar; // parallel comp. vi
+ vector<DIM,double> vjpar; // parallel comp. vj
+ vector<DIM,double> viperp; // perpendicular comp. vi
+ vector<DIM,double> vjperp; // perpendicular comp. vj
+
+ // make unit vector out of displacement vector
+ vector<DIM,double> dhat;
+ dhat = s[i].x - s[j].x - v.Double()*SIZE; // using image of j!!
+ double dhatmagnitude = sqrt(dhat.norm_squared());
+ dhat /= dhatmagnitude;
+
+ vipar = dhat*vector<DIM>::dot(s[i].v, dhat);
+ vjpar = dhat*vector<DIM>::dot(s[j].v, dhat);
+ viperp = s[i].v - vipar;
+ vjperp = s[j].v - vjpar;
+
+ s[i].v = vjpar + dhat*(s[i].gr+s[j].gr)*2 + viperp;
+ s[j].v = vipar - dhat*(s[i].gr+s[j].gr)*2 + vjperp;
+
+ // momentum exchange
+ double xvelocity; // exchanged velocity
+ xvelocity = vector<DIM>::dot(s[i].v - s[j].v, dhat) - (s[i].gr+s[j].gr);
+ xmomentum += xvelocity*dhatmagnitude*s[i].m*s[j].m*2/(s[i].m+s[j].m);
+}
+
+
+//==============================================================
+// Transfer, takes care of boundary events too
+//=============================================================
+void box::Transfer(event e)
+{
+ gtime = e.time;
+ int i = e.i;
+ int j = e.j;
+ int k=0; // dimension perpendicular to wall it crosses
+
+ // update position and lutime (velocity doesn't change)
+ s[i].x += s[i].v*(gtime-s[i].lutime);
+ s[i].lutime = gtime;
+
+ vector<DIM,int> celli; // new cell for i
+ celli = s[i].cell; // this is not redundant
+
+ // update cell
+ if (j>N+DIM+1) // right wall
+ {
+ k = j-N-DIM-2;
+ celli[k] = s[i].cell[k] + 1;
+
+ if (hardwallBC)
+ {
+ // if in right-most cell, reflect back
+ if (s[i].cell[k] == ngrids - 1)
+ s[i].v[k] = -s[i].v[k] - 4.*s[i].gr;
+ else
+ {
+ if (ngrids>1)
+ UpdateCell(i, celli);
+ }
+ }
+ else
+ {
+ // if in right-most cell, translate x and cell
+ if (s[i].cell[k] == ngrids - 1)
+ {
+ s[i].x[k] -= SIZE;
+ celli[k] -= ngrids;
+ }
+ if (ngrids>1)
+ UpdateCell(i, celli);
+ }
+ }
+ else if (j<N+DIM+1) // left wall
+ {
+ k = -j+N+DIM;
+ celli[k] = s[i].cell[k] - 1;
+
+ if (hardwallBC)
+ {
+ // if in left-most cell, reflect back
+ if (s[i].cell[k] == 0)
+ s[i].v[k] = -s[i].v[k] + 4.*s[i].gr;
+ else
+ {
+ if (ngrids>1)
+ UpdateCell(i, celli);
+ }
+ }
+ else
+ {
+ // if in left-most cell, translate x and cell
+ if (s[i].cell[k] == 0)
+ {
+ s[i].x[k] += SIZE;
+ celli[k] += ngrids;
+ }
+ if (ngrids>1)
+ UpdateCell(i, celli);
+ }
+ }
+ else
+ std::cout << "error in Transfer" << std::endl;
+
+}
+
+
+//==============================================================
+// Updates cell of a sphere to time
+//=============================================================
+void box::UpdateCell(int i, vector<DIM,int>& celli)
+{
+ if (celli == s[i].cell)
+ std::cout << "error in update cell..shouldn't be the same" << std::endl;
+
+ // delete i from cell array at cell
+
+ if (cells.get(s[i].cell) == i)
+ {
+ if (binlist[i] == -1)
+ cells.get(s[i].cell) = -1;
+ else
+ {
+ cells.get(s[i].cell) = binlist[i];
+ binlist[i] = -1;
+ }
+ }
+
+ else if (cells.get(s[i].cell) == -1)
+ {
+ std::cout << "error " << i << " not in claimed cell UpdateCell" << std::endl;
+ OutputCells();
+ }
+
+ else // if no, find i in binlist
+ {
+ int iterater = cells.get(s[i].cell);
+ int pointer = iterater;
+ while ((iterater != i)&&(iterater != -1))
+ {
+ pointer = iterater;
+ iterater = binlist[iterater];
+ }
+ if (iterater == -1) // got to end of list without finding i
+ {
+ std::cout << "problem " << i << " wasn't in claimed, cell iterater = -1" << std::endl;
+ OutputCells();
+ }
+ else // we found i!
+ {
+ binlist[pointer] = binlist[i];
+ binlist[i] = -1;
+ }
+ }
+
+ // now add i to cell array at celli
+ s[i].cell = celli;
+
+ //first check to see if entry at celli
+ if (cells.get(celli) == -1) //if yes, add i to cells gridfield
+ cells.get(celli) = i;
+ else // if no, add i to right place in binlist
+ {
+ int iterater = cells.get(celli); // now iterate through to end and add i
+ int pointer = iterater;
+ while (iterater != -1) // find the end of the list
+ {
+ pointer = iterater;
+ iterater = binlist[iterater];
+ }
+ binlist[pointer] = i;
+ binlist[i] = -1; // redundant
+ }
+}
+
+
+//==============================================================
+// Output event heap...purely used for debugging
+//==============================================================
+void box::OutputEvents()
+{
+ h.print();
+}
+
+
+//==============================================================
+// Output positions of spheres and their cells...purely used for debugging
+//==============================================================
+void box::OutputCells()
+{
+ for (int i=0; i<N; i++)
+ std::cout << i << " " << s[i].x << " " << s[i].v << " " << s[i].cell << std::endl;
+}
+
+
+//==============================================================
+// Update positions...purely for graphical display
+//==============================================================
+void box::TrackPositions()
+{
+ for (int i=0; i<N; i++)
+ x[i] = s[i].x + s[i].v*(gtime-s[i].lutime);
+}
+
+
+//==============================================================
+// Computes the total energy
+//==============================================================
+double box::Energy()
+{
+ double E=0;
+ for (int i=0; i<N; i++)
+ E += 0.5*s[i].m*s[i].v.norm_squared();
+
+ return E/N;
+}
+
+
+//==============================================================
+// Calculates the packing fraction
+//==============================================================
+double box::PackingFraction()
+{
+ double rfactor = 0.;
+ for (int i=0; i<N; i++)
+ {
+ rfactor += pow(s[i].r + gtime*s[i].gr, DIM);
+ }
+ double v = (rfactor*pow(sqrt(PI), DIM))/(exp(lgamma(1.+((double)(DIM))/2.)));
+ return v/(pow(SIZE, DIM));
+}
+
+
+//==============================================================
+// Checks to make sure all sphere diameters are less than dimension
+// of unit cell
+//==============================================================
+int box::CheckSphereDiameters()
+{
+ int offender = 0;
+ for (int i=0; i<N; i++)
+ {
+ if (s[i].r*2 > 1./ngrids){
+ offender = i;
+ break;
+ }
+ }
+ return offender;
+}
+
+
+//==============================================================
+// Change ngrids
+//==============================================================
+void box::ChangeNgrids(int newngrids)
+{
+ cells.set_size(newngrids);
+ cells.initialize(-1); // initialize cells to -1
+ for (int i=0; i<N; i++) // initialize binlist to -1
+ binlist[i] = -1;
+ AssignCells();
+ for (int i=0; i<N; i++)
+ s[i].nextevent = event(0., i, INF);
+ Process(N);
+}
+
+
+//==============================================================
+// Calculates the optimal ngrids for the initial configuration
+// and assumes that ngrids gets updated (reduced) as the packing
+// proceeds
+//==============================================================
+int box::Optimalngrids2(double currentradius)
+{
+ return (int)(1./(2.*currentradius));
+}
+
+
+//==============================================================
+// Calculates the optimal ngrids, assuming ngrids is not updated
+// automatically and is very conservative
+//==============================================================
+int box::Optimalngrids()
+{
+ double maxr;
+
+ maxr = pow(exp(lgamma(1.+((double)(DIM))/2.))*maxpf/
+ (N*(bidispersityfraction + (1.-bidispersityfraction)*
+ pow(bidispersityratio, DIM))),
+ 1./DIM)/sqrt(PI);
+
+ return (int)(1./(2.*maxr));
+}
+
+
+//==============================================================
+// Processes n events
+//==============================================================
+void box::Process(int n)
+{
+ double deltat = gtime;
+ for (int i=0; i<n; i++)
+ {
+ ProcessEvent();
+ }
+ pf = PackingFraction(); // packing fraction
+ deltat = gtime - deltat;
+ double oldenergy = energy;
+ energy = Energy(); // kinetic energy
+
+ energychange = ((oldenergy - energy)/oldenergy)*100; // percent change in energy
+
+ if (deltat != 0.)
+ {
+ pressure = 1+xmomentum/(2.*energy*N*deltat);
+ collisionrate = ((double)(ncollisions))/deltat;
+ }
+
+ // reset to 0
+ ncollisions = 0;
+ ntransfers = 0;
+ nchecks = 0;
+ xmomentum = 0.;
+ ncycles++;
+}
+
+
+//==============================================================
+// Prints statistics for n events
+//==============================================================
+void box::PrintStatistics()
+{
+ std::cout << "packing fraction = " << pf << std::endl;
+ std::cout << "gtime = " << gtime << std::endl;
+ std::cout << "total time = " << rtime+gtime << std::endl;
+ std::cout << "kinetic energy = " << energy << std::endl;
+ std::cout << "total # events = " << neventstot << std::endl;
+ std::cout << "# events = " << ncollisions+ntransfers+nchecks << ", # collisions = " << ncollisions << ", # transfers = " << ntransfers << ", # checks =" << nchecks << std::endl;
+ std::cout << "growthrate = " << growthrate << std::endl;
+ std::cout << "collisionrate = " << collisionrate << std::endl;
+ std::cout << "reduced pressure = " << pressure << std::endl;
+ std::cout << "-----------------" << std::endl;
+}
+
+
+//==============================================================
+// Updates spheres to gtime, synchronizes, and can change growth rate
+//==============================================================
+void box::Synchronize(bool rescale)
+{
+ double vavg = sqrt(2.*M*energy);
+
+ for (int i=0; i<N; i++)
+ {
+ s[i].x = s[i].x + s[i].v*(gtime-s[i].lutime);
+ s[i].nextevent.time -= gtime;
+
+ if (s[i].nextevent.time < 0.)
+ std::cout << "error, event times negative after synchronization" << std::endl;
+ if (rescale == true) // give everyone checks
+ {
+ s[i].nextevent = event(0., i, INF);
+ s[i].v /= vavg;
+ }
+
+ s[i].lutime = 0.;
+ s[i].r += gtime*s[i].gr;
+ }
+
+ //r += gtime*growthrate; // r defined at gtime = 0
+ rtime += gtime;
+ gtime = 0.;
+
+ if (rescale == true)
+ Process(N);
+}
+
+
+//==============================================================
+// Run time
+//==============================================================
+void box::RunTime()
+{
+ time(&end);
+ std::cout << "run time = " << difftime(end, start) << std::endl;
+}
+
+
+//==============================================================
+// Write configuration
+//==============================================================
+void box::WriteConfiguration(double *data)
+{
+ int count1;
+
+ if (gtime != 0.) // synchronize spheres if not currently synchronized
+ Synchronize(false);
+
+
+ for (int i=0; i<N; i++) // output radius and positions
+ {
+ for (int k=0; k<DIM; k++)
+ data[DIM * i + k] = s[i].x[k];
+ }
+
+}
diff --git a/src/spheres_poly/box.h b/src/spheres_poly/box.h
new file mode 100644
index 0000000..69418f4
--- /dev/null
+++ b/src/spheres_poly/box.h
@@ -0,0 +1,169 @@
+/*
+ Packing of hard spheres via molecular dynamics
+ Developed by Monica Skoge, 2006, Princeton University
+ Contact: Aleksandar Donev (adonev@math.princeton.edu) with questions
+ This code may be used, modified and distributed freely.
+ Please cite:
+
+ "Packing Hyperspheres in High-Dimensional Euclidean Spaces"
+ M. Skoge, A. Donev, F. H. Stillinger and S. Torquato, 2006
+
+ if you use these codes.
+*/
+
+//-----------------------------------------------------------------------------
+// Box maker
+//---------------------------------------------------------------------------
+
+#ifndef BOX_H
+#define BOX_H
+
+
+#include <vector>
+#include <math.h>
+
+#include "vector.h"
+#include "grid_field.h"
+#include "event.h"
+#include "sphere.h"
+#include "heap.h"
+
+
+#define PI 3.141592653589793238462643
+#define SIZE 1.0 // size of box
+#define VOLUMESPHERE pow(PI,((double)(DIM))/2.)/exp(lgamma(1+((double)(DIM))/2.)) // volume prefactor for sphere
+#define DBL_EPSILON 2.2204460492503131e-016 // smallest # such that 1.0+DBL_EPSILON!=1.0
+#define M 1.0
+
+//---------------------------------------------------------------------------
+// Class neighbor
+//---------------------------------------------------------------------------
+class neighbor
+{
+ public:
+ int i;
+
+ neighbor(int i_i);
+
+ public:
+ virtual void Operation(int j, vector<DIM, int>& pboffset) = 0;
+};
+
+
+class box {
+
+ public:
+
+ // constructor and destructor
+ box(int N_i, double r_i, double growthrate_i, double maxpf_i,
+ double bidispersityratio, double bidispersityfraction,
+ double massratio, int hardwallBC);
+ ~box();
+
+ // Creating configurations
+ int Optimalngrids();
+ int Optimalngrids2(double maxr);
+ void CreateSpheres(double temp);
+ void CreateSphere(int Ncurrent);
+ double Velocity(double temp);
+ void VelocityGiver(double temp);
+ void SetInitialEvents();
+ void RecreateSpheres(const char* filename, double temp);
+ void ReadPositions(const char* filename);
+ void AssignCells();
+
+ // Predicting next event
+ event FindNextEvent(int i);
+ void CollisionChecker(event c);
+ event FindNextTransfer(int i);
+ event FindNextCollision(int i);
+ void ForAllNeighbors(int, vector<DIM, int>, vector<DIM,int>, neighbor&);
+ void PredictCollision(int i, int j, vector<DIM, int> pboffset,
+ double& ctime, int& cpartner,
+ vector<DIM, int>& cpartnerpboffset);
+ double CalculateCollision(int i, int j, vector<DIM> pboffset);
+ double QuadraticFormula(double a, double b, double c);
+
+ // Processing an event
+ void Process(int n);
+ void ProcessEvent();
+ void Collision(event e);
+ void Transfer(event e);
+ void UpdateCell(int i, vector<DIM,int>& celli);
+ void Synchronize(bool rescale);
+ void ChangeNgrids(int newngrids);
+
+ // Debugging
+ void TrackPositions();
+ void OutputEvents();
+ void OutputCells();
+ void GetInfo();
+ int CheckSphereDiameters();
+
+ // Statistics
+ double Energy();
+ double PackingFraction();
+ void PrintStatistics();
+ void RunTime();
+ void WriteConfiguration(double* data);
+
+
+ //variables
+
+ const int N; // number of spheres
+
+ int ngrids; // number of cells in one direction
+ double maxpf;
+ double growthrate; // growth rate of the spheres
+ double r; // radius, defined at gtime = 0
+ double gtime; // this is global clock
+ double rtime; // reset time, total time = rtime + gtime
+ double collisionrate; // average rate of collision between spheres
+ double bidispersityratio; // ratio of sphere radii
+ double bidispersityfraction; // fraction of smaller spheres
+ double massratio; // ratio of sphere masses
+ int hardwallBC; // =0 for periodic BC, =1 for hard wall
+
+
+ // statistics
+ double pressure; // pressure
+ double xmomentum; // exchanged momentum
+ double pf; // packing fraction
+ double energy; // kinetic energy
+ double energychange;
+ int ncollisions; // number of collisions
+ int ntransfers; // number of transfers
+ int nchecks; // number of checks
+ int neventstot; // total number of events
+ int ncycles; // counts # cycles for output
+
+ time_t start, error, end; // run time of program
+
+ // arrays
+ sphere *s; // array of spheres
+ grid_field<DIM, int> cells; // array that keeps track of spheres in each cell
+ int *binlist; // linked-list for cells array
+ heap h; // event heap
+ vector<DIM> *x; // positions of spheres.used for graphics
+};
+
+
+//---------------------------------------------------------------------------
+// Predicts collisions, inherits neighbor operation
+//---------------------------------------------------------------------------
+class collision : public neighbor
+{
+ public:
+
+ box *b;
+ double ctime;
+ int cpartner;
+ vector<DIM,int> cpartnerpboffset;
+
+ public:
+ collision(int i_i, box *b);
+
+ virtual void Operation(int j, vector<DIM, int>& pboffset);
+};
+
+#endif
diff --git a/src/spheres_poly/event.cpp b/src/spheres_poly/event.cpp
new file mode 100644
index 0000000..c8c104d
--- /dev/null
+++ b/src/spheres_poly/event.cpp
@@ -0,0 +1,65 @@
+#include "event.h"
+
+//==============================================================
+//==============================================================
+// Class Event
+//==============================================================
+//==============================================================
+
+
+//==============================================================
+// Constructor
+//==============================================================
+event::event(double time_i, int i_i, int j_i, vector<DIM,int> v_i):
+ time(time_i),
+ i(i_i),
+ j(j_i),
+ v(v_i)
+{
+}
+
+event::event(double time_i, int i_i, int j_i):
+ time(time_i),
+ i(i_i),
+ j(j_i)
+{
+}
+
+event::event(const event& e)
+{
+ time = e.time;
+ i = e.i;
+ j = e.j;
+ v = e.v;
+}
+
+event::event()
+{
+}
+
+
+//==============================================================
+// Destructor
+//==============================================================
+event::~event()
+{
+}
+
+void event::erase()
+{
+ time = dblINF;
+ i = 0;
+ j = 0;
+}
+
+bool event::operator<(const event& e) const
+{
+ return e.time < time;
+}
+
+bool event::operator>(const event& e) const
+{
+ return e.time > time;
+}
+
+
diff --git a/src/spheres_poly/event.h b/src/spheres_poly/event.h
new file mode 100644
index 0000000..55dd1fc
--- /dev/null
+++ b/src/spheres_poly/event.h
@@ -0,0 +1,58 @@
+#ifndef EVENT_H
+#define EVENT_H
+
+#include "vector.h"
+#define INF 100000000
+#define dblINF 100000000.
+
+class event {
+
+ public:
+
+ // constructor and destructor
+ event(double time_i, int i_i, int j_i, vector<DIM,int> v_i);
+ event(double time_i, int i_i, int j_i);
+ event(const event& e);
+ event();
+
+ ~event();
+
+ bool operator<(const event&) const;
+ bool operator>(const event&) const;
+ void erase();
+
+ //variables
+
+
+ double time; // time of next collision
+ int i; // collision partner with lower number
+ int j; // collision partner with higher number
+ vector<DIM,int> v; // virtual image
+
+ /* 0<=j<=N binary collision between i and j
+ j=N+DIM+1+x transfer where x=-(k+1) for left wall
+ and x=k+1 for right wall
+ j=INF both check after event that did not altered motion of i and check after event that altered motion of i, i.e rescaling of velocities. I currently don't see need t o separate the two
+
+ j=-1 check after collision
+
+ Virtual identifiers as scalars...I think bad idea, but here's my work
+ there will be easily fixed problems if DIM>=10
+ -x<=v<=x x=k+1, image in k direction
+ v=xy x,y
+ =-xy -x,y
+ =-yx x,-y
+ =yx -x,-y
+ v=xyz x,y,z
+ =-xyz -x,y,z
+ =-yxz x,-y,z
+ =-zxy x,y,-z
+ =zyx -x,-y,z
+ =xzy x,-y,-z
+ =yzx -x,y,-z
+ =-zyx -x,-y,-z
+ */
+
+};
+
+#endif
diff --git a/src/spheres_poly/grid_field.h b/src/spheres_poly/grid_field.h
new file mode 100644
index 0000000..10430ed
--- /dev/null
+++ b/src/spheres_poly/grid_field.h
@@ -0,0 +1,148 @@
+/*
+ Packing of hard spheres via molecular dynamics
+ Developed by Monica Skoge, 2006, Princeton University
+ Contact: Aleksandar Donev (adonev@math.princeton.edu) with questions
+ This code may be used, modified and distributed freely.
+ Please cite:
+
+ "Packing Hyperspheres in High-Dimensional Euclidean Spaces"
+ M. Skoge, A. Donev, F. H. Stillinger and S. Torquato, 2006
+
+ if you use these codes.
+*/
+
+
+#ifndef GRID_FIELD_H
+#define GRID_FIELD_H
+
+#include "vector.h"
+
+// ======================================================================
+// grid_field
+// ======================================================================
+
+// A field of V-vectors on a D dimensional manifold
+
+template<int D, class T>
+class grid_field {
+
+ public:
+ int elements;
+
+ private:
+ T* f;
+ vector<D, int> size; // number of grid points for each dimension
+ vector<D, int> offset;
+
+ public:
+
+ grid_field();
+ grid_field(const vector<D, int>&);
+ ~grid_field();
+
+ T& get(const vector<D, int>&);
+
+ vector<D, int> get_size() const;
+ void set_size(const vector<D, int>&);
+ void set_size(const int);
+ void initialize(const int i);
+};
+
+
+// grid_field
+// ~~~~~~~~~~~~
+template<int D, class T>
+grid_field<D, T>::grid_field()
+ : f(0), elements(0)
+{
+}
+
+
+// grid_field
+// ~~~~~~~~~~~~
+template<int D, class T>
+grid_field<D, T>::grid_field(const vector<D, int>& s)
+ : f(0)
+{
+ set_size(s);
+}
+
+
+// ~grid_field
+// ~~~~~~~~~~~~~
+template <int D, class T>
+grid_field<D, T>::~grid_field()
+{
+ if(f != 0)
+ delete[] f;
+}
+
+
+// get_size
+// ~~~~~~~~
+template<int D, class T>
+inline vector<D, int> grid_field<D, T>::get_size() const
+{
+ return size;
+}
+
+
+// set_size
+// ~~~~~~~~
+template<int D, class T>
+void grid_field<D, T>::set_size(const vector<D, int>& s)
+{
+ if(f != 0)
+ delete[] f;
+
+ size = s;
+
+ elements = 1;
+ for(int i=0; i<D; i++) {
+ offset[i] = elements;
+ elements *= size.x[i];
+ }
+
+ f = new T[elements];
+}
+
+
+// set_size
+// ~~~~~~~~
+template<int D, class T>
+void grid_field<D, T>::set_size(const int s)
+{
+ vector<D, int> square;
+
+ for(int k=0; k<D; k++)
+ square[k] = s;
+
+ set_size(square);
+}
+
+
+// get
+// ~~~
+template<int D, class T>
+inline T& grid_field<D, T>::get(const vector<D, int>& pos)
+{
+ int p=0;
+ for(int i=0; i<D; i++)
+ p += pos.x[i]*offset[i];
+
+ return f[p];
+}
+
+
+// initialize
+// ~~~
+template<int D, class T>
+void grid_field<D, T>::initialize(const int value)
+{
+ for(int i=0; i<elements; i++)
+ f[i] = value;
+}
+
+
+
+#endif
diff --git a/src/spheres_poly/heap.cpp b/src/spheres_poly/heap.cpp
new file mode 100644
index 0000000..c83865c
--- /dev/null
+++ b/src/spheres_poly/heap.cpp
@@ -0,0 +1,186 @@
+#include "heap.h"
+#include "event.h"
+#include <iostream>
+#include "box.h"
+
+
+//==============================================================
+//==============================================================
+// Class heap: Event heap used in hard spheres calculation
+//==============================================================
+//==============================================================
+
+
+//==============================================================
+// Constructor
+//==============================================================
+heap::heap(int maxsize_i):
+ maxsize(maxsize_i)
+{
+ a = new int[maxsize];
+ index = new int[maxsize];
+
+ N = 0; // current number of events in heap
+}
+
+
+//==============================================================
+// Constructor
+//==============================================================
+heap::heap(const heap &h)
+{
+ maxsize = h.maxsize;
+ a = h.a;
+ index = h.index;
+ N = h.N; // current number of events in heap
+ s = h.s;
+}
+
+
+//==============================================================
+// Destructor
+//==============================================================
+heap::~heap()
+{
+ delete[] a;
+ delete[] index;
+}
+
+//==============================================================
+// Upheap
+//==============================================================
+void heap::upheap(int k)
+{
+ int i = a[k];
+
+ while ((k>1)&&(s[a[k/2]].nextevent.time > s[i].nextevent.time))
+ {
+ a[k] = a[k/2];
+ index[a[k/2]] = k;
+ k = k/2;
+ }
+ a[k] = i;
+ index[i] = k;
+}
+
+//==============================================================
+// Downheap
+//==============================================================
+void heap::downheap(int k)
+{
+ int j;
+ int i = a[k];
+
+ while(k <= N/2)
+ {
+ j = k+k;
+ if ((j < N)&&(s[a[j]].nextevent.time > s[a[j+1]].nextevent.time))
+ j++;
+ if (s[i].nextevent.time <= s[a[j]].nextevent.time)
+ break;
+ a[k] = a[j];
+ index[a[j]] = k;
+ k = j;
+ }
+ a[k] = i;
+ index[i] = k;
+}
+
+//==============================================================
+// Insert
+//==============================================================
+void heap::insert(int i)
+{
+ if (N >= maxsize)
+ std::cout << "error, N >= maxsize, cannot insert another event" << std::endl;
+ else
+ {
+ N++;
+ a[N] = i;
+ index[i] = N;
+ upheap(N);
+ }
+}
+
+
+//==============================================================
+// Extract max
+//==============================================================
+int heap::extractmax()
+{
+ return a[1];
+}
+
+/*
+//==============================================================
+// Replace
+//==============================================================
+void heap::replace(int i)
+{
+ a[1] = i;
+ s[i].nextevent = t;
+
+ if (!(e.time > s[a[1]].nextevent.time))
+ {
+ if (!(s[a[1]].nextevent.j == INF))// must be check i'm changing to coll. at same time
+ std::cout << "error replaced non check with earlier time" << std::endl;
+ a[1] = i;
+ index[i] = 1;
+ }
+ else
+ {
+ a[0] = i;
+ downheap(0);
+ }
+}
+*/
+
+//==============================================================
+// Search
+//==============================================================
+int heap::search(int j)
+{
+ for (int k=1; k<=N; k++)
+ {
+ if (a[k] == j)
+ return k;
+ }
+ return -1;
+}
+
+/*
+//==============================================================
+// Change
+//==============================================================
+void heap::change(int i)
+{
+ int iindex = index[i];
+
+ if (s[i].nextevent.time == s[a[iindex]].nextevent.time)
+ std::cout << "error changing an event to an equal time" << std::endl;
+ else if (s[i].nextevent.time > s[a[iindex]].nextevent.time)
+ std::cout << "error changing an event to a greater time" << std::endl;
+
+ a[iindex] = i;
+ upheap(iindex);
+}
+*/
+
+//==============================================================
+// Print
+//==============================================================
+void heap::print()
+{
+ for (int k=1; k<=N; k++)
+ std::cout << k << " " << a[k] << " " << s[a[k]].nextevent.j << " " << s[a[k]].nextevent.time << std::endl;
+}
+
+
+//==============================================================
+// Check index array
+//==============================================================
+void heap::checkindex()
+{
+ for (int k=1; k<=N; k++)
+ std::cout << k << " " << a[k] << " " << index[a[k]] << std::endl;
+}
diff --git a/src/spheres_poly/heap.h b/src/spheres_poly/heap.h
new file mode 100644
index 0000000..91180e3
--- /dev/null
+++ b/src/spheres_poly/heap.h
@@ -0,0 +1,42 @@
+//---------------------------------------------------------------------------
+// Event heap maker
+//---------------------------------------------------------------------------
+
+#ifndef HEAP_H
+#define HEAP_H
+
+#include "event.h"
+#include "sphere.h"
+
+class heap {
+
+ public:
+
+ // constructor and destructor
+ heap(int maxsize);
+ heap(const heap &h);
+ ~heap();
+
+ // variables
+ int maxsize; // max allowed number of events
+ int N; // current number of events
+ int *a;
+ sphere *s;
+ int *index; // array of indices for each sphere
+ //event minevent;
+
+
+ // functions which operate on a binary heap
+
+ void upheap(int k);
+ void downheap(int k);
+ void insert(int i);
+ void replace(int i);
+ int search(int j);
+ void change(int i);
+ int extractmax();
+ void print();
+ void checkindex();
+
+};
+#endif
diff --git a/src/spheres_poly/neighbor.cpp b/src/spheres_poly/neighbor.cpp
new file mode 100644
index 0000000..920c099
--- /dev/null
+++ b/src/spheres_poly/neighbor.cpp
@@ -0,0 +1,40 @@
+#include "vector.h"
+#include <iostream>
+#include "box.h"
+
+
+//==============================================================
+//==============================================================
+// Class neighbor
+//==============================================================
+//==============================================================
+
+neighbor::neighbor(int i_i)
+ : i(i_i)
+{
+}
+
+
+//==============================================================
+//==============================================================
+// Class collision
+//==============================================================
+//==============================================================
+
+collision::collision(int i_i, box *b_i)
+ : neighbor(i_i),
+ b(b_i)
+{
+ ctime = dblINF;
+ cpartner = i;
+}
+
+
+//==============================================================
+// Operation is finding the next collision from a given cell
+//==============================================================
+void collision::Operation(int j, vector<DIM, int>& pboffset)
+{
+ b->PredictCollision(i, j, pboffset, ctime, cpartner, cpartnerpboffset);
+}
+
diff --git a/src/spheres_poly/sphere.cpp b/src/spheres_poly/sphere.cpp
new file mode 100644
index 0000000..0fbbd68
--- /dev/null
+++ b/src/spheres_poly/sphere.cpp
@@ -0,0 +1,66 @@
+#include "box.h"
+#include <iostream>
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+
+#include "vector.h"
+
+//==============================================================
+//==============================================================
+// Class Sphere:
+//==============================================================
+//==============================================================
+
+
+//==============================================================
+// Constructor
+//==============================================================
+sphere::sphere()
+{
+}
+
+
+//==============================================================
+// Constructor
+//==============================================================
+sphere::sphere(const sphere& s)
+{
+ i = s.i;
+ x = s.x;
+ v = s.v;
+ cell = s.cell;
+ lutime = s.lutime;
+ nextevent = s.nextevent;
+ nextcollision = s.nextcollision;
+ r = s.r;
+ gr = s.gr;
+ m = s.m;
+ species=s.species;
+}
+
+//==============================================================
+// Constructor
+//==============================================================
+sphere::sphere(int i_i, vector<DIM> x_i, vector<DIM, int> cell_i,
+ double lutime_i, double r_i, double gr_i, double m_i, int species_i):
+ i(i_i),
+ x(x_i),
+ cell(cell_i),
+ lutime(lutime_i),
+ r(r_i),
+ gr(gr_i),
+ m(m_i),
+ species(species_i)
+{
+}
+
+//==============================================================
+// Destructor
+//==============================================================
+sphere::~sphere()
+{
+
+}
+
+
diff --git a/src/spheres_poly/sphere.h b/src/spheres_poly/sphere.h
new file mode 100644
index 0000000..28c64b4
--- /dev/null
+++ b/src/spheres_poly/sphere.h
@@ -0,0 +1,44 @@
+#ifndef SPHERE_H
+#define SPHERE_H
+
+
+#include "vector.h"
+
+
+class sphere {
+
+ public:
+
+ // constructor and destructor
+
+ sphere();
+ sphere(const sphere& s);
+ sphere(int i_i, vector<DIM> x, vector<DIM, int> cell_i, double lutime_i,
+ double r_i, double gr_i, double m_i, int species_i);
+ ~sphere();
+
+ //variables
+
+ int i; // sphere ID
+
+ // impending event
+ event nextevent; // next event...can be collision or transfer
+ event nextcollision; // next collision if next event is transfer
+ // maybe nextnext event
+
+ // past information
+ double lutime; // last update time
+ vector<DIM, int> cell; // cell that it belongs to
+ vector<DIM, double> x; // position
+ vector<DIM, double> v; // velocity
+ double r; // sphere radius
+ double gr; // sphere growth rate
+ double m; // sphere mass
+ int species; // species number (not used during the MD)
+ // make sure efficent in memory
+
+
+
+};
+
+#endif
diff --git a/src/spheres_poly/spheres.cpp b/src/spheres_poly/spheres.cpp
new file mode 100644
index 0000000..9b14b01
--- /dev/null
+++ b/src/spheres_poly/spheres.cpp
@@ -0,0 +1,64 @@
+//===========================================================
+//===========================================================
+//===========================================================
+//
+// Molecular dynamics simulation of hardspheres
+//
+//===========================================================
+//===========================================================
+//===========================================================
+
+#include <iostream>
+#include <math.h>
+#include <fstream>
+#include <vector>
+#include <time.h>
+#include <string.h>
+#include <stdlib.h>
+#include <cstdlib>
+
+#include "box.h"
+#include "sphere.h"
+#include "event.h"
+#include "heap.h"
+
+
+double *spherespp(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate)
+{
+ int eventspercycle = 20; // # events per sphere per cycle
+ double initialpf = 0.01; // initial packing fraction
+ double temp = 0.2; // initial temp., use 0 for zero velocities
+ double growthrate = 0.001; // growth rate
+ double massratio = 1.; // ratio of sphere masses
+ int hardwallBC = 0; // 0 for periodic, 1 for hard wall BC
+
+ double d, r; // initial diameter and radius of spheres
+
+ r = pow(initialpf*pow(SIZE, DIM)/(N*VOLUMESPHERE), 1.0/((double)(DIM)));
+
+ box b(N, r, growthrate, maxpf, bidispersityratio,
+ bidispersityfraction, massratio, hardwallBC);
+
+ b.CreateSpheres(temp);
+
+
+ while ((b.collisionrate < maxcollisionrate) && (b.pf < maxpf) && (b.pressure < maxpressure))
+ {
+ b.Process(eventspercycle*N);
+ b.Synchronize(true);
+ }
+
+ double *data = (double *)malloc(DIM * N * sizeof(double));
+
+ b.WriteConfiguration(data);
+
+ return data;
+}
+
+extern "C" {
+ double *spheres(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate) {
+ return spherespp(N, bidispersityratio, bidispersityfraction, maxpf, maxpressure, maxcollisionrate);
+ }
+}
+
+
diff --git a/src/spheres_poly/spheres.h b/src/spheres_poly/spheres.h
new file mode 100644
index 0000000..44d8052
--- /dev/null
+++ b/src/spheres_poly/spheres.h
@@ -0,0 +1,4 @@
+
+#pragma once
+
+double *spheres(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate);
diff --git a/src/spheres_poly/vector.h b/src/spheres_poly/vector.h
new file mode 100644
index 0000000..9ee481f
--- /dev/null
+++ b/src/spheres_poly/vector.h
@@ -0,0 +1,471 @@
+#ifndef VECTOR_H
+#define VECTOR_H
+
+#include <iostream>
+#include <fstream>
+
+#define DIM 2
+
+template <int D, typename T=double>
+class vector {
+
+ public:
+ T x[D];
+
+ public:
+
+ vector();
+ vector(const T[D]);
+ vector(const vector&);
+ ~vector();
+
+ vector<D, T>& operator+=(const vector<D, T>&);
+ vector<D, T>& operator-=(const vector<D, T>&);
+ vector<D, T>& operator*=(const T);
+ vector<D, T>& operator/=(const T);
+ vector<D, T> operator+(const vector<D, T>&) const;
+ vector<D, T> operator-(const vector<D, T>&) const;
+ vector<D, T> operator*(const T) const;
+ vector<D, T> operator/(const T) const;
+ vector<D, T> operator%(const T) const;
+ bool operator==(const vector<D, T> &a) const;
+
+ vector<D, int> integer() const;
+ vector<D, double> Double() const;
+ static vector<D, int> integer(const vector<D, T>&);
+ static vector<D, double> Double(const vector<D, T>&);
+
+ T& operator[](const unsigned int);
+
+ double dot(const vector<D, T>&) const;
+ static double dot(const vector<D, T>&, const vector<D, T>&);
+
+ double norm_squared() const;
+ static double norm_squared(const vector<D, T>&);
+
+ void read(std::ifstream&);
+ void write(std::ofstream&) const;
+};
+
+
+template <int D, typename T>
+std::ostream& operator<<(std::ostream&, const vector<D, T>&);
+
+
+// constructor
+// ~~~~~~~~~~~
+template <int D, typename T>
+vector<D, T>::vector()
+{
+ for(int k=0; k<D; k++)
+ x[k] = 0;
+}
+
+template <int D, typename T>
+vector<D, T>::vector(const T x_i[D])
+{
+ for(int k=0; k<D; k++)
+ x[k] = x_i[k];
+}
+
+template <int D, typename T>
+vector<D, T>::vector(const vector<D, T> &v)
+{
+ for(int k=0; k<D; k++)
+ x[k] = v.x[k];
+}
+
+
+// destructor
+// ~~~~~~~~~~
+template <int D, typename T>
+vector<D, T>::~vector()
+{
+}
+
+
+// +=
+// ~~
+template <int D, typename T>
+inline vector<D, T>& vector<D, T>::operator+=(const vector<D, T> &v)
+{
+ for(int k=0; k<D; k++)
+ x[k] += v.x[k];
+
+ return *this;
+}
+
+
+// -=
+// ~~
+template <int D, typename T>
+inline vector<D, T>& vector<D, T>::operator-=(const vector<D, T> &v)
+{
+ for(int k=0; k<D; k++)
+ x[k] -= v.x[k];
+
+ return *this;
+}
+
+
+// *=
+// ~~
+template <int D, typename T>
+inline vector<D, T>& vector<D, T>::operator*=(const T s)
+{
+ for(int k=0; k<D; k++)
+ x[k] *= s;
+
+ return *this;
+}
+
+// /=
+// ~~
+template <int D, typename T>
+inline vector<D, T>& vector<D, T>::operator/=(const T s)
+{
+ for(int k=0; k<D; k++)
+ x[k] /= s;
+
+ return *this;
+}
+
+
+// +
+// ~
+template <int D, typename T>
+inline vector<D, T> vector<D, T>::operator+(const vector<D, T> &a) const
+{
+ vector<D, T> c;
+
+ for(int k=0; k<D; k++)
+ c.x[k] = x[k] + a.x[k];
+
+ return c;
+}
+
+
+// -
+// ~
+template <int D, typename T>
+inline vector<D, T> vector<D, T>::operator-(const vector<D, T> &a) const
+{
+ vector<D, T> c;
+
+ for(int k=0; k<D; k++)
+ c[k] = x[k] - a.x[k];
+
+ return c;
+}
+
+
+// *
+// ~
+template <int D, typename T>
+inline vector<D, T> vector<D, T>::operator*(const T s) const
+{
+ vector<D, T> c;
+
+ for(int k=0; k<D; k++)
+ c[k] = x[k] * s;
+
+ return c;
+}
+
+
+// /
+// ~
+template <int D, typename T>
+inline vector<D, T> vector<D, T>::operator/(const T s) const
+{
+ vector<D, T> c;
+
+ for(int k=0; k<D; k++)
+ c[k] = x[k] / s;
+
+ return c;
+}
+
+
+// ==
+// ~
+template <int D, typename T>
+inline bool vector<D, T>::operator==(const vector<D, T> &a) const
+{
+ for(int k=0; k<D; k++)
+ {
+ if (!(x[k]==a.x[k]))
+ return false;
+ }
+ return true;
+}
+
+
+// %
+// ~
+template <int D, typename T>
+inline vector<D, T> vector<D, T>::operator%(const T s) const
+{
+ vector<D, T> c;
+
+ for(int k=0; k<D; k++)
+ c[k] = x[k] % s;
+
+ return c;
+}
+
+
+// integer
+// ~~~~~~~
+template <int D, typename T>
+inline vector<D, int> vector<D, T>::integer() const
+{
+ vector<D, int> c;
+
+ for(int k=0; k<D; k++)
+ c[k] = (int)x[k];
+
+ return c;
+}
+
+template <int D, typename T>
+inline vector<D, int> vector<D, T>::integer(const vector<D, T>& v)
+{
+ return v.integer();
+}
+
+
+// double
+// ~~~~~~~
+template <int D, typename T>
+inline vector<D, double> vector<D, T>::Double() const
+{
+ vector<D, double> c;
+
+ for(int k=0; k<D; k++)
+ c[k] = (double)x[k];
+
+ return c;
+}
+
+template <int D, typename T>
+inline vector<D, double> vector<D, T>::Double(const vector<D, T>& v)
+{
+ return v.Double();
+}
+
+
+
+// []
+// ~~
+template <int D, typename T>
+inline T& vector<D, T>::operator[](const unsigned int i)
+{
+ return x[i];
+}
+
+
+// Dot
+// ~~~
+template <int D, typename T>
+inline double vector<D, T>::dot(const vector<D, T> &a) const
+{
+ double d=0;
+
+ for(int k=0; k<D; k++)
+ d += x[k] * a.x[k];
+
+ return d;
+}
+
+template <int D, typename T>
+inline double vector<D, T>::dot(const vector<D, T> &a, const vector<D, T> &b)
+{
+ return a.dot(b);
+}
+
+
+// NormSquared
+// ~~~~~~~~~~~
+template <int D, typename T>
+inline double vector<D, T>::norm_squared() const
+{
+ return dot(*this, *this);
+}
+
+template <int D, typename T>
+inline double vector<D, T>::norm_squared(const vector<D, T>& v)
+{
+ return v.norm_squared();
+}
+
+
+// read
+// ~~~~
+template <int D, typename T>
+void vector<D, T>::read(std::ifstream& in)
+{
+ in.read((char*)x, sizeof(T)*D);
+}
+
+// write
+// ~~~~~
+template <int D, typename T>
+void vector<D, T>::write(std::ofstream& out) const
+{
+ out.write((const char*)x, sizeof(T)*D);
+}
+
+
+
+// Insertion
+// ~~~~~~~~~
+template <int D, typename T>
+std::ostream& operator<<(std::ostream& os, const vector<D, T>& v)
+{
+ os << "(";
+
+ for(int k=0; k<D-1; k++)
+ os << v.x[k] << ", ";
+
+ os << v.x[D-1] << ")";
+
+ return os;
+}
+
+
+// ======================================================================
+// vector_field
+// ======================================================================
+
+// A field of V-vectors on a D dimensional manifold
+
+template<int V, int D, typename T=double>
+class vector_field {
+
+ public:
+ int elements;
+
+ private:
+ vector<V, T>* f;
+ vector<D, int> size; // number of grid points for each dimension
+ vector<D, int> offset;
+
+ public:
+
+ vector_field();
+ vector_field(const vector<D, int>&);
+ ~vector_field();
+
+ vector<D, int> get_size() const;
+ void set_size(const vector<D, int>&);
+
+ vector<V, T>& get(const vector<D, int>&);
+
+ void read(std::ifstream&);
+ void write(std::ofstream&) const;
+
+ static void swap(vector_field<V, D, T>&, vector_field<V, D, T>&);
+};
+
+
+// vector_field
+// ~~~~~~~~~~~~
+template<int V, int D, typename T>
+vector_field<V, D, T>::vector_field()
+ : f(0), elements(0)
+{
+}
+
+
+// vector_field
+// ~~~~~~~~~~~~
+template<int V, int D, typename T>
+vector_field<V, D, T>::vector_field(const vector<D, int>& s)
+ : f(0)
+{
+ set_size(s);
+}
+
+// ~vector_field
+// ~~~~~~~~~~~~~
+template <int V, int D, typename T>
+vector_field<V, D, T>::~vector_field()
+{
+ if(f != 0)
+ delete[] f;
+}
+
+// get_size
+// ~~~~~~~~
+template<int V, int D, typename T>
+inline vector<D, int> vector_field<V, D, T>::get_size() const
+{
+ return size;
+}
+
+// set_size
+// ~~~~~~~~
+template<int V, int D, typename T>
+void vector_field<V, D, T>::set_size(const vector<D, int>& s)
+{
+ if(f != 0)
+ delete[] f;
+
+ size = s;
+
+ elements = 1;
+ for(int i=0; i<D; i++) {
+ offset[i] = elements;
+ elements *= size.x[i];
+ }
+
+ f = new vector<V, T>[elements];
+}
+
+// get
+// ~~~
+template<int V, int D, typename T>
+inline vector<V, T>& vector_field<V, D, T>::get(const vector<D, int>& pos)
+{
+ int p=0;
+ for(int i=0; i<D; i++)
+ p += pos.x[i]*offset[i];
+
+ return f[p];
+}
+
+
+// read
+// ~~~~
+template<int V, int D, typename T>
+void vector_field<V, D, T>::read(std::ifstream& in)
+{
+ in.read((char*)f, elements*sizeof(T)*V);
+}
+
+
+// write
+// ~~~~~
+template<int V, int D, typename T>
+void vector_field<V, D, T>::write(std::ofstream& out) const
+{
+ out.write((const char*)f, elements*sizeof(T)*V);
+}
+
+// swap
+// ~~~~
+template<int V, int D, typename T>
+void vector_field<V, D, T>::swap(vector_field<V, D, T>& v1,
+ vector_field<V, D, T>& v2)
+{
+ vector<V, T>* f;
+
+ f = v1.f;
+ v1.f = v2.f;
+ v2.f = f;
+}
+
+
+
+#endif