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

1 files changed, 169 insertions, 0 deletions

diff --git a/src/spheres_poly/box.h b/src/spheres_poly/box.h
new file mode 100644
index 0000000..69418f4
--- /dev/null
+++ b/src/spheres_poly/box.h
@@ -0,0 +1,169 @@
+/* 
+   Packing of hard spheres via molecular dynamics
+   Developed by Monica Skoge, 2006, Princeton University
+   Contact: Aleksandar Donev (adonev@math.princeton.edu) with questions
+   This code may be used, modified and distributed freely.
+   Please cite:
+   
+   "Packing Hyperspheres in High-Dimensional Euclidean Spaces"
+   	M. Skoge, A. Donev, F. H. Stillinger and S. Torquato, 2006
+	
+   if you use these codes.	
+*/
+
+//-----------------------------------------------------------------------------
+// Box maker
+//---------------------------------------------------------------------------
+
+#ifndef  BOX_H
+#define  BOX_H
+
+
+#include <vector>
+#include <math.h>
+
+#include "vector.h"
+#include "grid_field.h"
+#include "event.h"
+#include "sphere.h"
+#include "heap.h"
+
+
+#define PI     3.141592653589793238462643
+#define SIZE   1.0            // size of box
+#define VOLUMESPHERE pow(PI,((double)(DIM))/2.)/exp(lgamma(1+((double)(DIM))/2.)) // volume prefactor for sphere
+#define DBL_EPSILON  2.2204460492503131e-016 // smallest # such that 1.0+DBL_EPSILON!=1.0
+#define M 1.0
+
+//---------------------------------------------------------------------------
+// Class neighbor
+//---------------------------------------------------------------------------
+class neighbor
+{
+ public:
+  int i;
+  
+  neighbor(int i_i);
+
+ public:
+  virtual void Operation(int j, vector<DIM, int>& pboffset) = 0;
+};
+
+
+class box {
+  
+ public:
+  
+  // constructor and destructor
+  box(int N_i, double r_i, double growthrate_i, double maxpf_i, 
+      double bidispersityratio, double bidispersityfraction, 
+      double massratio, int hardwallBC);
+  ~box();
+
+  // Creating configurations
+  int Optimalngrids();
+  int Optimalngrids2(double maxr);
+  void CreateSpheres(double temp);
+  void CreateSphere(int Ncurrent);   
+  double Velocity(double temp);
+  void VelocityGiver(double temp);  
+  void SetInitialEvents();
+  void RecreateSpheres(const char* filename, double temp);
+  void ReadPositions(const char* filename);
+  void AssignCells();
+ 
+  // Predicting next event
+  event FindNextEvent(int i);
+  void CollisionChecker(event c);
+  event FindNextTransfer(int i);
+  event FindNextCollision(int i);
+  void ForAllNeighbors(int, vector<DIM, int>, vector<DIM,int>, neighbor&);
+  void PredictCollision(int i, int j, vector<DIM, int> pboffset, 
+			double& ctime, int& cpartner, 
+			vector<DIM, int>& cpartnerpboffset);
+  double CalculateCollision(int i, int j, vector<DIM>  pboffset);
+  double QuadraticFormula(double a, double b, double c);
+  
+  // Processing an event
+  void Process(int n);
+  void ProcessEvent();
+  void Collision(event e);
+  void Transfer(event e);
+  void UpdateCell(int i, vector<DIM,int>& celli);
+  void Synchronize(bool rescale);
+  void ChangeNgrids(int newngrids);
+
+  // Debugging
+  void TrackPositions();
+  void OutputEvents();
+  void OutputCells();
+  void GetInfo();  
+  int CheckSphereDiameters();
+
+  // Statistics
+  double Energy();
+  double PackingFraction();  
+  void PrintStatistics();
+  void RunTime();
+  void WriteConfiguration(double* data);
+  
+  
+  //variables
+
+  const int N;                   // number of spheres
+  
+  int ngrids;                    // number of cells in one direction
+  double maxpf;
+  double growthrate;             // growth rate of the spheres
+  double r;                      // radius, defined at gtime = 0
+  double gtime;                  // this is global clock
+  double rtime;                  // reset time, total time = rtime + gtime
+  double collisionrate;          // average rate of collision between spheres
+  double bidispersityratio;      // ratio of sphere radii
+  double bidispersityfraction;   // fraction of smaller spheres
+  double massratio;              // ratio of sphere masses
+  int hardwallBC;                // =0 for periodic BC, =1 for hard wall
+  
+
+  // statistics
+  double pressure;               // pressure
+  double xmomentum;              // exchanged momentum
+  double pf;                     // packing fraction
+  double energy;                 // kinetic energy
+  double energychange;
+  int ncollisions;               // number of collisions
+  int ntransfers;                // number of transfers
+  int nchecks;                   // number of checks
+  int neventstot;                // total number of events 
+  int ncycles;                   // counts # cycles for output
+
+  time_t start, error, end;      // run time of program
+
+  // arrays
+  sphere *s;                      // array of spheres
+  grid_field<DIM, int> cells; // array that keeps track of spheres in each cell
+  int *binlist;                   // linked-list for cells array
+  heap h;                         // event heap
+  vector<DIM> *x;                 // positions of spheres.used for graphics
+};
+
+
+//---------------------------------------------------------------------------
+// Predicts collisions, inherits neighbor operation
+//---------------------------------------------------------------------------
+class collision : public neighbor 
+{
+ public:
+  
+  box *b; 
+  double ctime;
+  int cpartner;
+  vector<DIM,int> cpartnerpboffset;
+
+ public:
+  collision(int i_i, box *b);
+
+  virtual void Operation(int j, vector<DIM, int>& pboffset);
+};
+
+#endif