ParM  parm
A molecular dynamics library
box.hpp
Go to the documentation of this file.
1 #include "vecrand.hpp"
2 
3 #ifndef BOX_H
4 #define BOX_H
5 
6 #include <vector>
7 #include <boost/shared_ptr.hpp>
8 
9 #define sptr boost::shared_ptr
10 
11 using namespace boost; // required for SWIG for some reason
12 
13 typedef const unsigned int cuint;
14 
15 inline bool to_buffer(vector<Vec*> arr, double* buffer, size_t sizet) {
16  if(sizet < NDIM * arr.size()){return false;};
17  for(uint i=0; i < arr.size(); i++)
18  for(uint j=0; j < NDIM; j++)
19  {
20  buffer[i*NDIM + j] = (double) ((*arr[i])[j]);
21  }
22  return true;
23 };
24 
25 class AtomGroup;
26 
50 class Box {
51  public:
53 
56  virtual Vec diff(Vec r1, Vec r2)=0;
58  virtual flt V()=0;
59  virtual ~Box(){};
60 };
61 
62 /***********************************************************************
63  * Boxes
64  */
65 
67 
70 #ifdef VEC3D
71 inline Vec vec_mod(Vec r1, Vec r2){
72  return Vec(remainder(r1[0], r2[0]), remainder(r1[1], r2[1]), remainder(r1[2], r2[2]));
73 };
74 #endif
75 #ifdef VEC2D
76 inline Vec vec_mod(Vec r1, Vec r2){
77  return Vec(remainder(r1[0], r2[0]), remainder(r1[1], r2[1]));
78 };
79 #endif
80 
82 
86 class InfiniteBox : public Box {
87  public:
89  Vec diff(Vec r1, Vec r2){return r1-r2;};
91  flt V(){return NAN;};
92 };
93 
95 
98 class OriginBox : public Box {
99  protected:
101  public:
102  OriginBox(Vec size) : boxsize(size){};
103  Vec diff(Vec r1, Vec r2){
104  return vec_mod((r1-r2), boxsize);
105  };
106  virtual Vec diff(Vec r1, Vec r2, boost::array<int,NDIM> boxes){
107  Vec dr = r1 - r2;
108  for(uint i=0; i<NDIM; i++) dr[i] -= boxsize[i] * boxes[i];
109  return dr;
110  };
112  virtual boost::array<int,NDIM> box_round(Vec r1, Vec r2){
113  boost::array<int,NDIM> boxes;
114  Vec dr = r1 - r2;
115  for(uint i=0; i<NDIM; i++) boxes[i] = (int) round(dr[i] / boxsize[i]);
116  return boxes;
117  };
118  #ifdef VEC3D
119  OriginBox(flt L) : boxsize(L,L,L){};
120  flt V(){return boxsize[0] * boxsize[1] * boxsize[2];};
121  flt L(){return (boxsize[0] + boxsize[1] + boxsize[2])/3.0;};
122  #endif
123  #ifdef VEC2D
124  OriginBox(flt L) : boxsize(L,L){};
125  flt V(){return boxsize[0] * boxsize[1];};
126  flt L(){return (boxsize[0] + boxsize[1])/2.0;};
127  #endif
128  flt resize(flt factor);
131  flt resize(flt factor, AtomGroup &atoms);
133  flt resize_to(Vec newsize);
135  flt resize_to(Vec newsize, AtomGroup &atoms);
137  flt resize_to_V(flt newV);
139  flt resize_to_V(flt newV, AtomGroup &atoms);
140 
142  flt resize_to_L(flt newL);
144  flt resize_to_L(flt newL, AtomGroup &atoms);
147  Vec v = rand_vec_boxed();
148  for(uint i=0; i<NDIM; i++){
149  v[i] *= boxsize[i];
150  }
151  return diff(v, Vec::Zero());
152  };
153  Vec box_shape(){return boxsize;};
154 
156  void pure_shear_to(flt epsilon);
158  void pure_shear_to(flt epsilon, AtomGroup &atoms);
159 };
160 
162 
165 class LeesEdwardsBox : public OriginBox {
166  protected:
168  public:
169  LeesEdwardsBox(Vec size, flt gamma=0.0) : OriginBox(size), gamma(gamma){};
170  LeesEdwardsBox(flt L, flt gamma=0.0) : OriginBox(L), gamma(gamma){};
171  Vec diff(Vec r1, Vec r2);
172  virtual Vec diff(Vec r1, Vec r2, boost::array<int,NDIM> boxes);
173  virtual boost::array<int,NDIM> box_round(Vec r1, Vec r2);
174 
176  flt get_gamma(){return gamma;};
177 
179  void shear(flt dgamma);
181  void shear(flt dgamma, AtomGroup &atoms);
183  void shear_to(flt gamma);
185  void shear_to(flt gamma, AtomGroup &atoms);
187  v[0] -= gamma * v[1];
188  return v;
189  }
190 
192  v[0] += gamma * v[1];
193  return v;
194  }
195 };
196 
198 
208 class SCBox : public Box {
209  protected:
210  flt L, R;
211  public:
212  SCBox(flt L, flt R);
213  Vec diff(Vec r1, Vec r2){return r1-r2;};
214  flt V();
215  Vec dist(Vec r1);
216  Vec edge_dist(Vec r1);
217  bool inside(Vec r1, flt buffer=0.0);
218  Vec rand_loc(flt min_dist_to_wall=0.0);
219  flt length(){return L;};
220  flt radius(){return R;};
221 };
222 
225 
229 struct Atom {
232 
235 
238 
241 
244 };
245 
247 class AtomRef {
248  private:
249  Atom *ptr;
250  public:
251  inline AtomRef() : ptr(NULL){};
252  inline AtomRef(Atom *a) : ptr(a){};
253  inline Atom& operator *() const {return *ptr;};
254  inline Atom *operator->() const{ return ptr;}
255  inline bool operator==(const AtomRef &other) const {return other.ptr == ptr;};
256  inline bool operator==(const Atom* other) const {return other == ptr;};
257  inline bool operator!=(const AtomRef &other) const {return other.ptr != ptr;};
258  inline bool operator<(const AtomRef &other) const {return ptr < other.ptr;};
259  inline bool operator<=(const AtomRef &other) const {return ptr <= other.ptr;};
260  inline bool operator>=(const AtomRef &other) const {return ptr >= other.ptr;};
261  inline bool operator>(const AtomRef &other) const {return ptr > other.ptr;};
262  inline bool is_null(){return ptr==NULL;};
263 };
264 
266 
270 class AtomID : public AtomRef {
271  private:
272  uint num; // note that these are generally only in reference to
273  // a specific AtomGroup
274  public:
276  inline AtomID() : AtomRef(), num(UINT_MAX){};
278  inline AtomID(Atom *a, uint n) : AtomRef(a), num(n){};
279  inline uint n() const {return num;};
280 };
281 
282 class IDPair {
283  private:
284  AtomID id1, id2;
285  public:
286  IDPair() : id1(), id2(){};
287  IDPair(AtomID a, AtomID b) : id1(a), id2(b){};
288  inline AtomID first() const {return id1;};
289  inline AtomID last() const {return id2;};
290 };
291 
292 class AtomGroup;
293 
295 class AtomIter{
296  private:
297  uint i;
298  AtomGroup &g;
299  public:
300  AtomIter (AtomGroup& g, uint i): i(i), g(g){};
301  bool operator!=(const AtomIter& other) const{return i != other.i;};
302  Atom& operator* () const;
303  inline const AtomIter& operator++(){++i; return *this;};
304 };
305 
306 class AtomVec;
307 
309 
312 class AtomGroup {
313  public:
314  // access individual atoms
315  virtual AtomVec& vec()=0;
316  virtual Atom& operator[](cuint n)=0;
317  virtual Atom& operator[](cuint n) const=0;
318  virtual Atom& get(cuint n){return ((*this)[n]);};
319  virtual AtomID get_id(cuint n)=0;
321  virtual uint size() const=0;
323  virtual AtomIter begin(){return AtomIter(*this, 0);};
324  virtual AtomIter end(){return AtomIter(*this, (uint) size());};
325 
327  Vec com() const;
329  Vec com_force() const;
331  Vec com_velocity() const;
332 
334  flt mass() const;
336 
340  flt kinetic_energy(const Vec originvelocity=Vec::Zero()) const;
342  Vec momentum() const;
344  flt gyradius() const;
345  #ifdef VEC3D
346  Vec torque(const Vec loc) const;
348  Vec torque() const{return torque(com());};
350  flt moment_about(const Vec axis, const Vec loc) const;
351  flt moment_about(const Vec axis) const{return moment_about(axis, com());};
353  Vec angular_momentum(const Vec loc) const;
354  Vec angular_momentum() const {return angular_momentum(com());};
356  Matrix moment(const Vec loc) const;
357  Matrix moment() const{return moment(com());};
359  Vec omega(const Vec loc) const;
360  Vec omega() const {return omega(com());};
362  void add_omega(Vec w, Vec origin);
363  void add_omega(Vec w){return add_omega(w, com());};
365  inline void reset_L(){
366  Vec c = com(), w = omega(c);
367  if (w.squaredNorm() == 0) return;
368  add_omega(-w, c);
369  }
370  #elif defined VEC2D
371  flt torque(const Vec loc) const;
374  flt moment(const Vec loc) const;
375  flt moment() const{return moment(com());};
377  flt angular_momentum(const Vec loc) const;
378  flt angular_momentum() const {return angular_momentum(com());};
380  flt omega(const Vec loc) const{return angular_momentum(loc) / moment(loc);};
381  flt omega() const{return omega(com());};
383  void add_omega(flt w, Vec origin);
384  void add_omega(flt w){add_omega(w, com());}
386  inline void reset_L(){
387  Vec c = com();
388  flt w = omega(c);
389  if (w == 0) return;
390  add_omega(-w, c);
391  }
392  #endif
393 
395  void add_velocity(Vec v);
397  void reset_com_velocity(){add_velocity(-com_velocity());};
399  void randomize_velocities(flt T);
400 
402  void reset_forces();
403  virtual ~AtomGroup(){};
404 };
405 
406 inline Atom& AtomIter::operator*() const{return g[i];};
407 
412 class AtomVec : public virtual AtomGroup {
413  private:
414  Atom* atoms;
415  uint sz;
416  void zero() {
417  for(uint i=0; i < sz; i++){
418  atoms[i].x = Vec::Zero();
419  atoms[i].v = Vec::Zero();
420  atoms[i].f = Vec::Zero();
421  atoms[i].a = Vec::Zero();
422  }
423  }
424  public:
425  AtomVec(vector<double> masses) : sz((uint) masses.size()){
426  atoms = new Atom[sz];
427  for(uint i=0; i < sz; i++) atoms[i].m = masses[i];
428  zero();
429  };
430  AtomVec(uint N, flt mass) : sz(N){
431  atoms = new Atom[sz];
432  for(uint i=0; i < sz; i++) atoms[i].m = mass;
433  zero();
434  };
435  AtomVec(AtomVec& other) : sz(other.size()){
436  atoms = new Atom[sz];
437  for(uint i=0; i < sz; i++) atoms[i] = other.atoms[i];
438  };
439  AtomVec& vec(){return *this;};
440  inline Atom& operator[](cuint n){return atoms[n];};
441  inline Atom& operator[](cuint n) const {return atoms[n];};
442  inline AtomID get_id(cuint n) {
443  if (n > sz) return AtomID(); return AtomID(atoms + n,n);};
444  inline uint size() const {return sz;};
445 
446  ~AtomVec(){ delete [] atoms;};
447 };
448 
454 class SubGroup : public AtomGroup {
455  protected:
456  sptr<AtomVec> atoms;
457  vector<AtomID> ids;
458  public:
459  SubGroup(sptr<AtomVec> atoms) : atoms(atoms){};
460  AtomVec &vec(){return *atoms;};
461  inline Atom& operator[](cuint n){return *ids[n];};
462  inline Atom& operator[](cuint n) const{return *ids[n];};
463  inline Atom& get(cuint n){return *ids[n];};
464  inline void add(AtomID a){
465  std::vector<AtomID>::iterator it = std::find(ids.begin(), ids.end(), a);
466  if(it != ids.end())
467  throw std::invalid_argument("Cannot add AtomID to SubGroup: it already exists.");
468  return ids.push_back(a);
469  };
470  inline AtomID get_id(cuint n) {return ids[n];};
471  inline uint size() const {return (uint) ids.size();};
472 };
473 
474 #endif
Vec torque() const
Definition: box.hpp:348
OriginBox(flt L)
Definition: box.hpp:119
virtual AtomIter begin()
For use in a for loop.
Definition: box.hpp:323
sptr< AtomVec > atoms
Definition: box.hpp:456
Vec vec()
A one-liner for creating a Vec object, occasionally useful from within Python.
Definition: vecrand.hpp:72
uint n() const
Definition: box.hpp:279
Vec rand_loc()
Get a random point in the box.
Definition: box.hpp:146
A spheocylinder box, also known as a capsule.
Definition: box.hpp:208
The basic class for representing each particle.
Definition: box.hpp:229
SubGroup(sptr< AtomVec > atoms)
Definition: box.hpp:459
bool operator<=(const AtomRef &other) const
Definition: box.hpp:259
Atom * operator->() const
Definition: box.hpp:254
Vec affine(Vec v)
Definition: box.hpp:191
bool operator<(const AtomRef &other) const
Definition: box.hpp:258
Atom & operator[](cuint n) const
Definition: box.hpp:462
Vec non_affine(Vec v)
Definition: box.hpp:186
AtomVec(vector< double > masses)
Definition: box.hpp:425
~AtomVec()
Definition: box.hpp:446
flt m
mass
Definition: box.hpp:243
a group of atoms, such as all of them (AtomVec), or a smaller group such as a molecule, sidebranch, etc.
Definition: box.hpp:312
Atom & operator[](cuint n)
Definition: box.hpp:461
IDPair(AtomID a, AtomID b)
Definition: box.hpp:287
The main class for representing particles.
Definition: box.hpp:412
unsigned int uint
Definition: vec.hpp:20
Vec x
location.
Definition: box.hpp:231
OriginBox(Vec size)
Definition: box.hpp:102
virtual AtomIter end()
Definition: box.hpp:324
double flt
The basic floating point type used in the simulations.
Definition: vecrand.hpp:45
void reset_L()
Reset angular momentum to 0.
Definition: box.hpp:365
Vec v
velocity
Definition: box.hpp:234
Vec angular_momentum() const
Definition: box.hpp:354
flt length()
Definition: box.hpp:219
const AtomIter & operator++()
Definition: box.hpp:303
An infinite Box, for use with, e.g., sticky conglomerations or proteins.
Definition: box.hpp:86
A rectilinear Box, with periodic boundary conditions.
Definition: box.hpp:98
bool operator==(const Atom *other) const
Definition: box.hpp:256
Atom & operator[](cuint n)
Definition: box.hpp:440
For iterating through an AtomGroup.
Definition: box.hpp:295
AtomRef()
Definition: box.hpp:251
AtomIter(AtomGroup &g, uint i)
Definition: box.hpp:300
void add(AtomID a)
Definition: box.hpp:464
Vec diff(Vec r1, Vec r2)
Distance between two points, given boundary conditions.
Definition: box.hpp:103
Lees-Edwards boundary conditions, with shear in the x-direction, relative to y.
Definition: box.hpp:165
flt V()
Returns NaN.
Definition: box.hpp:91
A pointer to an Atom, that also knows its own index in an AtomVec.
Definition: box.hpp:270
Vec a
acceleration
Definition: box.hpp:237
Vec diff(Vec r1, Vec r2)
Simply r1-r2.
Definition: box.hpp:89
virtual ~AtomGroup()
Definition: box.hpp:403
flt V()
Volume. Can return NaN.
Definition: box.hpp:120
uint size() const
Number of atoms in the group.
Definition: box.hpp:444
flt gamma
Definition: box.hpp:167
LeesEdwardsBox(Vec size, flt gamma=0.0)
Definition: box.hpp:169
void add_omega(Vec w)
Definition: box.hpp:363
AtomVec & vec()
Definition: box.hpp:460
AtomVec(AtomVec &other)
Definition: box.hpp:435
Atom & operator[](cuint n) const
Definition: box.hpp:441
vector< AtomID > ids
Definition: box.hpp:457
virtual Vec diff(Vec r1, Vec r2, boost::array< int, NDIM > boxes)
Definition: box.hpp:106
Vec f
forces
Definition: box.hpp:240
virtual boost::array< int, NDIM > box_round(Vec r1, Vec r2)
The div function to go with diff.
Definition: box.hpp:112
AtomID()
Constructor not recommended: uses UINT_MAX for number, i.e. unknown.
Definition: box.hpp:276
virtual ~Box()
Definition: box.hpp:59
bool operator!=(const AtomIter &other) const
Definition: box.hpp:301
AtomVec(uint N, flt mass)
Definition: box.hpp:430
flt L()
Definition: box.hpp:121
bool to_buffer(vector< Vec * > arr, double *buffer, size_t sizet)
Definition: box.hpp:15
AtomID first() const
Definition: box.hpp:288
AtomID get_id(cuint n)
Definition: box.hpp:442
bool operator==(const AtomRef &other) const
Definition: box.hpp:255
AtomID get_id(cuint n)
Definition: box.hpp:470
AtomID(Atom *a, uint n)
Recommended constructor.
Definition: box.hpp:278
AtomID last() const
Definition: box.hpp:289
A 3x3 matrix, with methods for adding, subtracting, dot product, etc.
Definition: vec.hpp:360
void reset_com_velocity()
Subtracts the center of mass velocity from all atoms.
Definition: box.hpp:397
Matrix moment() const
Definition: box.hpp:357
Vec vec_mod(Vec r1, Vec r2)
The modulus function for Vec.
Definition: box.hpp:71
AtomRef(Atom *a)
Definition: box.hpp:252
Atom & operator*() const
Definition: box.hpp:406
Vec box_shape()
Definition: box.hpp:153
Vec boxsize
Definition: box.hpp:100
IDPair()
Definition: box.hpp:286
bool operator>(const AtomRef &other) const
Definition: box.hpp:261
bool operator!=(const AtomRef &other) const
Definition: box.hpp:257
flt get_gamma()
The current shear amount.
Definition: box.hpp:176
bool operator>=(const AtomRef &other) const
Definition: box.hpp:260
#define NDIM
Definition: vecrand.hpp:12
Definition: box.hpp:282
uint size() const
Number of atoms in the group.
Definition: box.hpp:471
The virtual interface for the shape of the space and its boundaries.
Definition: box.hpp:50
AtomVec & vec()
Definition: box.hpp:439
Eigen::Matrix< flt, NDIM, 1 > Vec
The basic physics vector.
Definition: vecrand.hpp:53
Definition: box.hpp:454
Vec rand_vec_boxed()
Generate a random vector inside a box with sides of length 1.
Definition: vecrand.cpp:19
Vec omega() const
Definition: box.hpp:360
flt radius()
Definition: box.hpp:220
LeesEdwardsBox(flt L, flt gamma=0.0)
Definition: box.hpp:170
flt moment_about(const Vec axis) const
Definition: box.hpp:351
bool is_null()
Definition: box.hpp:262
A pointer to an Atom.
Definition: box.hpp:247
flt R
Definition: box.hpp:210
Vec diff(Vec r1, Vec r2)
Distance between two points, given boundary conditions.
Definition: box.hpp:213
const unsigned int cuint
Definition: box.hpp:13