changed code to rely on jst

13 files changed, 0 insertions, 2514 deletions

diff --git a/src/spheres_poly/box.cpp b/src/spheres_poly/box.cpp
deleted file mode 100644
index 716090d..0000000
--- a/src/spheres_poly/box.cpp
+++ /dev/null
@@ -1,1157 +0,0 @@
-/* 
-   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
deleted file mode 100644
index 69418f4..0000000
--- a/src/spheres_poly/box.h
+++ /dev/null
@@ -1,169 +0,0 @@
-/* 
-   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
deleted file mode 100644
index c8c104d..0000000
--- a/src/spheres_poly/event.cpp
+++ /dev/null
@@ -1,65 +0,0 @@
-#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
deleted file mode 100644
index 55dd1fc..0000000
--- a/src/spheres_poly/event.h
+++ /dev/null
@@ -1,58 +0,0 @@
-#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
deleted file mode 100644
index 10430ed..0000000
--- a/src/spheres_poly/grid_field.h
+++ /dev/null
@@ -1,148 +0,0 @@
-/* 
-   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
deleted file mode 100644
index c83865c..0000000
--- a/src/spheres_poly/heap.cpp
+++ /dev/null
@@ -1,186 +0,0 @@
-#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
deleted file mode 100644
index 91180e3..0000000
--- a/src/spheres_poly/heap.h
+++ /dev/null
@@ -1,42 +0,0 @@
-//---------------------------------------------------------------------------
-// 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
deleted file mode 100644
index 920c099..0000000
--- a/src/spheres_poly/neighbor.cpp
+++ /dev/null
@@ -1,40 +0,0 @@
-#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
deleted file mode 100644
index 0fbbd68..0000000
--- a/src/spheres_poly/sphere.cpp
+++ /dev/null
@@ -1,66 +0,0 @@
-#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
deleted file mode 100644
index 28c64b4..0000000
--- a/src/spheres_poly/sphere.h
+++ /dev/null
@@ -1,44 +0,0 @@
-#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
deleted file mode 100644
index 9b14b01..0000000
--- a/src/spheres_poly/spheres.cpp
+++ /dev/null
@@ -1,64 +0,0 @@
-//===========================================================
-//===========================================================
-//===========================================================
-//
-//  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
deleted file mode 100644
index 44d8052..0000000
--- a/src/spheres_poly/spheres.h
+++ /dev/null
@@ -1,4 +0,0 @@
-
-#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
deleted file mode 100644
index 9ee481f..0000000
--- a/src/spheres_poly/vector.h
+++ /dev/null
@@ -1,471 +0,0 @@
-#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