summaryrefslogtreecommitdiff
path: root/src/spheres_poly
diff options
context:
space:
mode:
Diffstat (limited to 'src/spheres_poly')
-rw-r--r--src/spheres_poly/box.cpp1157
-rw-r--r--src/spheres_poly/box.h169
-rw-r--r--src/spheres_poly/event.cpp65
-rw-r--r--src/spheres_poly/event.h58
-rw-r--r--src/spheres_poly/grid_field.h148
-rw-r--r--src/spheres_poly/heap.cpp186
-rw-r--r--src/spheres_poly/heap.h42
-rw-r--r--src/spheres_poly/neighbor.cpp40
-rw-r--r--src/spheres_poly/sphere.cpp66
-rw-r--r--src/spheres_poly/sphere.h44
-rw-r--r--src/spheres_poly/spheres.cpp64
-rw-r--r--src/spheres_poly/spheres.h4
-rw-r--r--src/spheres_poly/vector.h471
13 files changed, 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