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

1 files changed, 64 insertions, 0 deletions

diff --git a/src/spheres_poly/spheres.cpp b/src/spheres_poly/spheres.cpp
new file mode 100644
index 0000000..9b14b01
--- /dev/null
+++ b/src/spheres_poly/spheres.cpp
@@ -0,0 +1,64 @@
+//===========================================================
+//===========================================================
+//===========================================================
+//
+//  Molecular dynamics simulation of hardspheres
+//
+//===========================================================
+//===========================================================
+//===========================================================
+
+#include <iostream>
+#include <math.h>
+#include <fstream>
+#include <vector>
+#include <time.h>
+#include <string.h>
+#include <stdlib.h>
+#include <cstdlib>
+
+#include "box.h"
+#include "sphere.h"
+#include "event.h"
+#include "heap.h"
+
+
+double *spherespp(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate)
+{
+	int eventspercycle = 20;              // # events per sphere per cycle
+	double initialpf = 0.01;              // initial packing fraction
+	double temp = 0.2;                    // initial temp., use 0 for zero velocities
+	double growthrate = 0.001;            // growth rate
+	double massratio = 1.;               // ratio of sphere masses
+	int hardwallBC = 0;                   // 0 for periodic, 1 for hard wall BC
+
+  double d, r;   // initial diameter and radius of spheres
+
+  r = pow(initialpf*pow(SIZE, DIM)/(N*VOLUMESPHERE), 1.0/((double)(DIM)));
+
+  box b(N, r, growthrate, maxpf, bidispersityratio, 
+	bidispersityfraction, massratio, hardwallBC);
+
+	b.CreateSpheres(temp);
+  
+  
+  while ((b.collisionrate < maxcollisionrate) && (b.pf < maxpf) && (b.pressure < maxpressure)) 
+    {
+      b.Process(eventspercycle*N);
+      b.Synchronize(true);
+    }
+  
+	double *data = (double *)malloc(DIM * N * sizeof(double));
+
+  b.WriteConfiguration(data);
+  
+  return data;
+}
+
+extern "C" {
+	double *spheres(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate) {
+		return spherespp(N, bidispersityratio, bidispersityfraction, maxpf, maxpressure, maxcollisionrate);
+	}
+}
+
+