1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
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
|
