added support for hyperuniform lattices using existing hard-sphere jamming routine

1 files changed, 1157 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];
+   }
+ 
+}