ParM  parm
A molecular dynamics library
interaction.hpp
Go to the documentation of this file.
1 #include "trackers.hpp"
2 
3 #ifndef INTERACTION_H
4 #define INTERACTION_H
5 
6 #define sptr boost::shared_ptr
7 using namespace boost; // required for SWIG for some reason
8 
9 /***********************************************************************
10  * L-J Note:
11  * The convention here is ε(σ⁶/r⁶ - 1)², which has its minimum at r = σ.
12  * This is what to use with Alice's parameters.
13  */
14 
15 /*
16  * New plan of attack:
17  *
18  * Interaction<N> - blindly takes N vectors,
19  * and calculates energy or forces.
20  *
21  * AtomGroup(N) (or (N,M) with array?)-
22  * has a list of atoms, can tell mass, vel, loc.
23  * keeps forces?
24  * possibly has M groups of groups?
25  * iterator over atoms necessary
26  * perhaps class iterator over all of them?
27  * perhaps class iterator over all atoms?
28  * has "uint size()" function, says number of atoms
29  *
30  * InteractGroup - has a single Interaction<N>, and pointers to however
31  * many AtomGroups it needs, and returns energy and/or sets forces.
32  * one example: neighbors; yields pairs of i,i+1.
33  * second example: neighbor list, keeps track of nearby atoms.
34  * perhaps class iterator over all of them?
35  *
36  * perhaps split: interactiongroup has a group, and an Interaction, and
37  * applies Interaction to each pair in each group.
38  * Not incompatible.
39  *
40  * integrator - keeps a list of Interaction groups, can get energy.
41  * Also goes through list to integrate over time.
42  *
43  * statistic - keeps a list of Interaction groups, derives a statistic
44  * from them.
45  *
46  * statkeeper - keeps a list of statistics, can write to file?
47  *
48  * constants - keeps track of constants; in map, or as properties?
49  * perhaps a struct.
50 */
51 
52 //typedef double flt; defined in vecrand.hpp
53 //typedef unsigned int uint;
58 class Interaction {
59  public:
63  virtual flt energy(Box &box)=0;
64  virtual void set_forces(Box &box)=0;
70  std::string s = std::string("set_forces_get_pressure not defined for class ");
71  s.append(typeid(*this).name());
72  throw std::runtime_error(s);
73  };
84  virtual flt pressure(Box &box)=0;
85 
93  virtual Matrix stress(Box &box){
94  std::string s = std::string("stress not defined for class ");
95  s.append(typeid(*this).name());
96  throw std::runtime_error(s);
97  }
98  virtual ~Interaction(){};
99 };
100 
101 /***********************************************************************
102  * Interaction Basics
103  */
105  public:
106  virtual flt energy(const Vec diff)=0;
107  virtual Vec forces(const Vec diff)=0;
108  virtual ~InteractPair(){};
109 };
110 
111 static const flt LJr0 = pow(2.0, 1.0/6.0);
112 static const flt LJr0sq = pow(2.0, 1.0/3.0);
113 
114 class LJRepulsive {
115  protected:
118  public:
119  LJRepulsive(const flt epsilon, const flt sigma):
120  epsilon(epsilon), sigma(sigma){};
121  inline static flt energy(const Vec diff, const flt eps, const flt sig){
122  flt rsq = diff.squaredNorm()/(sig*sig);
123  if(rsq > 1) return 0;
124  flt rsix = rsq*rsq*rsq;
125  //~ return eps*(4*(1/(rsix*rsix) - 1/rsix) + 1);
126  flt mid = (1-1/rsix);
127  return eps*(mid*mid);
128  };
129  inline flt energy(const Vec& diff){return energy(diff, epsilon, sigma);};
130  inline static Vec forces(const Vec diff, const flt eps, const flt sig){
131  flt dsq = diff.squaredNorm();
132  flt rsq = dsq/(sig*sig);
133  if(rsq > 1) return Vec::Zero();
134  flt rsix = rsq*rsq*rsq; //r^6 / sigma^6
135  //~ flt fmagTimesR = eps*(4*(12/(rsix*rsix) - 6/rsix));
136  flt fmagTimesR = 12*eps/rsix*(1/rsix - 1);
137  //~ cout << "Repulsing " << diff << "with force"
138  //~ << diff * (fmagTimesR / dsq) << '\n';
139  //~ cout << "mag: " << diff.norm() << " sig:" << sig << " eps:" << eps
140  //~ << "F: " << (diff * (fmagTimesR / dsq)).norm() << '\n';
141  //~ cout << "E: " << energy(diff, sig, eps) << " LJr0: " << LJr0 << ',' << LJr0sq << "\n";
142  return diff * (fmagTimesR / dsq);
143  };
144  inline Vec forces(const Vec& diff){return forces(diff, epsilon, sigma);};
145 };
146 
147 class LJAttract {
148  protected:
151  public:
152  LJAttract(const flt epsilon, const flt sigma):
153  epsilon(epsilon), sigma(sigma){};
154  inline static flt energy(const Vec diff, const flt eps, const flt sig){
155  flt rsq = diff.squaredNorm()/(sig*sig);
156  if(rsq < 1) return -eps;
157  flt rsix = rsq*rsq*rsq;
158  //~ return eps*(4*(1/(rsix*rsix) - 1/rsix) + 1);
159  flt mid = (1-1/rsix);
160  return eps*(mid*mid-1);
161  };
162  inline static flt energy(const flt rsig){
163  if(rsig < 1) return -1;
164  flt rsq = rsig * rsig;
165  flt rsix = rsq*rsq*rsq;
166  flt mid = (1-1/rsix);
167  return mid*mid-1;
168  // Minimum of -1 occurs at r = 1
169  };
170  inline flt energy(const Vec& diff){return energy(diff, epsilon, sigma);};
171  inline static Vec forces(const Vec diff, const flt eps, const flt sig){
172  flt dsq = diff.squaredNorm();
173  flt rsq = dsq/(sig*sig);
174  if(rsq < 1) return Vec::Zero();
175  flt rsix = rsq*rsq*rsq; //r^6 / sigma^6
176  flt fmagTimesR = 12*eps/rsix*(1/rsix - 1);
177  return diff * (fmagTimesR / dsq);
178  };
179  inline static flt forces(const flt rsig){
180  if(rsig < 1) return 0;
181  flt rsq = rsig*rsig;
182  flt rsix = rsq*rsq*rsq; //r^6 / sigma^6
183  flt fmagTimesR = 12/rsix*(1/rsix - 1);
184  return fmagTimesR / rsig;
185  };
186  inline Vec forces(const Vec& diff){return forces(diff, epsilon, sigma);};
187 };
188 
190  // Purely attractive
191  protected:
194  flt cut_distance, cut_energy;
195  public:
196  inline LJAttractCut(const flt epsilon, const flt sigma, const flt cutsig):
197  epsilon(epsilon), sigma(sigma), cut_distance(cutsig),
198  cut_energy(LJAttract::energy(cut_distance) * epsilon){};
199  inline static flt energy(const Vec diff, const flt eps,
200  const flt sig, const flt cutsig){
201  if(eps == 0) return 0;
202  if(diff.squaredNorm() > (cutsig*cutsig*sig*sig)) return 0;
203  return (LJAttract::energy(diff, eps, sig) -
204  eps*LJAttract::energy(cutsig));
205  };
206  inline flt energy(const Vec& diff){
207  if(epsilon == 0) return 0;
208  if(diff.squaredNorm() > (cut_distance*cut_distance*sigma*sigma)) return 0;
209  return LJAttract::energy(diff, epsilon, sigma) - cut_energy;
210  };
211  inline static Vec forces(const Vec diff, const flt eps,
212  const flt sig, const flt cutsig){
213  if(eps == 0) return Vec::Zero();
214  flt dsq = diff.squaredNorm();
215  flt rsq = dsq/(sig*sig);
216  if(rsq < 1 or rsq > (cutsig*cutsig)) return Vec::Zero();
217  flt rsix = rsq*rsq*rsq; //r^6 / sigma^6
218  flt fmagTimesR = 12*eps/rsix*(1/rsix - 1);
219  return diff * (fmagTimesR / dsq);
220  };
221  inline Vec forces(const Vec& diff){
222  return forces(diff, epsilon, sigma, cut_distance);};
223 };
224 
226  // uses ε*((1-1/r⁶)² - 1)
227  // sigma is thus at the minimum
228  // cut_distance is unitless; that times sigma is the distance we cut the potential
229  protected:
232  flt cut_distance, cut_energy;
233  public:
234  LennardJonesCut(const flt epsilon, const flt sigma, const flt cutsig):
235  epsilon(epsilon), sigma(sigma), cut_distance(cutsig){
236  flt rsix = pow(cut_distance, 6);
237  flt mid = (1-1/rsix);
238  cut_energy = epsilon*(mid*mid-1);
239  };
240  inline flt energy(const Vec& diff){
241  flt rsq = diff.squaredNorm()/(sigma*sigma);
242  if(rsq > (cut_distance*cut_distance)) return 0;
243  flt rsix = rsq*rsq*rsq;
244  flt mid = (1-1/rsix);
245  return epsilon*(mid*mid-1) - cut_energy;
246  };
247  inline static Vec forces(const Vec diff, const flt eps, const flt sig, const flt cutsig){
248  flt dsq = diff.squaredNorm();
249  flt rsq = dsq/(sig*sig);
250  if(rsq > (cutsig*cutsig)) return Vec::Zero();
251  flt rsix = rsq*rsq*rsq; //r^6 / sigma^6
252  flt fmagTimesR = 12*eps/rsix*(1/rsix - 1);
253  return diff * (fmagTimesR / dsq);
254  };
255  inline static flt forces(const flt rsig, const flt cutsig){
256  if(rsig > cutsig) return 0;
257  flt rsq = rsig*rsig;
258  flt rsix = rsq*rsq*rsq; //r^6 / sigma^6
259  flt fmagTimesR = 12/rsix*(1/rsix - 1);
260  return fmagTimesR / rsig;
261  };
262  inline Vec forces(const Vec& diff){return forces(diff, epsilon, sigma, cut_distance);};
263 };
264 
265 class Spring : public InteractPair {
266  protected:
269  public:
270  Spring(const flt k, const flt x0) : springk(k),x0(x0){};
271  flt energy(const Vec diff);
272  Vec forces(const Vec diff);
273  ~Spring(){};
274 };
275 
276 class BondAngle {
277  // two vectors for E and F are from the central one to the other 2
278  protected:
281  bool usecos;
282 
283  public:
284  BondAngle(const flt k, const flt theta, const bool cosine=false)
285  :springk(k), theta0(theta), usecos(cosine){};
286  inline static flt get_angle(const Vec& r1, const Vec& r2){
287  flt costheta = r1.dot(r2) / r1.norm() / r2.norm();
288  if(costheta > 1) costheta = 1;
289  else if(costheta < -1) costheta = -1;
290  return acos(costheta);
291  };
292  flt energy(const Vec& diff1, const Vec& diff2);
293  array<Vec,3> forces(const Vec& diff1, const Vec& diff2);
295 };
296 
297 #ifdef VEC3D
299  array<Vec,4> derivs;
301 };
302 
303 class Dihedral {
304  // vectors are from 1 to 2, 2 to 3, 3 to 4
305  protected:
306  vector<flt> cos_coefficients;
307  vector<flt> sincoeffs;
308  bool usepow;
309  flt dU_dcostheta_cos(const flt costheta) const;
310  public:
311  Dihedral(const vector<flt> cosvals,
312  const vector<flt> sinvals = vector<flt>(),
313  bool usepow = true);
314  static flt get_cos(const Vec& diff1, const Vec& diff2, const Vec& diff3);
315  static flt get_angle(const Vec& diff1, const Vec& diff2, const Vec& diff3);
316 
318  static DihedralDerivs dr_dcostheta(const Vec& diff1, const Vec& diff2, const Vec& diff3);
319 
320  flt dU_dcostheta(const flt theta) const;
321 
322  inline flt energy(const Vec& diff1, const Vec& diff2, const Vec& diff3) const {
323  return energy(get_angle(diff1, diff2, diff3));
324  };
325  flt energy(flt ang) const;
326  array<Vec,4> forces(const Vec& diff1, const Vec& diff2, const Vec& diff3) const;
327 };
328 #endif
329 
331  protected:
333  flt q1, q2;
336  public:
337  ElectricScreened(const flt screenLength, const flt q1,
338  const flt q2, const flt cutoff);
339  inline flt energy(const Vec r){return energy(r.norm(),q1*q2,screen, 0) - cutoffE;};
340  static flt energy(const flt r, const flt qaqb, const flt screen, const flt cutoff=0);
341  inline Vec forces(const Vec r){return forces(r,q1*q2,screen, cutoff);};
342  static Vec forces(const Vec r, const flt qaqb, const flt screen, const flt cutoff=0);
343 };
344 
345 //~ class InteractGroup : public Interaction {
346  //~ protected:
347  //~ vector<Interaction*> inters;
348  //~ public:
349  //~ InteractGroup(vector<Interaction*> inters=vector<Interaction*>())
350  //~ : inters(inters){};
351  //~ void add(Interaction* a){inters.push_back(a);};
352  //~ uint size() const{ return inters.size();};
353  //~ flt energy(Box *box);
354  //~ void set_forces(Box *box);
355 //~ };
356 
360  FixedForceAtom(Vec F, AtomID a) : F(F), a(a) {};
361  flt energy(Box &box){return -F.dot(a->x);};
362  void set_force(Box &box){a->f += F;};
363 };
364 
365 class FixedForce : public Interaction {
366  protected:
367  vector<FixedForceAtom> atoms;
368  public:
369  FixedForce(vector<FixedForceAtom> atoms = vector<FixedForceAtom>()) : atoms(atoms){};
370  void add(FixedForceAtom a){atoms.push_back(a);};
371  void add(Vec F, AtomID a){add(FixedForceAtom(F,a));};
372  #ifdef VEC3D
373  void add(flt x, flt y, flt z, AtomID a){add(FixedForceAtom(Vec(x,y,z),a));};
374  #elif defined VEC2D
375  void add(flt x, flt y, AtomID a){add(FixedForceAtom(Vec(x,y),a));};
376  #endif
377  uint size() const{ return (uint) atoms.size();};
378  flt energy(Box &box){
379  flt E=0;
380  for(vector<FixedForceAtom>::iterator it = atoms.begin(); it < atoms.end(); ++it)
381  E += it->energy(box);
382  return E;
383  };
384  void set_forces(Box &box){
385  for(vector<FixedForceAtom>::iterator it = atoms.begin(); it < atoms.end(); ++it)
386  it->set_force(box);
387  };
388  flt pressure(Box &box){return NAN;};
389 };
390 
391 
392 struct FixedForceRegionAtom : public AtomID {
393  Vec direction; // will be normalized
394  vector<flt> boundaries; // should be 1 smaller than Fs
395  // each boundary is at (direction * b),
396  // where b is in boundaries
397  vector<flt> Fs;
398  FixedForceRegionAtom(AtomID a, Vec direction, vector<flt> boundaries, vector<flt> Fs);
399  flt energy(Box &box);
400  void set_force(Box &box);
401 };
402 
404 
405  protected:
406  vector<FixedForceRegionAtom> atoms;
407  public:
408  FixedForceRegion(vector<FixedForceRegionAtom> atoms = vector<FixedForceRegionAtom>())
409  : atoms(atoms){};
410 
411  void add(FixedForceRegionAtom a){atoms.push_back(a);};
412  void add(AtomID a, Vec dir, vector<flt> bound, vector<flt> F){
413  add(FixedForceRegionAtom(a, dir, bound, F));};
414  uint size() const{ return (uint) atoms.size();};
415 
416  flt energy(Box &box){
417  flt E=0;
418  for(vector<FixedForceRegionAtom>::iterator it = atoms.begin(); it < atoms.end(); ++it)
419  E += it->energy(box);
420  return E;
421  };
422  void set_forces(Box &box){
423  for(vector<FixedForceRegionAtom>::iterator it = atoms.begin(); it < atoms.end(); ++it)
424  it->set_force(box);
425  };
426  flt pressure(Box &box){return NAN;};
427 };
428 
432  bool usecoord[3];
434  FixedSpringAtom(AtomID a, Vec loc, flt k, bool usex=true, bool usey=true, bool usez=true) :
435  loc(loc), k(k), a(a) {
436  usecoord[0] = usex;
437  usecoord[1] = usey;
438  usecoord[2] = usez;
439  };
440  flt energy(Box &box){
441  Vec diffx = a->x - loc;
442  for(uint i=0; i<2; ++i){
443  if(!usecoord[i]) diffx[i] = 0;
444  }
445  return k*diffx.squaredNorm()/2;
446  };
447  void set_force(Box &box){
448  Vec diffx = a->x - loc;
449  for(uint i=0; i<2; ++i){
450  if(!usecoord[i]) diffx[i] = 0;
451  }
452  a->f -= diffx * k;
453  };
454 };
455 
456 class FixedSpring : public Interaction{
457  protected:
458  vector<FixedSpringAtom> atoms;
459  public:
460  FixedSpring(vector<FixedSpringAtom> atoms = vector<FixedSpringAtom>()) : atoms(atoms){};
461  void add(FixedSpringAtom a){atoms.push_back(a);};
462  void add(AtomID a, Vec loc, flt k, bool usex=true, bool usey=true, bool usez=true){
463  add(FixedSpringAtom(a, loc, k, usex, usey, usez));};
464  //void add(Vec F, AtomID a){add(FixedForceAtom(F,a));};
465  //void add(flt x, flt y, flt z, AtomID a){add(FixedForceAtom(Vec(x,y,z),a));};
466  uint size() const{ return (uint) atoms.size();};
467  flt energy(Box &box){
468  flt E=0;
469  for(vector<FixedSpringAtom>::iterator it = atoms.begin(); it < atoms.end(); ++it)
470  E += it->energy(box);
471  return E;
472  };
473  void set_forces(Box &box){
474  for(vector<FixedSpringAtom>::iterator it = atoms.begin(); it < atoms.end(); ++it)
475  it->set_force(box);
476  };
477  flt pressure(Box &box){return NAN;};
478 };
479 
480 class COMSpring : public Interaction{
481  protected:
484  flt k, x0;
485  flt m1, m2;
486  public:
487  COMSpring(AtomGroup *g1, AtomGroup *g2, flt k, flt x0=0) :
488  g1(g1), g2(g2), k(k), x0(x0), m1(g1->mass()), m2(g2->mass()){};
489  flt energy(Box &box){
490  flt dx = (g1->com() - g2->com()).norm() - x0;
491  return k/2 * dx * dx;
492  };
493  void set_forces(Box &box){
494  Vec comvec = g1->com() - g2->com();
495  flt comdist = comvec.norm();
496  flt fmag = -k * (comdist - x0);
497  Vec a1 = comvec * (fmag / m1 / comdist);
498  for(uint i=0; i < g1->size(); ++i){
499  Atom &atm = g1->get(i);
500  atm.f += a1 * atm.m;
501  }
502 
503  Vec a2 = comvec * (-fmag / m2 / comdist);
504  for(uint i=0; i < g2->size(); ++i){
505  Atom &atm = g2->get(i);
506  atm.f += a2 * atm.m;
507  }
508  };
509  flt pressure(Box &box){
510  //~ return NAN;
511  // I think this is right, but I haven't checked it
512  Vec comvec = g1->com() - g2->com();
513  flt comdist = comvec.norm();
514  flt fmag = -k * (comdist - x0);
515 
516  Vec a12 = comvec * (fmag / m1 / m2 / comdist);
517 
518  flt P = 0;
519  for(uint i=0; i < g1->size(); ++i){
520  for(uint j=0; j < g2->size(); ++j){
521  Atom &atm1 = g1->get(i);
522  Atom &atm2 = g2->get(i);
523  Vec fij = a12 * (atm1.m * atm2.m);
524  Vec rij = box.diff(atm1.x, atm2.x);
525  P += fij.dot(rij);
526  }
527  }
528  return P;
529  };
530 };
531 
533 // Random Force
534 
536  FIXED, // Always the same magnitude
537  UNIFORM, // uniform probability distribution for magnitude, from 0 to force_mag
538  GAUSSIAN // Gaussian probability distribution for magnitude
539 };
540 
541 struct RandomForceAtom : public AtomRef {
542  public:
544  flt freq; // in general, how many timesteps on average between kicks
546  public:
547  RandomForceAtom(AtomID a, flt force_mag, flt freq, RandomForceType force_type=UNIFORM) :
548  AtomRef(a), force_mag(force_mag), freq(freq), force_type(force_type){};
549 };
550 
551 class RandomForce : public Interaction {
552  public:
553  vector<RandomForceAtom> group;
554 
555  public:
557  RandomForce(AtomGroup& agroup, flt force_mag, flt freq, RandomForceType force_type=UNIFORM){
558  for(uint i=0; i<agroup.size(); ++i){
559  group.push_back(RandomForceAtom(agroup.get_id(i), force_mag, freq, force_type));
560  }
561  };
562 
563  uint size() const{ return (uint) group.size();};
564  RandomForceAtom get(uint i) const{ return group[i];};
565 
568  bool add(RandomForceAtom a, bool replace=true);
569 
570  // No potential energy
571  flt energy(Box &box){return 0;};
572  void set_forces(Box &box);
573 
574  // no pressure from this Interaction
575  //TODO: maybe this should be average pressure
576  flt pressure(Box &box){return 0;};
577 
578  // no pressure from this Interaction
579  //TODO: maybe this should actually work; it could
580  flt set_forces_get_pressure(Box &box){set_forces(box); return 0;};
581 
582 };
583 
585 // Bonds
586 
588  BOXED, // use box.diff(r1, r2)
589  UNBOXED, // use r2 - r1
590  FIXEDBOX // use r2 - r1 - (original box separation)
591 };
592 
593 struct BondGrouping {
594  flt k, x0;
595  AtomID a1, a2;
597  array<int,NDIM> fixed_box;
598  BondGrouping(flt k, flt x0, AtomID a1, AtomID a2,
599  BondDiffType diff=UNBOXED, OriginBox *box=NULL);
600  Vec diff(Box &box) const;
601  int get_fixed(uint i){return fixed_box[i];};
602  inline bool same_atoms(BondGrouping &other){
603  return ((a1 == other.a1) and (a2 == other.a2)) or ((a1 == other.a2) and (a2 == other.a1));
604  };
605 };
606 
607 class BondPairs : public Interaction {
608  protected:
610  vector<BondGrouping> pairs;
611  //inline static Vec diff(Box &box, Vec r1, Vec r2){return r1-r2;};
612  //inline static Vec diff(Box &box, Vec r1, Vec r2){return box.diff(r1, r2);};
613  public:
614  BondPairs(vector<BondGrouping> pairs, bool zeropressure=true);
615  BondPairs(bool zeropressure=true);
619  bool add(BondGrouping b, bool replace=true);
620  inline bool add(flt k, flt x0, AtomID a1, AtomID a2, bool replace=true){
621  return add(BondGrouping(k,x0,a1,a2), replace);};
622  void add_forced(BondGrouping b){pairs.push_back(b);};
624  inline bool add(flt k, AtomID a1, AtomID a2, bool replace=true){
625  flt x0 = (a1->x - a2->x).norm();
626  return add(BondGrouping(k,x0,a1,a2), replace);
627  };
628 
629  uint size() const{ return (uint) pairs.size();};
630  BondGrouping get(uint i) const{ return pairs[i];};
631  flt mean_dists(Box &box) const;
632  flt std_dists(Box &box) const;
633  flt energy(Box &box);
634  void set_forces(Box &box);
635  flt pressure(Box &box);
636  flt set_forces_get_pressure(Box &box);
637 };
638 
640  flt k, x0;
641  AtomID a1, a2, a3;
642  AngleGrouping(flt k, flt x0, AtomID a1, AtomID a2, AtomID a3) :
643  k(k),x0(x0), a1(a1), a2(a2), a3(a3){};
644  inline bool same_atoms(AngleGrouping &other){
645  if(a2 != other.a2) return false;
646  return ((a1 == other.a1) and (a3 == other.a3)) or ((a1 == other.a3) and (a3 == other.a1));
647  };
648 };
649 
650 class AngleTriples : public Interaction {
651  protected:
652  vector<AngleGrouping> triples;
653  inline static Vec diff(Vec r1, Vec r2){return r1-r2;};
654  public:
655  AngleTriples(vector<AngleGrouping> triples = vector<AngleGrouping>());
659  bool add(AngleGrouping b, bool replace=true);
660  inline bool add(flt k, flt x0, AtomID a1, AtomID a2, AtomID a3, bool replace=true){
661  return add(AngleGrouping(k,x0,a1,a2,a3), replace);};
663  bool add(flt k, AtomID a1, AtomID a2, AtomID a3, bool replace=true);
664  void add_forced(AngleGrouping b){triples.push_back(b);};
665 
666  flt energy(Box &box);
667  inline flt pressure(Box &box){return 0;};
668  void set_forces(Box &box);
669  inline flt set_forces_get_pressure(Box &box){set_forces(box); return 0;};
670  uint size() const {return (uint) triples.size();};
671  flt mean_dists() const;
672  flt std_dists() const;
673 };
674 
675 #ifdef VEC3D
677  inline static Vec diff(Vec r1, Vec r2){return r1-r2;};
678  Dihedral dih;
679  AtomID a1, a2, a3, a4;
680  DihedralGrouping(vector<flt> cos_coefficients, vector<flt> sincoeffs,
681  AtomID a1, AtomID a2, AtomID a3, AtomID a4, bool usepow=true) :
682  dih(cos_coefficients, sincoeffs, usepow), a1(a1),
683  a2(a2), a3(a3), a4(a4){};
684 };
685 
686 class Dihedrals : public Interaction {
687  protected:
688  vector<DihedralGrouping> groups;
689  public:
690  Dihedrals(vector<DihedralGrouping> pairs = vector<DihedralGrouping>());
691  void add(DihedralGrouping b){groups.push_back(b);};
692  inline void add(vector<flt> nums, AtomID a1, AtomID a2, AtomID a3, AtomID a4){
693  add(DihedralGrouping(nums, vector<flt>(), a1,a2,a3,a4));};
694  inline void add(vector<flt> cos_coefficients, vector<flt> sincoeffs,
695  AtomID a1, AtomID a2, AtomID a3, AtomID a4, bool usepow=true){
696  add(DihedralGrouping(cos_coefficients, sincoeffs,a1,a2,a3,a4, usepow));};
698  inline void add(flt k, flt theta0, AtomID a1, AtomID a2, AtomID a3, AtomID a4){
699  vector<flt> cos_coefficients(2,k);
700  cos_coefficients[1] = -k*cos(theta0);
701  vector<flt> sincoeffs(2,0);
702  sincoeffs[1] = -k*sin(theta0);
703  add(DihedralGrouping(cos_coefficients, sincoeffs, a1,a2,a3,a4, true));
704  }
709  inline void add(flt k, AtomID a1, AtomID a2, AtomID a3, AtomID a4){
710  Vec r1 = DihedralGrouping::diff(a2->x, a1->x);
711  Vec r2 = DihedralGrouping::diff(a3->x, a2->x);
712  Vec r3 = DihedralGrouping::diff(a4->x, a3->x);
713  add(k, Dihedral::get_angle(r1, r2, r3), a1, a2, a3, a4);
714  };
715  uint size() const{ return uint(groups.size());};
716  flt mean_dists() const;
717  //~ flt std_dists() const;
718  flt energy(Box &box);
719  void set_forces(Box &box);
720  inline flt pressure(Box &box){return 0;};
721  inline flt set_forces_get_pressure(Box &box){set_forces(box); return 0;};
722 };
723 #endif
724 
725 
727 struct ForcePair {
728  Atom a1, a2;
730 };
731 
732 struct ForcePairX {
733  AtomID a1, a2;
736 };
737 
738 //
739 class FPairXFunct {
740  public:
741  virtual void run(ForcePairX*) = 0;
742  virtual ~FPairXFunct() = 0;
743 };
744 
745 inline FPairXFunct::~FPairXFunct() { } // defined even though it's pure virtual; it's faster this way; trust me
746 
748  public:
750  virtual void set_forces(Box &box, FPairXFunct*)=0;
751  virtual ~InteractionPairsX(){};
752 };
753 
754 struct Charged : public AtomID {
756  Charged() : AtomID(), q(0){};
757  Charged(flt q, AtomID a) : AtomID(a), q(q){};
758 };
759 
760 struct ChargePair {
762  AtomID atom1, atom2;
763  ChargePair(Charged a1, Charged a2) : q1q2(a1.q*a2.q){};
764 };
765 
766 struct EpsSigAtom : public AtomID {
767  flt epsilon, sigma;
769  EpsSigAtom(AtomID a, flt epsilon, flt sigma) : AtomID(a),
770  epsilon(epsilon), sigma(sigma){};
772  epsilon(other.epsilon), sigma(other.sigma){};
773  flt max_size(){return sigma;};
774 };
775 
784  flt epsilon, sigma;
785  AtomID atom1, atom2;
787  epsilon(sqrt(LJ1.epsilon * LJ2.epsilon)),
788  sigma((LJ1.sigma + LJ2.sigma) / 2),
789  atom1(LJ1), atom2(LJ2){};
790  inline flt energy(Box &box){
791  return LJRepulsive::energy(box.diff(atom1->x,atom2->x), epsilon, sigma);};
792  inline Vec forces(Box &box){
793  return LJRepulsive::forces(box.diff(atom1->x,atom2->x), epsilon, sigma);};
794 };
795 
797 // Purely attractive LJ, with ε = √(ε₁ ε₂) and σ = (σ₁+σ₂)/2
798 // cutoff at some sigcut
799 // Minimum at σ
800 struct EpsSigCutAtom : public EpsSigAtom {
801  flt sigcut; // sigma units
803  EpsSigCutAtom(AtomID a, flt epsilon, flt sigma, flt cut) :
804  EpsSigAtom(a, epsilon, sigma), sigcut(cut){};
806  sigcut(other.sigcut){};
807  flt max_size(){return sigma*sigcut;};
808 };
809 
810 
811 
813 // Purely attractive LJ, with ε = ε₁₂ and σ = σ₁₂ (both indexed)
814 // cutoff at some sigcut (and r < σ)
815 
816 struct IEpsISigCutAtom : public AtomID {
817  vector<flt> epsilons; // for finding epsilons
818  vector<flt> sigmas; // for finding epsilons
819  uint indx; // for finding this one in other atoms' epsilon lists
820  // note that for two HydroAtoms a1, a2, the epsilon for the pair
821  // is then either a1.epsilons[a2.indx] or a2.epsilons[a1.indx]; same for sigma
822  flt sigcut; // sigma units
824  IEpsISigCutAtom(AtomID a, vector<flt> epsilons, vector<flt> sigmas, uint indx, flt cut) :
825  AtomID(a), epsilons(epsilons), sigmas(sigmas), indx(indx),
826  sigcut(cut){
827  assert(sigmas.size() == epsilons.size());
828  };
829  IEpsISigCutAtom(AtomID a, IEpsISigCutAtom other) : AtomID(a), epsilons(other.epsilons),
830  sigmas(other.sigmas), indx(other.indx), sigcut(other.sigcut){};
832  assert(other.indx < epsilons.size());
833  flt myeps = epsilons[other.indx];
834  //~ assert(indx < other.epsilons.size());
835  //~ assert(other.epsilons[indx] == myeps);
836  return myeps;
837  }
839  assert(other.indx < sigmas.size());
840  flt myeps = sigmas[other.indx];
841  //~ assert(indx < other.sigmas.size());
842  //~ assert(other.sigmas[indx] == myeps);
843  return myeps;
844  }
846  flt sigma = sigmas[0];
847  for(uint i=1; i<sigmas.size(); ++i){
848  if(sigma < sigmas[i]) sigma = sigmas[i];
849  }
850  return sigma*sigcut;
851  };
852 };
853 
863  AtomID atom1, atom2;
865  inter(sqrt(a1.epsilon * a2.epsilon),
866  (a1.sigma + a2.sigma) / 2,
867  max(a1.sigcut, a2.sigcut)),
868  atom1(a1), atom2(a2){};
870  inter(a1.get_epsilon(a2),
871  a1.get_sigma(a2),
872  max(a1.sigcut, a2.sigcut)),
873  atom1(a1), atom2(a2){};
874 
875  inline flt energy(Box &box){return inter.energy(box.diff(atom1->x, atom2->x));};
876  inline Vec forces(Box &box){return inter.forces(box.diff(atom1->x, atom2->x));};
877 };
878 
879 struct IEpsSigCutAtom : public AtomID {
880  vector<flt> epsilons; // for finding epsilons
881  uint indx; // for finding this one in other atoms' epsilon lists
882  // note that for two HydroAtoms a1, a2, the epsilon for the pair
883  // is then either a1.epsilons[a2.indx] or a2.epsilons[a1.indx]
885  flt sigcut; // sigma units
887  IEpsSigCutAtom(AtomID a, vector<flt> epsilons, uint indx, flt sigma, flt cut) :
888  AtomID(a), epsilons(epsilons), indx(indx), sigma(sigma),
889  sigcut(cut){};
890  IEpsSigCutAtom(AtomID a, IEpsSigCutAtom other) : AtomID(a), epsilons(other.epsilons),
891  indx(other.indx), sigma(other.sigma), sigcut(other.sigcut){};
893  assert(other.indx < epsilons.size());
894  flt myeps = epsilons[other.indx];
895  assert(indx < other.epsilons.size());
896  assert(other.epsilons[indx] == myeps);
897  return myeps;
898  }
899  flt max_size(){return sigma*sigcut;};
900 };
901 
904  AtomID atom1, atom2;
906  inter(sqrt(a1.epsilon * a2.epsilon),
907  (a1.sigma + a2.sigma) / 2,
908  max(a1.sigcut, a2.sigcut)),
909  atom1(a1), atom2(a2){};
911  inter(a1.get_epsilon(a2),
912  a1.get_sigma(a2),
913  max(a1.sigcut, a2.sigcut)),
914  atom1(a1), atom2(a2){};
916  inter(a1.get_epsilon(a2),
917  (a1.sigma + a2.sigma) / 2,
918  max(a1.sigcut, a2.sigcut)),
919  atom1(a1), atom2(a2){};
920  inline flt energy(Box &box){return inter.energy(box.diff(atom1->x, atom2->x));};
921  inline Vec forces(Box &box){return inter.forces(box.diff(atom1->x, atom2->x));};
922 };
923 
925 // Full LJish, with ε = ε₁₂ (indexed) and σ = (σ₁+σ₂)
926 // Potential is V(r) = ε (σ^n/r^n - 1)² - ε (r₀^-n - 1)²
927 // cutoff at some sigcut r₀
928 
929 struct IEpsRepsSigExpCutAtom : public AtomID {
930  vector<flt> epsilons; // for finding epsilons
931  flt repeps, sigma;
932  flt exponent; // power
933  uint indx; // for finding this one in other atoms' epsilon lists
934  // note that for two atoms a1, a2, the epsilon for the pair
935  // is then either a1.epsilons[a2.indx] or a2.epsilons[a1.indx]
936  // these should be the same
937  flt sigcut; // sigma units
939  IEpsRepsSigExpCutAtom(AtomID a, vector<flt> epsilons, flt repeps, flt sigma,
940  flt n, uint indx, flt cut) :
941  AtomID(a), epsilons(epsilons), repeps(repeps), sigma(sigma),
942  exponent(n), indx(indx), sigcut(cut){
943  assert(indx < epsilons.size());
944  //~ assert(sigma > 0.01);
945  };
947  AtomID(a), epsilons(other.epsilons), repeps(other.repeps),
948  sigma(other.sigma), exponent(other.exponent), indx(other.indx),
949  sigcut(other.sigcut){
950  //~ assert(other.sigma > 0.01);
951  //~ assert(sigma > 0.01);
952  };
954  assert(other.indx < epsilons.size());
955  flt myeps = epsilons[other.indx];
956  return myeps;
957  }
959  return (sigma + other.sigma)/2.0;
960  }
961  flt max_size(){return sigma*sigcut;};
962 
963 };
964 
965 struct LJishPair {
966  flt epsilon, repeps, sigma, n, cut_distance, cut_energy;
967  AtomID atom1, atom2;
969  epsilon(LJ1.get_epsilon(LJ2)),
970  repeps(sqrt(LJ1.repeps * LJ2.repeps)),
971  sigma(LJ1.get_sigma(LJ2)),
972  n((LJ1.exponent + LJ2.exponent) / 2),
973  cut_distance(max(LJ1.sigcut, LJ2.sigcut)),
974  atom1(LJ1), atom2(LJ2){
975  if(epsilon <= 0){
976  cut_distance = 1;
977  cut_energy = 0;
978  epsilon = 0;
979  return;
980  }
981  flt mid = (1-pow(cut_distance,-n));
982  cut_energy = epsilon*(mid*mid);
983  };
984  inline flt energy(Box &box){
985  Vec rij = box.diff(atom1->x, atom2->x);
986  flt rsq = rij.squaredNorm()/(sigma*sigma);
987  if(rsq > cut_distance * cut_distance){
988  //~ printf("LJish: dx=%.2f, σ=%.2f, rsq=%.2f, cutR²=%.2f\n", rij.norm(), sigma, rsq, cut_distance);
989  return 0;
990  }
991 
992  flt mid = (1-pow(rsq,-n/2));
993  //~ printf("LJish: rsq=%.2f, mid %.2f, epsilon %.2f, repeps %.2f\n", rsq, mid, epsilon, repeps);
994  if (rsq > 1) return epsilon*(mid*mid) - cut_energy;
995  return repeps*(mid*mid) - cut_energy;
996  };
997  inline Vec forces(Box &box){
998  Vec rij = box.diff(atom1->x, atom2->x);
999  flt dsq = rij.squaredNorm();
1000  flt rsq = dsq/(sigma*sigma);
1001  if(rsq > cut_distance * cut_distance) return Vec::Zero();
1002  flt rmid = pow(rsq,-n/2); // σ^n/r^n
1003  flt fmagTimesR = 2*n*rmid*(rmid - 1); // 2 n r^-n(r^-n - 1)
1004  if (rsq < 1) return rij * (repeps * fmagTimesR / dsq);
1005  return rij * (epsilon * fmagTimesR / dsq); // r⃗ * 2 n * r^-n(r^-n - 1) / r² = 2 n r^-(n-1) * (r^-n - 1) r̂
1006  }
1007 };
1008 
1010 // Full LJ, with ε = ε₁₂ (indexed) and σ = (σ₁ + σ₂)/2
1011 // Potential is V(r) = ε (σ⁶/r⁶ - 1)² - ε (r₀⁻⁶ - 1)²
1012 // cutoff at some sigcut r₀
1013 // if ε₁₂ < 0, purely repulsive with ε = -ε₁₂
1014 
1016  flt eps, sig;
1017  flt cut_distance, cut_energy;
1018  AtomID atom1, atom2;
1020  eps(a1.get_epsilon(a2)), sig((a1.sigma + a2.sigma)/2.0),
1021  cut_distance(max(a1.sigcut, a2.sigcut)),
1022  atom1(a1), atom2(a2){
1023  if(eps <= 0){
1024  cut_distance = 1;
1025  cut_energy = 0;
1026  eps = abs(eps);
1027  return;
1028  }
1029  flt mid = (1-pow(cut_distance,-6));
1030  cut_energy = eps*(mid*mid);
1031  };
1032  inline flt energy(Box &box){
1033  Vec rij = box.diff(atom1->x, atom2->x);
1034  flt rsq = rij.squaredNorm()/(sig*sig);
1035  if(rsq > cut_distance*cut_distance) {
1036  //~ printf("Distance: %.2f Energy: %.2f (ε: %.2f σ: %.2f cut: %.2f cut_energy: %.2f)\n",
1037  //~ sqrt(rij.squaredNorm()), 0.0, eps, sig, cut_distance, cut_energy);
1038  return 0;
1039  }
1040  flt mid = (1-pow(rsq,-3));
1041  //~ if(eps*(mid*mid) - cut_energy < 0) {
1042  //~ printf("Distance: %.2f Energy: %.2f (ε: %.2f σ: %.2f cut: %.2f cut_energy: %.2f)\n",
1043  //~ rij.norm(), eps*(mid*mid) - cut_energy, eps, sig, cut_distance, cut_energy);
1044  //~ }
1045  return eps*(mid*mid) - cut_energy;
1046  };
1047  inline Vec forces(Box &box){
1048  if(eps == 0) return Vec::Zero();
1049  Vec rij = box.diff(atom1->x, atom2->x);
1050  flt dsq = rij.squaredNorm();
1051  flt rsq = dsq/(sig*sig);
1052  if(rsq > (cut_distance*cut_distance)) return Vec::Zero();
1053  flt rsix = pow(rsq,-3); // σ⁶/r⁶
1054  flt fmagTimesR = 12*eps*rsix*(rsix - 1);
1055  return rij * (fmagTimesR / dsq);
1056  };
1057 };
1058 
1060 // Full LJ, with ε_A = ε₁₂ (indexed) and σ = (σ₁ + σ₂)/2
1061 // Attractive Potential is V(r) = ε (σ⁶/r⁶ - 1)² - ε (r₀⁻⁶ - 1)²
1062 // cutoff at some sigcut r₀
1063 // For repulsive, ε_R = repeps
1064 
1065 struct IEpsRepsSigCutAtom : public AtomID {
1066  vector<flt> epsilons; // for finding epsilons
1067  flt repeps, sig;
1069  flt sigcut; // sigma units
1071  IEpsRepsSigCutAtom(AtomID a, vector<flt> epsilons, flt repeps,
1072  flt sigma, uint indx, flt cut) :
1073  AtomID(a), epsilons(epsilons), repeps(repeps), sig(sigma),
1074  indx(indx), sigcut(cut){
1075  assert(indx < epsilons.size());
1076  };
1078  AtomID(a), epsilons(other.epsilons),
1079  repeps(other.repeps), sig(other.sig), indx(other.indx), sigcut(other.sigcut){};
1081  assert(other.indx < epsilons.size());
1082  flt myeps = epsilons[other.indx];
1083  return myeps;
1084  }
1085  flt max_size(){return sig*sigcut;};
1086 };
1087 
1088 // r < 1 εr (1 - σ⁶/r⁶)² - εa (1 - σ⁶/R⁶)²
1089 // 1 < r < R εa (1 - σ⁶/r⁶)² - εa (1 - σ⁶/R⁶)²
1090 // R < r 0
1091 
1092 // so cut_energy = εa (1 - σ⁶/R⁶)²
1094  flt eps, repeps, sig;
1095  flt cut_distance, cut_energy;
1096  bool attract;
1097  AtomID atom1, atom2;
1100  eps(abs(a1.get_epsilon(a2))), repeps(sqrt(a1.repeps * a2.repeps)),
1101  sig((a1.sig + a2.sig)/2.0),
1102  cut_distance(max(a1.sigcut, a2.sigcut)), attract(a1.get_epsilon(a2) > 0),
1103  atom1(a1), atom2(a2){
1104  if(!attract or eps == 0){
1105  cut_distance = 1;
1106  cut_energy = 0;
1107  eps = 0;
1108  attract = false;
1109  return;
1110  }
1111  flt mid = (1-pow(cut_distance,-6.0));
1112  cut_energy = eps*(mid*mid);
1113  };
1114  inline flt energy(Box &box){
1115  Vec rij = box.diff(atom1->x, atom2->x);
1116  flt rsq = rij.squaredNorm()/(sig*sig);
1117  if(rsq > cut_distance*cut_distance) {
1118  //~ printf("Distance: %.2f Energy: %.2f (ε: %.2f σ: %.2f cut: %.2f cut_energy: %.2f)\n",
1119  //~ sqrt(rij.squaredNorm()), 0.0, eps, sig, cut_distance, cut_energy);
1120  return 0;
1121  }
1122  flt mid = (1-pow(rsq,-3)); // # 1 - σ⁶/r⁶
1123  //~ printf("Distance: %.2f Energy: %.2f (ε: %.2f σ: %.2f cut: %.2f cut_energy: %.2f)\n",
1124  //~ sqrt(rij.squaredNorm()), eps*(mid*mid) - cut_energy, eps, sig, cut_distance, cut_energy);
1125  if (rsq > 1) return eps*(mid*mid) - cut_energy;
1126  return repeps*(mid*mid) - cut_energy;
1127  //~ flt E;
1128  //~ if (rsq > 1) E = eps*(mid*mid) - cut_energy;
1129  //~ else E = repeps*(mid*mid) - cut_energy;
1130  //~ if(E > 1e4)
1131  //~ printf("Distance: %.2f Energy: %.2f (ε: %.2f,%.2f σ: %.2f cut: %.2f cut_energy: %.2f)\n",
1132  //~ sqrt(rij.squaredNorm()), E, eps, repeps, sig, cut_distance, cut_energy);
1133  //~ return E;
1134  };
1135  inline Vec forces(Box &box){
1136  Vec rij = box.diff(atom1->x, atom2->x);
1137  flt dsq = rij.squaredNorm();
1138  flt rsq = dsq/(sig*sig);
1139  if(rsq > (cut_distance*cut_distance)) return Vec::Zero();
1140  flt rsix = pow(rsq,-3); // σ⁶/r⁶
1141  flt fmagTimesR = 12*rsix*(rsix - 1);
1142  if (rsq < 1) return rij * (repeps * fmagTimesR / dsq);
1143  return rij * (eps * fmagTimesR / dsq);
1144  //~ flt fmag;
1145  //~ if (rsq < 1) fmag = repeps * fmagTimesR / dsq;
1146  //~ else fmag = eps * fmagTimesR / dsq;
1147  //~ if(fmag * rij.norm() > 1e4)
1148  //~ printf("Distance: %.2f Force: %.2f (ε: %.2f,%.2f σ: %.2f cut: %.2f cut_energy: %.2f)\n",
1149  //~ sqrt(rij.squaredNorm()), fmag * rij.norm(), eps, repeps, sig, cut_distance, cut_energy);
1150  //~ return rij * fmag;
1151  };
1152 };
1153 
1154 struct EisMclachlanAtom : public AtomID {
1155  flt dist, sigmai; // dist is radius of Atom + radius of water (1.4 Å)
1157  EisMclachlanAtom(AtomID a, flt dist, flt sigmai) : AtomID(a),
1158  dist(dist), sigmai(sigmai){};
1160  dist(other.dist), sigmai(other.sigmai){};
1161  flt max_size(){return dist;};
1162 };
1163 
1165  flt c0,c1,c2; // 1/R term, const term, R term in E (coefficients)
1166  // E = c₀/R + c₁ + c₂ R
1167  flt cutoff; // r1 + r2 (includes two water radii)
1168  AtomID atom1, atom2;
1170  c0(-M_PI*(a1.sigmai*a2.dist - a2.sigmai*a1.dist)
1171  * (a1.dist*a1.dist - a2.dist*a2.dist)),
1172  c1(-2*M_PI*(a1.sigmai*a2.dist*a2.dist + a2.sigmai*a1.dist*a1.dist)),
1173  c2(M_PI*(a1.sigmai*a2.dist + a2.sigmai*a1.dist)),
1174  cutoff(a1.dist + a2.dist),
1175  atom1(a1), atom2(a2){};
1176  inline flt energy(Box &box){
1177  Vec rij = box.diff(atom1->x, atom2->x);
1178  flt R = rij.norm();
1179  if (R > cutoff) return 0;
1180  return c0/R + c1 + c2*R;
1181  }
1182  inline Vec forces(Box &box){
1183  Vec rij = box.diff(atom1->x, atom2->x);
1184  flt dsq = rij.squaredNorm();
1185  if(dsq > (cutoff*cutoff)) return Vec::Zero();
1186  flt R = sqrt(dsq);
1187  return rij * ((c0/dsq-c2)/R);
1188  }
1189 };
1190 
1191 struct EpsSigExpAtom : public AtomID {
1192  flt eps, sigma, exponent;
1194  EpsSigExpAtom(AtomID a, flt eps, flt sigma, flt exponent) : AtomID(a),
1195  eps(eps), sigma(sigma), exponent(exponent){};
1197  eps(other.eps), sigma(other.sigma), exponent(other.exponent){};
1198  flt max_size(){return sigma;};
1199 };
1200 
1201 struct EnergyForce {
1204  EnergyForce(Vec f, flt E) : f(f), E(E){};
1205 };
1206 
1207 //----------------------------------------------------------------------
1208 // Repulsion potential, as above, with ε = ε₁₂ and σ = σ₁₂ (both indexed)
1209 // exponent is n = (n₁ + n₂)/2
1210 
1211 
1212 struct IEpsISigExpAtom : public AtomID {
1213  vector<flt> epsilons; // for finding epsilons
1214  vector<flt> sigmas; // for finding epsilons
1216  uint indx; // for finding this one in other atoms' epsilon lists
1217  // note that for two IEpsISigExpAtom atoms a1, a2, the epsilon for the pair
1218  // is then either a1.epsilons[a2.indx] or a2.epsilons[a1.indx]; same for sigma
1220  IEpsISigExpAtom(AtomID a, vector<flt> epsilons, vector<flt> sigmas,
1221  uint indx, flt exponent=2.5) :
1222  AtomID(a), epsilons(epsilons), sigmas(sigmas), exponent(exponent),
1223  indx(indx){
1224  assert(sigmas.size() == epsilons.size());
1225  };
1226  IEpsISigExpAtom(AtomID a, IEpsISigCutAtom other) : AtomID(a), epsilons(other.epsilons),
1227  sigmas(other.sigmas), indx(other.indx){};
1229  assert(other.indx < epsilons.size());
1230  flt myeps = epsilons[other.indx];
1231  //~ assert(indx < other.epsilons.size());
1232  //~ assert(other.epsilons[indx] == myeps);
1233  return myeps;
1234  }
1236  assert(other.indx < sigmas.size());
1237  flt myeps = sigmas[other.indx];
1238  //~ assert(indx < other.sigmas.size());
1239  //~ assert(other.sigmas[indx] == myeps);
1240  return myeps;
1241  }
1243  flt sigma = sigmas[0];
1244  for(uint i=1; i<sigmas.size(); ++i){
1245  if(sigma < sigmas[i]) sigma = sigmas[i];
1246  }
1247  return sigma;
1248  };
1249 };
1250 
1255 
1257  flt eps, sig, exponent;
1258  AtomID atom1, atom2;
1260  eps(sqrt(a1.eps * a2.eps)), sig((a1.sigma + a2.sigma)/2.0),
1261  exponent((a1.exponent + a2.exponent)/2.0), atom1(a1), atom2(a2){};
1263  eps(a1.get_epsilon(a2)), sig(a1.get_sigma(a2)),
1264  exponent((a1.exponent + a2.exponent)/2.0), atom1(a1), atom2(a2){};
1265  inline flt energy(Box &box){
1266  Vec rij = box.diff(atom1->x, atom2->x);
1267  flt dsq = rij.squaredNorm();
1268  if(dsq > sig*sig) return 0.0;
1269  flt R = sqrt(dsq);
1270  return eps * pow(1.0 - (R/sig), exponent) / exponent;
1271  }
1272  inline Vec forces(Box &box){
1273  Vec rij = box.diff(atom1->x, atom2->x);
1274  flt dsq = rij.squaredNorm();
1275  if(dsq > sig*sig) return Vec::Zero();
1276  flt R = sqrt(dsq);
1277  return rij * (eps * pow(1.0 - (R/sig), exponent-1) /sig/R);
1278  }
1280  Vec rij = box.diff(atom1->x, atom2->x);
1281  flt dsq = rij.squaredNorm();
1282  if(dsq > sig*sig) return EnergyForce(Vec::Zero(),0);
1283  flt R = sqrt(dsq);
1284 
1285  Vec f = rij * (eps * pow(1.0 - (R/sig), exponent-1) /sig/R);
1286  flt E = eps * pow(1.0 - (R/sig), exponent) / exponent;
1287  return EnergyForce(f, E);
1288  }
1289  //~ inline flt xrij(Box &box){
1290  //~ Vec rij = box.diff(atom1->x, atom2->x);
1291  //~ flt dsq = rij.squaredNorm();
1292  //~ if(dsq > sig*sig) return 0.0;
1293  //~ flt R = sqrt(dsq);
1294  //~ return (R*eps*(exponent-1)/sig/sig) * pow(1.0 - (R/sig), exponent-2);
1295  //~ }
1296  inline void fill(Box &box, ForcePairX &fpair){
1297  fpair.a1 = atom1;
1298  fpair.a2 = atom2;
1299  Vec rij = box.diff(atom1->x, atom2->x);
1300  flt dsq = rij.squaredNorm();
1301  if(dsq > sig*sig){
1302  fpair.fij = Vec::Zero();
1303  fpair.xij = 0.0;
1304  return;
1305  }
1306  flt R = sqrt(dsq);
1307  fpair.fij = rij * (eps*pow(1.0 - (R/sig), exponent-1) /sig/R);
1308  flt Rs = R / sig;
1309  fpair.xij = -Rs*eps*pow(1.0 - (R/sig), exponent-2)*(1-exponent*Rs);
1310  }
1311 };
1312 
1318 
1319 struct EpsSigExpDragAtom : public AtomID {
1320  flt eps, sigma, exponent, gamma;
1322  EpsSigExpDragAtom(AtomID a, flt eps, flt sigma, flt gamma, flt exponent=2.5) : AtomID(a),
1323  eps(eps), sigma(sigma), exponent(exponent), gamma(gamma){};
1325  eps(other.eps), sigma(other.sigma), exponent(other.exponent), gamma(other.gamma){};
1326  flt max_size(){return sigma;};
1327 };
1328 
1330  flt eps, sig, exponent, gamma;
1331  AtomID atom1, atom2;
1333  eps(sqrt(a1.eps * a2.eps)), sig((a1.sigma + a2.sigma)/2.0),
1334  exponent((a1.exponent + a2.exponent)/2.0),
1335  gamma((a1.gamma + a2.gamma)/2), atom1(a1), atom2(a2){};
1336  inline flt energy(Box &box){
1337  Vec rij = box.diff(atom1->x, atom2->x);
1338  flt dsq = rij.squaredNorm();
1339  if(dsq > sig*sig) return 0.0;
1340  flt R = sqrt(dsq);
1341  return eps * pow(1.0 - (R/sig), exponent) / exponent;
1342  }
1343  inline Vec forces(Box &box){
1344  Vec rij = box.diff(atom1->x, atom2->x);
1345  Vec vij = atom1->v - atom2->v;
1346  flt dsq = rij.squaredNorm();
1347  if(dsq > sig*sig) return Vec::Zero();
1348  flt R = sqrt(dsq);
1349  Vec v_perp = rij * (vij.dot(rij)) / dsq;
1350  return rij * (eps * pow(1.0 - (R/sig), exponent-1) /sig/R) -
1351  (v_perp * gamma);
1352  }
1353 };
1354 
1355 /***********************************************************************
1356  * Lois-O'Hern potential (that's what I'm calling it) for approximating
1357  * sticky spheres. Force is linear repulsive, linear attractive, with
1358  * a repulsive epsilon, a sigma, and an attractive epsilon and
1359  * lengthscale.
1360  * See PRL 100, 028001 (2008)
1361  * They used C = 10^-2, l = 0. How do you have l=0? Well, that's a
1362  * discontinuity in the force, but not in the energy, so it works.
1363  *
1364  * More specifically:
1365 F=\begin{cases}
1366 Y\left(1-\frac{r}{\sigma}\right) & \frac{r}{\sigma}<1+C\\
1367 \frac{CY}{\ell}\left(\frac{r}{\sigma}-\left(1+C+\ell\right)\right) & 1+C\leq\frac{r}{\sigma}\leq1+C+\ell
1368 \end{cases}
1369 
1370 U\left(r\right)=\begin{cases}
1371 -\frac{Y\sigma}{2}\left(C\left(C+\ell\right)-\left(\frac{r}{\sigma}-1\right)^{2}\right) & \frac{r}{\sigma}<1+C\\
1372 -\frac{CY\sigma}{2\ell}\left(C+\ell+1-\frac{r}{\sigma}\right)^{2} & 1+C\leq\frac{r}{\sigma}\leq1+C+\ell
1373 \end{cases}
1374 
1375 The minimum is at σ, its normal harmonic from 0 to (1+C)σ, and is then
1376 "rounded" between (1+C)σ < r < (1+C+l)σ.
1377 * With l = 0, its a harmonic potential, depth C²σ/2, minimum at σ, attractive
1378 * out to (1+C)σ.
1379 
1380 Parameters:
1381 * C: how wide the "first half" of the well is
1382 * l: how wide the "second half" of the well is
1383 * l/C: how steep the "second half" of the well is; l == C means perfectly harmonic
1384 * C(C+l)σ/2: depth of the well
1385 
1386 */
1387 
1388 struct LoisOhernAtom : public AtomID {
1389  flt eps, sigma, C, l;
1391  LoisOhernAtom(AtomID a, flt eps, flt sigma, flt C, flt l) : AtomID(a),
1392  eps(eps), sigma(sigma), C(C), l(l){};
1394  eps(other.eps), sigma(other.sigma), C(other.C), l(other.l){};
1395  flt max_size(){return sigma*(1+C+l);};
1396 };
1397 
1399  flt eps, sig, C, l, sigcut;
1400  AtomID atom1, atom2;
1402  eps(sqrt(a1.eps * a2.eps)), sig((a1.sigma + a2.sigma)/2.0),
1403  C((a1.C + a2.C)/2.0), l((a1.l + a2.l)/2.0), sigcut(sig*(1+C+l)),
1404  atom1(a1), atom2(a2){};
1406  eps(eps), sig(sig), C(C), l(l), sigcut(sig*(1+C+l)),
1407  atom1(a1), atom2(a2){};
1408  inline flt energy(Box &box){
1409  Vec rij = box.diff(atom1->x, atom2->x);
1410  flt dsq = rij.squaredNorm();
1411  // using >= to prevent NaNs when l = 0
1412  if(dsq >= sigcut*sigcut) return 0.0;
1413  flt R = sqrt(dsq)/sig;
1414  // using <= to prevent NaNs when l = 0
1415  if(R <= 1 + C){
1416  flt dR = R-1;
1417  return -eps*sig/2*(C*(C+l) - dR*dR);
1418  }
1419 
1420  flt dR2 = C+l+1 - R;
1421  return -C*eps*sig/2/l*dR2*dR2;
1422  }
1423 
1424  inline Vec forces(Box &box){
1425  Vec rij = box.diff(atom1->x, atom2->x);
1426  flt dsq = rij.squaredNorm();
1427  if(dsq >= sigcut*sigcut) return Vec::Zero();
1428  flt R = sqrt(dsq);
1429  flt rsig = R/sig;
1430 
1431  if(rsig <= 1 + C){
1432  flt dR = rsig-1;
1433  return rij*(-eps*dR/R);
1434  }
1435 
1436  flt dR2 = rsig-(C+l+1);
1437  return rij*(C*eps/l*dR2/R);
1438  }
1439 };
1440 
1441 
1442 
1445  LoisOhernPair(a1, a2, sqrt(a1.eps * a2.eps),
1446  (a1.sigma + a2.sigma)/2.0,
1447  (a1.C < a2.C ? a1.C : a2.C),
1448  (a1.l < a2.l ? a1.l : a2.l)){};
1449 };
1450 
1452 /***********************************************************************
1453  * LoisLin atoms
1454  * Harmonic repulsive, fixed-force attractive
1455  * epsilon is the strength of the harmonic Interaction
1456  * f is the fixed-force amount (real units)
1457  * l is the width of the fixed-force region (real units)
1458  * To combine, we use a geometric average of epsilon, and everything else
1459  * is averaged.
1460  */
1461 
1462 struct LoisLinAtom : public AtomID {
1463  flt eps, sigma, f, l;
1465  LoisLinAtom(AtomID a, flt eps, flt sigma, flt depth, flt width) : AtomID(a),
1466  eps(eps), sigma(sigma), f(width > 0 ? depth/width : 0), l(width){};
1468  eps(other.eps), sigma(other.sigma), f(other.f), l(other.l){};
1469  flt max_size(){return sigma*(1+l);};
1470 };
1471 
1472 struct LoisLinPair {
1473  flt eps, sig, f, l, sigcut;
1474  AtomID atom1, atom2;
1476  eps(sqrt(a1.eps * a2.eps)), sig((a1.sigma + a2.sigma)/2.0),
1477  f((a1.f + a2.f)/2.0), l((a1.l + a2.l)/2.0), sigcut(sig +l),
1478  atom1(a1), atom2(a2){};
1479  LoisLinPair(LoisLinAtom a1, LoisLinAtom a2, flt eps, flt sig, flt f, flt l) :
1480  eps(eps), sig(sig), f(f), l(l), sigcut(sig+l),
1481  atom1(a1), atom2(a2){};
1482  inline flt energy(Box &box){
1483  Vec rij = box.diff(atom1->x, atom2->x);
1484  flt dsq = rij.squaredNorm();
1485  // using >= to prevent NaNs when l = 0
1486  if(dsq >= sigcut*sigcut) return 0.0;
1487  flt R = sqrt(dsq);
1488  // using <= to prevent NaNs when l = 0
1489  if(R <= sig){
1490  flt dR = 1.0 - (R/sig);
1491  return eps/2 * dR*dR - f*l;
1492  }
1493 
1494  return -f*(sig + l - R);
1495  }
1496 
1497  inline Vec forces(Box &box){
1498  Vec rij = box.diff(atom1->x, atom2->x);
1499  flt dsq = rij.squaredNorm();
1500  if(dsq >= sigcut*sigcut) return Vec::Zero();
1501  flt R = sqrt(dsq);
1502 
1503  if(R <= sig){
1504  flt dR = 1.0 - (R/sig);
1505  return rij * (eps * dR /sig/R);
1506  }
1507 
1508  return -rij*(f/R);
1509  }
1510 };
1511 
1512 struct LoisLinPairMin : public LoisLinPair {
1514  LoisLinPair(a1, a2, sqrt(a1.eps * a2.eps),
1515  (a1.sigma + a2.sigma)/2.0,
1516  (a1.f < a2.f ? a1.f : a2.f),
1517  (a1.l < a2.l ? a1.l : a2.l)){};
1518 };
1519 
1520 template <class A, class P>
1521 class SCBoxed : public Interaction {
1522  protected:
1523  sptr<SCBox> box;
1524  sptr<AtomVec> atoms;
1525  vector<A> group; // data on which atoms interact with the wall
1526  public:
1527  SCBoxed(sptr<AtomVec> atomv, sptr<SCBox> box)
1528  : box(box), atoms(atomv){};
1529  inline void add(A atm){group.push_back(atm);};
1530  flt energy(Box &box);
1531  flt pressure(Box &box);
1532  void set_forces(Box &box);
1533  flt set_forces_get_pressure(Box &box);
1534  inline vector<A> &atom_list(){return group;};
1535 };
1536 
1537 template <class A, class P>
1538 class SimpleListed : public Interaction {
1539  protected:
1540  vector<A> atoms;
1541 
1542  public:
1544  inline void add(A atm){atoms.push_back(atm);};
1545  //flt energy(Box &box, IDPair &pair);
1546  flt energy(Box &box);
1547  flt pressure(Box &box);
1548  uint size(){return ((uint) (atoms.size()));};
1549  //inline flt energy_pair(P pair, Box &box){return pair.energy(box);}; // This may need to be written!
1550  void set_forces(Box &box);
1551  flt set_forces_get_pressure(Box &box);
1552  //inline Vec forces_pair(P pair, Box &box){return pair.forces(box);}; // This may need to be written!
1553  inline vector<A> &atom_list(){return atoms;};
1555 };
1556 
1557 template <class A, class P>
1558 class NListed : public Interaction {
1559  // NeighborList keeps track of Atom* pairs within a particular distance.
1560  // NListed keeps track of atoms with additional properties
1561  // that interact through a particular Interaction.
1562  // class A needs to inherit from AtomID, and also have an initialization
1563  // method A(AtomID a, A other)
1564  // you also need to implement a couple methods below
1565 
1566  // Implementation detail:
1567  // note that NeighborList maintains its own MetaGroup, so that when
1568  // a member of A is created and passed to this group, the AtomID in that
1569  // member has an A.n() value referring to its place in the *original*
1570  // AtomGroup, not this one. This is why we need an A(AtomID, A) method;
1571  // so we can make a new A with all the same properties as before,
1572  // but with the n() referring to the NeighborList maintained AtomGroup.
1573  protected:
1574  vector<A> atoms; // NOT all full.
1575  // This is the length of the AtomVec, and if
1576  // Atom i has been added, then atoms[i] is
1577  // filled in this vector
1578  vector<P> pairs;
1579  sptr<NeighborList> neighbors;
1581  public:
1582  NListed(sptr<AtomVec> vec, sptr<NeighborList> neighbors) :
1583  atoms(vec->size()), neighbors(neighbors), last_update(0){}; //group(vec),
1584  inline void add(A atm){
1585  neighbors->add(atm, atm.max_size());
1586  // assert(atoms.size() > atm.n());
1587  atoms[atm.n()] = atm;
1588  }
1589  NListed(sptr<Box> box, sptr<AtomVec> atomv, const flt skin) :
1590  atoms(atomv->size()),
1591  neighbors(new NeighborList(box, atomv, skin)), last_update(0){};
1592  void update_pairs();
1593  P get_pair(IDPair &pair){
1594  return P(atoms[pair.first().n()], atoms[pair.last().n()]);}
1595  A& getatom(uint n){return atoms[n];}
1596  flt energy(Box &box, IDPair &pair);
1597  flt energy(Box &box);
1598 
1600  unsigned long long contacts(Box &box);
1602  unsigned long long overlaps(Box &box);
1603  flt pressure(Box &box);
1604  inline vector<P> &pair_iter(){return pairs;};
1605  uint size(){return ((uint) (atoms.size()));};
1606  inline flt energy_pair(P pair, Box &box){return pair.energy(box);}; // This may need to be written!
1607  void set_forces(Box &box);
1608  flt set_forces_get_pressure(Box &box);
1609  Matrix stress(Box &box);
1610  Matrix set_forces_get_stress(Box &box);
1611  inline Vec forces_pair(P pair, Box &box){return pair.forces(box);}; // This may need to be written!
1612  inline vector<A> &atom_list(){return atoms;};
1613  inline sptr<NeighborList> neighbor_list(){return neighbors;};
1614  //~ flt energy_test(flt dist);
1615  //~ flt force_test(flt dist);
1617 };
1618 
1619 template <class A, class P>
1621  private:
1622  NListed<A,P> nlisted;
1623  public:
1624  NListedVirial(sptr<AtomVec> vec, sptr<NeighborList> neighbors) :
1625  nlisted(vec, neighbors){};
1626  void set_forces(Box &box){nlisted.set_forces(box);};
1627  void set_forces(Box &box, FPairXFunct*);
1628  virtual inline flt set_forces_get_pressure(Box &box){return nlisted.set_forces_get_pressure(box);};
1629  virtual flt setForcesGetEnergy(Box &box);
1630  virtual inline flt energy(Box &box){return nlisted.energy(box);};
1631  virtual inline flt pressure(Box &box){return nlisted.pressure(box);};
1632  inline void add(A atm){nlisted.add(atm);}
1633  inline sptr<NeighborList> neighbor_list(){return nlisted.neighbor_list();};
1634 };
1635 
1636 template <class A, class P>
1638  flt E = 0;
1639  Atom wallatom;
1640  wallatom.m = NAN;
1641  AtomID wallid = AtomID(&wallatom, atoms->size());
1642  typename vector<A>::iterator it;
1643  for(it = group.begin(); it != group.end(); ++it){
1644  Vec r0 = (*it)->x;
1645  Vec edger = box->edge_dist(r0);
1646  wallatom.x = r0 + (edger*2);
1647  E += P(*it, EpsSigExpAtom(wallid, *it)).energy(newbox) / 2; // no energy from fake Atom
1648  }
1649  return E;
1650 };
1651 
1652 template <class A, class P>
1654  Atom wallatom;
1655  wallatom.m = NAN;
1656  AtomID wallid = AtomID(&wallatom, atoms->size());
1657  typename vector<A>::iterator it;
1658  for(it = group.begin(); it != group.end(); ++it){
1659  Vec r0 = (*it)->x;
1660  Vec edger = box->edge_dist(r0);
1661  wallatom.x = r0 + (edger*2);
1662  P pair = P(*it, EpsSigExpAtom(wallid, *it));
1663  Vec f = pair.forces(*box);
1664  (*it)->f += f;
1665  }
1666 };
1667 
1668 template <class A, class P>
1670  flt p = 0;
1671  Atom wallatom;
1672  wallatom.m = NAN;
1673  AtomID wallid = AtomID(&wallatom, atoms->size());
1674  typename vector<A>::iterator it;
1675  for(it = group.begin(); it != group.end(); ++it){
1676  Vec r0 = (*it)->x;
1677  Vec dr = box->edge_dist(r0) * 2;
1678  wallatom.x = r0 + dr;
1679  P pair = P(*it, EpsSigExpAtom(wallid, *it));
1680  Vec f = pair.forces(*box);
1681  (*it)->f += f;
1682  p += f.dot(dr);
1683  }
1684  return p;
1685 };
1686 
1687 template <class A, class P>
1689  flt p = 0;
1690  Atom wallatom;
1691  wallatom.m = NAN;
1692  AtomID wallid = AtomID(&wallatom, atoms->size());
1693  typename vector<A>::iterator it;
1694  for(it = group.begin(); it != group.end(); ++it){
1695  Vec r0 = (*it)->x;
1696  Vec dr = box->edge_dist(r0) * 2;
1697  wallatom.x = r0 + dr;
1698  P pair = P(*it, EpsSigExpAtom(wallid, *it));
1699  Vec f = pair.forces(*box);
1700  p += f.dot(dr);
1701  }
1702  return p;
1703 };
1704 
1705 template <class A, class P>
1707  flt E = 0;
1708  typename vector<A>::iterator it1;
1709  typename vector<A>::iterator it2;
1710  for(it1 = atoms.begin(); it1 != atoms.end(); ++it1){
1711  for(it2 = it1+1; it2 != atoms.end(); ++it2){
1712  E += P(*it1, *it2).energy(box);
1713  }
1714  }
1715  return E;
1716 };
1717 
1718 template <class A, class P>
1720  typename vector<A>::iterator it1;
1721  typename vector<A>::iterator it2;
1722  for(it1 = atoms.begin(); it1 != atoms.end(); ++it1){
1723  for(it2 = it1+1; it2 != atoms.end(); ++it2){
1724  P pair = P(*it1, *it2);
1725  Vec f = pair.forces(box);
1726  (*it1)->f += f;
1727  (*it2)->f -= f;
1728  }
1729  }
1730 };
1731 
1732 template <class A, class P>
1734  flt p=0;
1735  typename vector<A>::iterator it1;
1736  typename vector<A>::iterator it2;
1737  for(it1 = atoms.begin(); it1 != atoms.end(); ++it1){
1738  for(it2 = it1+1; it2 != atoms.end(); ++it2){
1739  P pair = P(*it1, *it2);
1740  Vec f = pair.forces(box);
1741  (*it1)->f += f;
1742  (*it2)->f -= f;
1743  Vec r = box.diff((*it1)->x, (*it2)->x);
1744  p += r.dot(f);
1745  }
1746  }
1747  //~ cout << "Set forces, got pressure" << p << '\n';
1748  return p;
1749 };
1750 
1751 template <class A, class P>
1753  flt p=0;
1754  typename vector<A>::iterator it1;
1755  typename vector<A>::iterator it2;
1756  for(it1 = atoms.begin(); it1 != atoms.end(); ++it1){
1757  for(it2 = it1+1; it2 != atoms.end(); ++it2){
1758  P pair = P(*it1, *it2);
1759  Vec f = pair.forces(box);
1760  Vec r = box.diff((*it1)->x, (*it2)->x);
1761  p += r.dot(f);
1762  }
1763  }
1764  return p;
1765 };
1766 
1767 template <class A, class P>
1769  if(last_update == neighbors->which()) return; // already updated
1770 
1771  last_update = neighbors->which();
1772  pairs.clear();
1773  vector<IDPair>::iterator pairit;
1774  for(pairit = neighbors->begin(); pairit != neighbors->end(); ++pairit){
1775  // assert(atoms.size() > pairit->first().n());
1776  // assert(atoms.size() > pairit->last().n());
1777  A firstatom = atoms[pairit->first().n()];
1778  A secondatom = atoms[pairit->last().n()];
1779  pairs.push_back(P(firstatom, secondatom));
1780  }
1781 }
1782 
1783 template <class A, class P>
1785  update_pairs(); // make sure the LJpairs match the neighbor list ones
1786  P Epair = P(atoms[pair.first().n()],atoms[pair.last().n()]);
1787  //~ Vec dist = box.diff(Epair.atom1->x, Epair.atom2->x);
1788  return energy_pair(Epair, box);
1789 };
1790 
1791 template <class A, class P>
1792 unsigned long long NListed<A, P>::contacts(Box &box){
1793  unsigned long long Nc = 0;
1794  update_pairs(); // make sure the LJpairs match the neighbor list ones
1795  typename vector<P>::iterator it;
1796  for(it = pairs.begin(); it != pairs.end(); ++it){
1797  flt E = energy_pair(*it, box);
1798  if(E != 0.0){Nc++;};
1799  }
1800  return Nc;
1801 };
1802 
1803 template <class A, class P>
1804 unsigned long long NListed<A, P>::overlaps(Box &box){
1805  unsigned long long Nc = 0;
1806  update_pairs(); // make sure the LJpairs match the neighbor list ones
1807  typename vector<P>::iterator it;
1808  for(it = pairs.begin(); it != pairs.end(); ++it){
1809  flt E = energy_pair(*it, box);
1810  if(E > 0.0){Nc++;};
1811  }
1812  return Nc;
1813 };
1814 
1815 template <class A, class P>
1817  update_pairs(); // make sure the LJpairs match the neighbor list ones
1818  flt E = 0;
1819  typename vector<P>::iterator it;
1820  for(it = pairs.begin(); it != pairs.end(); ++it){
1821  //~ Vec dist = box.diff(it->atom1->x, it->atom2->x);
1822  E += energy_pair(*it, box);
1823  }
1824  return E;
1825 };
1826 
1827 template <class A, class P>
1829  update_pairs(); // make sure the LJpairs match the neighbor list ones
1830  typename vector<P>::iterator it;
1831  for(it = pairs.begin(); it != pairs.end(); ++it){
1832  Vec f = forces_pair(*it, box);
1833  it->atom1->f += f;
1834  it->atom2->f -= f;
1835  //~ assert(f.squaredNorm() < 1000000);
1836  }
1837 };
1838 
1839 template <class A, class P>
1841  nlisted.update_pairs(); // make sure the LJpairs match the neighbor list ones
1842  flt E = 0;
1843  vector<P> & pairs = nlisted.pair_iter();
1844  for(typename vector<P>::iterator it = pairs.begin(); it != pairs.end(); ++it){
1845  EnergyForce EF = it->get_energy_forces(box);
1846  it->atom1->f += EF.f;
1847  it->atom2->f -= EF.f;
1848  E += EF.E;
1849  }
1850  return E;
1851 }
1852 
1853 //~ template <class A, class P>
1854 //~ flt NListed<A, P>::energy_test(flt dist){
1855  //~ typename vector<P>::iterator it = pairs.begin();
1856  //~ Atom a1 = Atom(it->first());
1857  //~ Atom a2 = Atom(it->second());
1858  //~ A atm1 = A(it->first(), &a1);
1859  //~ A atm2 = A(it->second(), &a1);
1860  //~ a1.x = Vec::Zero();
1861  //~ a1.v = Vec::Zero();
1862  //~ a1.a = Vec::Zero();
1863  //~ a1.f = Vec::Zero();
1864  //~
1865  //~ a2.x = Vec::Zero();
1866  //~ a2.v = Vec::Zero();
1867  //~ a2.a = Vec::Zero();
1868  //~ a2.f = Vec::Zero();
1869  //~ a2.x.setx(dist);
1870  //~ P Epair = P(atm1,atm2);
1871  //~ return energy_pair(Epair, infbox);
1872 //~ };
1873 
1874 template <class A, class P>
1876  nlisted.update_pairs(); // make sure the LJpairs match the neighbor list ones
1877  ForcePairX myfpair;
1878  typename vector<P>::iterator it;
1879  for(it = nlisted.pair_iter().begin(); it != nlisted.pair_iter().end(); ++it){
1880  it->fill(box, myfpair);
1881  //~ assert(!isnan(myfpair.fij[0]));
1882  funct->run(&myfpair);
1883  it->atom1->f += myfpair.fij;
1884  it->atom2->f -= myfpair.fij;
1885  //~ assert(f.squaredNorm() < 1000000);
1886  }
1887 };
1888 
1889 template <class A, class P>
1891  update_pairs(); // make sure the LJpairs match the neighbor list ones
1892  flt p=0;
1893  typename vector<P>::iterator it;
1894  for(it = pairs.begin(); it != pairs.end(); ++it){
1895  Vec f = forces_pair(*it, box);
1896  it->atom1->f += f;
1897  it->atom2->f -= f;
1898  Vec r = box.diff(it->atom1->x, it->atom2->x);
1899  p += r.dot(f);
1900  }
1901  //~ cout << "Set forces, got pressure" << p << '\n';
1902  return p;
1903 };
1904 
1905 template <class A, class P>
1907  update_pairs(); // make sure the LJpairs match the neighbor list ones
1908  flt p=0;
1909  //~ printf("Updating pressure...\n");
1910  typename vector<P>::iterator it;
1911  for(it = pairs.begin(); it != pairs.end(); ++it){
1912  Vec f = forces_pair(*it, box);
1913  Vec r = box.diff(it->atom1->x, it->atom2->x);
1914  //Vec r = it->atom1->x - it->atom2->x;
1915  p += r.dot(f);
1916  //~ assert(!isnan(p));
1917  }
1918  return p;
1919 };
1920 
1921 template <class A, class P>
1923  update_pairs(); // make sure the LJpairs match the neighbor list ones
1924  Matrix stress = Matrix::Zero();
1925  typename vector<P>::iterator it;
1926  for(it = pairs.begin(); it != pairs.end(); ++it){
1927  Vec f = forces_pair(*it, box);
1928  it->atom1->f += f;
1929  it->atom2->f -= f;
1930  Vec r = box.diff(it->atom1->x, it->atom2->x);
1931  stress += r * f.transpose();
1932  }
1933  //~ cout << "Set forces, got pressure" << p << '\n';
1934  return stress;
1935 };
1936 
1937 template <class A, class P>
1939  update_pairs(); // make sure the LJpairs match the neighbor list ones
1940  Matrix stress = Matrix::Zero();
1941  typename vector<P>::iterator it;
1942  for(it = pairs.begin(); it != pairs.end(); ++it){
1943  Vec f = forces_pair(*it, box);
1944  Vec r = box.diff(it->atom1->x, it->atom2->x);
1945  stress += r * f.transpose();
1946  }
1947  //~ cout << "Set forces, got pressure" << p << '\n';
1948  return stress;
1949 };
1950 
1951 class Charges : public Interaction {
1952  protected:
1953  vector<Charged> atoms;
1957  AtomID get_id(Atom* a);
1958 
1959  public:
1960  Charges(flt screenlength, flt k=1, vector<Charged> atms=vector<Charged>());
1961  // cutoffdist in sigma units
1962 
1963  inline void add(Charged a){atoms.push_back(a); return;};
1964  inline void add(AtomID a, flt q){add(Charged(q,a));};
1965  inline void ignore(AtomID a, AtomID b){ignorepairs.add_pair(a,b);};
1966  inline void ignore(Atom* a, Atom* b){
1967  ignore(get_id(a),get_id(b));};
1968  inline uint ignore_size() const{return ignorepairs.size();};
1969  inline uint size() const{return (uint) atoms.size();};
1970  flt energy(Box &box);
1971  flt pressure(Box &box);
1972  void set_forces(Box &box);
1973  //~ ~LJsimple(){};
1974 };
1975 
1977 // A wall
1978 
1979 struct WallAtom : AtomRef {
1980  public:
1983  public:
1984  WallAtom(AtomID a, flt sigma, flt epsilon=1.0) :
1985  AtomRef(a), sigma(sigma), epsilon(epsilon){};
1986 };
1987 
1988 class SoftWall : public Interaction {
1989  protected:
1993  vector<WallAtom> group;
1995  public:
1996  SoftWall(Vec loc, Vec norm, flt expt=2.0) :
1997  loc(loc), norm(norm.normalized()), expt(expt), lastf(NAN){};
1998  void add(WallAtom a){group.push_back(a);};
1999  flt energy(Box &box);
2000  void set_forces(Box &box);
2001  flt set_forces_get_pressure(Box &box);
2002  flt pressure(Box &box);
2003 
2004  void set_location(Vec newloc){loc = newloc;};
2005  Vec get_location(){return loc;};
2006  void set_normal(Vec newNorm){norm = newNorm.normalized();};
2007  Vec get_normal(){return norm;};
2008 
2009  flt get_last_f(){return lastf;};
2010 };
2011 
2013  protected:
2019  vector<WallAtom> group;
2020  public:
2021  SoftWallCylinder(Vec loc, Vec axis, flt radius, flt expt=2.0) :
2022  loc(loc), axis(axis.normalized()), radius(radius), expt(expt), lastf(NAN){};
2023  void add(WallAtom a){
2024  if(a.sigma > radius*2)
2025  throw std::invalid_argument("SoftWallCylinder::add: sigma must be less than cylinder diameter");
2026  group.push_back(a);
2027  };
2028  flt energy(Box &box);
2029  void set_forces(Box &box);
2030  flt set_forces_get_pressure(Box &box);
2031  flt pressure(Box &box);
2032 
2033  void set_location(Vec new_loc){loc = new_loc;};
2034  Vec get_location(){return loc;};
2035  void set_axis(Vec new_axis){axis = new_axis.normalized();};
2036  Vec get_axis(){return axis;};
2037  flt get_last_f(){return lastf;};
2038 };
2039 
2040 #ifdef VEC2D
2041 class WalledBox2D : public OriginBox {
2042  protected:
2043  bool xwalls, ywalls;
2044  vector<SoftWall*> walls;
2045  public:
2046  WalledBox2D(Vec size, bool xwalled, bool ywalled, flt expt=2.0) :
2047  OriginBox(size), xwalls(xwalled), ywalls(ywalled){
2048  if(xwalls){
2049  walls.push_back(new SoftWall(
2050  Vec(-size[0]/2.0, 0), Vec(1, 0), expt));
2051  walls.push_back(new SoftWall(
2052  Vec(size[0]/2.0, 0), Vec(-1, 0), expt));
2053  }
2054  if(ywalls){
2055  walls.push_back(new SoftWall(
2056  Vec(0, -size[0]/2.0), Vec(0, 1), expt));
2057  walls.push_back(new SoftWall(
2058  Vec(0, size[0]/2.0), Vec(0, -1), expt));
2059  }
2060  };
2061  Vec diff(Vec r1, Vec r2){
2062  Vec dr = r1 - r2;
2063  if(!xwalls) dr[0] = remainder(r1[0], boxsize[0]);
2064  if(!ywalls) dr[1] = remainder(r1[1], boxsize[1]);
2065  return dr;
2066  }
2067  Vec diff(Vec r1, Vec r2, array<int,NDIM> boxes){
2068  Vec dr = r1 - r2;
2069  if(!xwalls) dr[0] -= boxes[0]*boxsize[0];
2070  if(!ywalls) dr[1] -= boxes[1]*boxsize[1];
2071  return dr;
2072  }
2073  flt V(){return boxsize[0] * boxsize[1];};
2074  flt L(){return (boxsize[0] + boxsize[1])/2.0;};
2075 
2076  flt resize(flt factor){
2077  boxsize *= factor;
2078  vector<SoftWall*>::iterator it;
2079  for(it = walls.begin(); it != walls.end(); ++it){
2080  (*it)->set_location((*it)->get_location() * factor);
2081  };
2082  return V();
2083  }
2085  flt curV = V();
2086  return resize(pow(newV/curV, OVERNDIM));
2087  }
2088  Vec rand_loc(flt walldist){
2089  Vec bxvec = boxsize;
2090  if(xwalls) bxvec[0] -= walldist/2.0;
2091  if(ywalls) bxvec[1] -= walldist/2.0;
2092  Vec v = rand_vec_boxed();
2093  for(uint i=0; i<NDIM; ++i){
2094  v[i] *= bxvec[i];
2095  }
2096  return diff(v, Vec::Zero());
2097  };
2098  vector<SoftWall*> get_walls() {return walls;};
2100  for(vector<SoftWall*>::iterator it = walls.begin(); it != walls.end(); ++it){
2101  delete *it;
2102  }
2103  };
2104 };
2105 #endif
2106 
2107 inline flt confine_range(flt minimum, flt val, flt maximum){
2108  if(val <= minimum) return minimum;
2109  if(val >= maximum) return maximum;
2110  return val;
2111 }
2112 
2113 class SCAtomVec : public virtual AtomGroup {
2114  // this is an AtomGroup which actually owns the atoms, which are
2115  // arranged in pairs.
2116  private:
2117  AtomVec atoms;
2118  public:
2119  SCAtomVec(vector<double> masses) : atoms((uint) masses.size()*2, 0.0){
2120  for(uint i=0; i < masses.size(); ++i){
2121  atoms[2*i].m = masses[i]/2;
2122  atoms[2*i+1].m = masses[i]/2;
2123  }
2124  };
2125  SCAtomVec(uint N, flt mass) : atoms(N*2, mass/2.0){};
2126  inline AtomVec &vec(){return atoms;};
2127  inline Atom& operator[](cuint n){return atoms[n];};
2128  inline Atom& operator[](cuint n) const {return atoms[n];};
2129  inline AtomID get_id(cuint n){return atoms.get_id(n);};
2130  inline IDPair pair(cuint n){return IDPair(atoms.get_id(n*2), atoms.get_id(n*2 + 1));};
2131  inline uint size() const {return atoms.size();};
2132  inline uint pairs() const {return atoms.size()/2;};
2136  static flt volume(flt diameter, flt length, uint dim=NDIM);
2138 };
2139 
2141  Vec delta, r;
2142  flt lambda1, lambda2;
2143 };
2144 
2145 struct SCPair {
2148  flt l1, l2;
2149  SCPair(IDPair &p1, IDPair &p2, flt l1, flt l2) :
2150  p1(p1), p2(p2), l1(l1), l2(l2){};
2151  SCPair(IDPair &p1, IDPair &p2, flt l) :
2152  p1(p1), p2(p2), l1(l), l2(l){};
2153  SCPair(const SCPair &other) : p1(other.p1), p2(other.p2),
2154  l1(other.l1), l2(other.l2){}
2155  SpheroCylinderDiff nearest_location(Box &box);
2156  void apply_force(Box &box, Vec f, SpheroCylinderDiff diff, flt I1, flt I2);
2157  inline void apply_force(Box &box, Vec f, SpheroCylinderDiff diff, flt I){
2158  return apply_force(box, f, diff, I, I);
2159  }
2160 };
2161 
2162 struct SCSpringPair : public SCPair {
2164  flt eps, sig;
2165 
2166  SCSpringPair(IDPair &p1, IDPair &p2, flt eps, flt sig, flt l1, flt l2) :
2167  SCPair(p1, p2, l1, l2), eps(eps), sig(sig){};
2168  SCSpringPair(IDPair &p1, IDPair &p2, flt eps, flt sig, flt l) :
2169  SCPair(p1, p2, l), eps(eps), sig(sig){};
2170 
2171  inline flt max_dist(){return sig + (l1+l2)/2;};
2172  inline flt max_delta(){return sig;};
2173 
2175  flt dsq = diff.delta.squaredNorm();
2176  if(dsq > sig*sig) return 0;
2177  flt d = sqrt(dsq);
2178  flt dsig = d-sig;
2179  return dsig*dsig*eps/2;
2180  }
2182  flt dsq = diff.delta.squaredNorm();
2183  if(dsq > sig*sig) return Vec::Zero();
2184  flt dmag = sqrt(dsq);
2185  Vec dhat = diff.delta / dmag;
2186 
2187  return dhat * (eps * (sig - dmag));
2188  };
2189 };
2190 
2191 class SCSpringList : public Interaction {
2193  private:
2194  SCAtomVec *scs;
2195  flt eps, sig;
2196  vector<flt> ls;
2197  set<array<uint, 2> > ignore_list;
2198  public:
2201  SCSpringList(SCAtomVec *scs, flt eps, flt sig, flt l) :
2202  scs(scs), eps(eps), sig(sig), ls(scs->pairs(), l){};
2203  SCSpringList(SCAtomVec *scs, flt eps, flt sig, vector<flt> ls) :
2204  scs(scs), eps(eps), sig(sig), ls(ls){};
2205  flt energy(Box &box);
2206  void set_forces(Box &box);
2207  flt set_forces_get_pressure(Box &box);
2208  flt pressure(Box &box);
2209  Matrix set_forces_get_stress(Box &box);
2210  Matrix stress(Box &box);
2211  void ignore(uint n1, uint n2){
2212  if(n1 > n2){uint n3=n1; n1=n2; n2=n3;}
2213  array<uint, 2> pair;
2214  pair[0] = n1;
2215  pair[1] = n2;
2216  ignore_list.insert(pair);
2217  }
2218  flt volume();
2219  flt phi(Box &box){return volume() / box.V();};
2220  void ignore(AtomID a1, AtomID a2){ignore(a1.n() / 2, a2.n() / 2);}
2221 
2223 };
2224 
2225 #endif
void ignore(AtomID a, AtomID b)
Definition: interaction.hpp:1965
LJAttractCutPair(EpsSigCutAtom a1, EpsSigCutAtom a2)
Definition: interaction.hpp:905
Vec com() const
center of mass
Definition: box.cpp:222
virtual void run(ForcePairX *)=0
flt get_epsilon(IEpsISigExpAtom &other)
Definition: interaction.hpp:1228
void set_axis(Vec new_axis)
Definition: interaction.hpp:2035
flt radius
Definition: interaction.hpp:2016
void add(A atm)
Definition: interaction.hpp:1529
AtomVec & vec()
Definition: interaction.hpp:2126
Vec rand_loc(flt walldist)
Definition: interaction.hpp:2088
Definition: interaction.hpp:686
void add_forced(AngleGrouping b)
Definition: interaction.hpp:664
Definition: interaction.hpp:392
Vec vec()
A one-liner for creating a Vec object, occasionally useful from within Python.
Definition: vecrand.hpp:72
uint indx
Definition: interaction.hpp:819
uint n() const
Definition: box.hpp:279
flt l2
Definition: interaction.hpp:2148
flt sig
Definition: interaction.hpp:1016
vector< A > atoms
Definition: interaction.hpp:1574
vector< flt > epsilons
Definition: interaction.hpp:817
~SCSpringList()
Definition: interaction.hpp:2222
The basic class for representing each particle.
Definition: box.hpp:229
Vec forces(Box &box, SpheroCylinderDiff diff)
Definition: interaction.hpp:2181
static Vec forces(const Vec diff, const flt eps, const flt sig, const flt cutsig)
Definition: interaction.hpp:247
flt pressure(Box &box)
Partial pressure due to this Interaction.
Definition: interaction.hpp:1688
uint size()
Definition: interaction.hpp:1605
AtomID get_id(cuint n)
Definition: interaction.hpp:2129
SCSpringPair(IDPair &p1, IDPair &p2, flt eps, flt sig, flt l1, flt l2)
Definition: interaction.hpp:2166
void set_normal(Vec newNorm)
Definition: interaction.hpp:2006
flt cut_energy
Definition: interaction.hpp:1017
flt confine_range(flt minimum, flt val, flt maximum)
Definition: interaction.hpp:2107
Definition: interaction.hpp:536
LJAttractRepulsePair(IEpsSigCutAtom a1, IEpsSigCutAtom a2)
Definition: interaction.hpp:1019
bool attract
Definition: interaction.hpp:1096
IEpsISigCutAtom(AtomID a, IEpsISigCutAtom other)
Definition: interaction.hpp:829
Definition: interaction.hpp:365
LennardJonesCut(const flt epsilon, const flt sigma, const flt cutsig)
Definition: interaction.hpp:234
virtual flt pressure(Box &box)
Partial pressure due to this Interaction.
Definition: interaction.hpp:1631
IEpsRepsSigExpCutAtom(AtomID a, IEpsRepsSigExpCutAtom other)
Definition: interaction.hpp:946
IDPair & p1
Definition: interaction.hpp:2146
void add(FixedForceAtom a)
Definition: interaction.hpp:370
LennardJonesCut inter
Definition: interaction.hpp:862
LJAttractCut(const flt epsilon, const flt sigma, const flt cutsig)
Definition: interaction.hpp:196
void set_force(Box &box)
Definition: interaction.hpp:447
flt epsilon
Definition: interaction.hpp:192
vector< WallAtom > group
Definition: interaction.hpp:2019
Vec axis
Definition: interaction.hpp:2015
flt energy(Box &box)
Potential energy due to this Interaction.
Definition: interaction.hpp:1637
flt sigcut
Definition: interaction.hpp:1399
RandomForceAtom(AtomID a, flt force_mag, flt freq, RandomForceType force_type=UNIFORM)
Definition: interaction.hpp:547
flt resize(flt factor)
Definition: interaction.hpp:2076
flt exponent
Definition: interaction.hpp:1215
Repulsion potential, with ε = √(ε₁ ε₂) and σ = (σ₁ + σ₂)/2 Potential is V(r) = ε/n (1 - r/σ)^n...
Definition: interaction.hpp:1256
IEpsISigExpAtom()
Definition: interaction.hpp:1219
vector< flt > sigmas
Definition: interaction.hpp:818
flt energy(Box &box)
Definition: interaction.hpp:984
vector< flt > epsilons
Definition: interaction.hpp:930
FixedForceAtom(Vec F, AtomID a)
Definition: interaction.hpp:360
void set_location(Vec new_loc)
Definition: interaction.hpp:2033
flt m
mass
Definition: box.hpp:243
void add(flt k, flt theta0, AtomID a1, AtomID a2, AtomID a3, AtomID a4)
Add 4 atoms with the potential .
Definition: interaction.hpp:698
IEpsSigCutAtom(AtomID a, IEpsSigCutAtom other)
Definition: interaction.hpp:890
uint size() const
Definition: interaction.hpp:377
LennardJonesCutPair(EpsSigCutAtom a1, EpsSigCutAtom a2)
Definition: interaction.hpp:864
a group of atoms, such as all of them (AtomVec), or a smaller group such as a molecule, sidebranch, etc.
Definition: box.hpp:312
void add(flt k, AtomID a1, AtomID a2, AtomID a3, AtomID a4)
Add 4 atoms with the potential , where is set to match the current positions of the atoms...
Definition: interaction.hpp:709
EnergyForce(Vec f, flt E)
Definition: interaction.hpp:1204
void add_forced(BondGrouping b)
Definition: interaction.hpp:622
Definition: interaction.hpp:114
Definition: interaction.hpp:1558
A & getatom(uint n)
Definition: interaction.hpp:1595
uint indx
Definition: interaction.hpp:933
EpsSigExpAtom()
Definition: interaction.hpp:1193
uint last_update
Definition: interaction.hpp:1580
SCPair(IDPair &p1, IDPair &p2, flt l1, flt l2)
Definition: interaction.hpp:2149
flt cut_energy
Definition: interaction.hpp:1095
AtomID a1
Definition: interaction.hpp:595
Atom a2
Definition: interaction.hpp:728
Definition: interaction.hpp:590
void set_forces(Box &box)
Definition: interaction.hpp:493
void ignore(uint n1, uint n2)
Definition: interaction.hpp:2211
void ignore(Atom *a, Atom *b)
Definition: interaction.hpp:1966
void add(FixedForceRegionAtom a)
Definition: interaction.hpp:411
flt cutoffE
Definition: interaction.hpp:335
Definition: interaction.hpp:593
Maintains a Verlet list of "neighbors": molecules within a 'skin radius' of each other.
Definition: trackers.hpp:118
flt sigma
Definition: interaction.hpp:1320
FixedSpringAtom(AtomID a, Vec loc, flt k, bool usex=true, bool usey=true, bool usez=true)
Definition: interaction.hpp:434
Repulsion potential with drag, with ε = √(ε₁ ε₂) and σ = (σ₁ + σ₂)/2 Potential is V(r) = ε/n (1 ...
Definition: interaction.hpp:1319
void set_forces(Box &box)
Definition: interaction.hpp:473
NListed(sptr< Box > box, sptr< AtomVec > atomv, const flt skin)
Definition: interaction.hpp:1589
flt pressure(Box &box)
Partial pressure due to this Interaction.
Definition: interaction.hpp:1906
flt L()
Definition: interaction.hpp:2074
virtual ~InteractPair()
Definition: interaction.hpp:108
uint size() const
Definition: interaction.hpp:715
flt max_size()
Definition: interaction.hpp:1395
LennardJonesCutPair(IEpsISigCutAtom a1, IEpsISigCutAtom a2)
Definition: interaction.hpp:869
Vec fij
Definition: interaction.hpp:729
flt sigcut
Definition: interaction.hpp:1069
Vec forces(const Vec &diff)
Definition: interaction.hpp:144
The main class for representing particles.
Definition: box.hpp:412
unsigned int uint
Definition: vec.hpp:20
array< Vec, 4 > derivs
Definition: interaction.hpp:299
flt max_delta()
Definition: interaction.hpp:2172
LoisOhernAtom(AtomID a, flt eps, flt sigma, flt C, flt l)
Definition: interaction.hpp:1391
vector< SoftWall * > get_walls()
Definition: interaction.hpp:2098
Vec x
location.
Definition: box.hpp:231
Definition: interaction.hpp:480
void set_forces(Box &box)
Definition: interaction.hpp:384
flt max_size()
Definition: interaction.hpp:899
EpsSigExpAtom(AtomID a, flt eps, flt sigma, flt exponent)
Definition: interaction.hpp:1194
flt energy(Box &box)
Definition: interaction.hpp:875
Vec delta
Definition: interaction.hpp:2141
flt epsilon
Definition: interaction.hpp:149
double flt
The basic floating point type used in the simulations.
Definition: vecrand.hpp:45
AtomID atom2
Definition: interaction.hpp:1018
flt k
Definition: interaction.hpp:1956
sptr< NeighborList > neighbor_list()
Definition: interaction.hpp:1633
AtomID a2
Definition: interaction.hpp:595
SCPair(IDPair &p1, IDPair &p2, flt l)
Definition: interaction.hpp:2151
Definition: interaction.hpp:303
BondAngle(const flt k, const flt theta, const bool cosine=false)
Definition: interaction.hpp:284
Definition: interaction.hpp:2041
Definition: interaction.hpp:1472
LJAttractFixedRepulsePair()
Definition: interaction.hpp:1098
uint indx
Definition: interaction.hpp:1068
LoisLinPair(LoisLinAtom a1, LoisLinAtom a2)
Definition: interaction.hpp:1475
Vec forces(const Vec &diff)
Definition: interaction.hpp:221
LoisLinAtom()
Definition: interaction.hpp:1464
WallAtom(AtomID a, flt sigma, flt epsilon=1.0)
Definition: interaction.hpp:1984
void add(AtomID a, flt q)
Definition: interaction.hpp:1964
bool add(flt k, flt x0, AtomID a1, AtomID a2, bool replace=true)
Definition: interaction.hpp:620
Definition: interaction.hpp:456
IEpsRepsSigCutAtom()
Definition: interaction.hpp:1070
Vec loc
Definition: interaction.hpp:430
sptr< NeighborList > neighbors
Definition: interaction.hpp:1579
LoisOhernPair(LoisOhernAtom a1, LoisOhernAtom a2, flt eps, flt sig, flt C, flt l)
Definition: interaction.hpp:1405
vector< FixedForceAtom > atoms
Definition: interaction.hpp:367
IEpsISigCutAtom()
Definition: interaction.hpp:823
~SimpleListed()
Definition: interaction.hpp:1554
flt pressure(Box &box)
Partial pressure due to this Interaction.
Definition: interaction.hpp:667
vector< A > & atom_list()
Definition: interaction.hpp:1534
Definition: interaction.hpp:298
flt cutoff
Definition: interaction.hpp:1167
BondDiffType
Definition: interaction.hpp:587
Definition: interaction.hpp:1398
flt max_size()
Definition: interaction.hpp:961
flt phi(Box &box)
Definition: interaction.hpp:2219
flt energy(Box &box)
Definition: interaction.hpp:1176
LoisOhernAtom(AtomID a, LoisOhernAtom other)
Definition: interaction.hpp:1393
uint size() const
Definition: interaction.hpp:629
A rectilinear Box, with periodic boundary conditions.
Definition: box.hpp:98
virtual flt V()=0
Volume. Can return NaN.
Atom & operator[](cuint n)
Definition: interaction.hpp:2127
sptr< AtomVec > atoms
Definition: interaction.hpp:1524
Definition: interaction.hpp:650
Definition: interaction.hpp:1212
AtomID a4
Definition: interaction.hpp:679
Definition: interaction.hpp:929
flt energy(Box &box)
Definition: interaction.hpp:1482
Vec forces(const Vec &diff)
Definition: interaction.hpp:186
uint size() const
Number of atoms in the group.
Definition: interaction.hpp:2131
flt x0
Definition: interaction.hpp:640
flt max_size()
Definition: interaction.hpp:773
Definition: interaction.hpp:589
Definition: interaction.hpp:1164
Definition: interaction.hpp:760
Definition: interaction.hpp:1191
void add(vector< flt > cos_coefficients, vector< flt > sincoeffs, AtomID a1, AtomID a2, AtomID a3, AtomID a4, bool usepow=true)
Definition: interaction.hpp:694
flt epsilon
Definition: interaction.hpp:1982
flt max_dist()
Definition: interaction.hpp:2171
AtomID atom2
Definition: interaction.hpp:1331
Vec f
Definition: interaction.hpp:1202
flt sigma
Definition: interaction.hpp:117
flt max_size()
Definition: interaction.hpp:1242
A mapping of Atom -> [list of Atom], used by NeighborList to keep track of what atoms are near what o...
Definition: trackers.hpp:29
Definition: interaction.hpp:1093
SCAtomVec(vector< double > masses)
Definition: interaction.hpp:2119
flt energy(const Vec r)
Definition: interaction.hpp:339
Matrix set_forces_get_stress(Box &box)
Definition: interaction.hpp:1922
flt pressure(Box &box)
Partial pressure due to this Interaction.
Definition: interaction.hpp:426
void add(AtomID a, Vec dir, vector< flt > bound, vector< flt > F)
Definition: interaction.hpp:412
static flt energy(const Vec diff, const flt eps, const flt sig)
Definition: interaction.hpp:121
flt energy(const Vec &diff1, const Vec &diff2, const Vec &diff3) const
Definition: interaction.hpp:322
flt expt
Definition: interaction.hpp:2017
AtomID atom2
Definition: interaction.hpp:1474
bool add(flt k, AtomID a1, AtomID a2, bool replace=true)
Add a pair of atoms with the current distance.
Definition: interaction.hpp:624
static Vec diff(Vec r1, Vec r2)
Definition: interaction.hpp:677
~SCAtomVec()
Definition: interaction.hpp:2137
Vec get_location()
Definition: interaction.hpp:2005
SCPair(const SCPair &other)
Definition: interaction.hpp:2153
void update_pairs()
Definition: interaction.hpp:1768
EisMclachlanAtom()
Definition: interaction.hpp:1156
static Vec diff(Vec r1, Vec r2)
Definition: interaction.hpp:653
Vec forces(Box &box)
Definition: interaction.hpp:1497
~WalledBox2D()
Definition: interaction.hpp:2099
flt sigma
Definition: interaction.hpp:1463
Definition: interaction.hpp:766
LJAttractCut inter
Definition: interaction.hpp:903
virtual uint size() const =0
Number of atoms in the group.
flt E
Definition: interaction.hpp:1203
flt screen
Definition: interaction.hpp:332
A pointer to an Atom, that also knows its own index in an AtomVec.
Definition: box.hpp:270
Definition: interaction.hpp:538
vector< WallAtom > group
Definition: interaction.hpp:1993
flt max_size()
Definition: interaction.hpp:845
IEpsISigCutAtom(AtomID a, vector< flt > epsilons, vector< flt > sigmas, uint indx, flt cut)
Definition: interaction.hpp:824
Vec fij
Definition: interaction.hpp:735
static Vec forces(const Vec diff, const flt eps, const flt sig)
Definition: interaction.hpp:171
flt energy(Box &box)
Potential energy due to this Interaction.
Definition: interaction.hpp:571
unsigned long long overlaps(Box &box)
number of Atom pairs with E > 0
Definition: interaction.hpp:1804
Definition: interaction.hpp:1065
flt max_size()
Definition: interaction.hpp:807
static flt forces(const flt rsig)
Definition: interaction.hpp:179
RepulsionPair(IEpsISigExpAtom a1, IEpsISigExpAtom a2)
Definition: interaction.hpp:1262
~NListed()
Definition: interaction.hpp:1616
flt energy(Box &box)
Potential energy due to this Interaction.
Definition: interaction.hpp:467
uint indx
Definition: interaction.hpp:881
vector< flt > boundaries
Definition: interaction.hpp:394
Vec forces(Box &box)
Definition: interaction.hpp:921
flt energy(Box &box)
Definition: interaction.hpp:920
Vec r
Definition: interaction.hpp:2141
flt energy(const Vec &diff)
Definition: interaction.hpp:129
~Spring()
Definition: interaction.hpp:273
EpsSigExpDragAtom()
Definition: interaction.hpp:1321
AtomID a
Definition: interaction.hpp:433
RepulsionDragPair(EpsSigExpDragAtom a1, EpsSigExpDragAtom a2)
Definition: interaction.hpp:1332
void set_forces(Box &box)
Definition: interaction.hpp:1626
void add(Vec F, AtomID a)
Definition: interaction.hpp:371
AtomID atom2
Definition: interaction.hpp:863
RandomForce()
Definition: interaction.hpp:556
flt energy_pair(P pair, Box &box)
Definition: interaction.hpp:1606
void add(DihedralGrouping b)
Definition: interaction.hpp:691
vector< flt > sincoeffs
Definition: interaction.hpp:307
void add(A atm)
Definition: interaction.hpp:1632
flt epsilon
Definition: interaction.hpp:116
Charged(flt q, AtomID a)
Definition: interaction.hpp:757
Definition: interaction.hpp:816
LJAttract(const flt epsilon, const flt sigma)
Definition: interaction.hpp:152
bool zeropressure
Definition: interaction.hpp:609
Vec forces(Box &box)
Definition: interaction.hpp:1343
flt freq
Definition: interaction.hpp:544
flt springk
Definition: interaction.hpp:267
AtomID a2
Definition: interaction.hpp:641
Atom & operator[](cuint n) const
Definition: interaction.hpp:2128
flt sigcut
Definition: interaction.hpp:801
LoisOhernPair(LoisOhernAtom a1, LoisOhernAtom a2)
Definition: interaction.hpp:1401
static flt get_angle(const Vec &r1, const Vec &r2)
Definition: interaction.hpp:286
Vec forces(Box &box)
Definition: interaction.hpp:876
uint pairs() const
Definition: interaction.hpp:2132
uint size() const
Number of atoms in the group.
Definition: box.hpp:444
void add(flt x, flt y, flt z, AtomID a)
Definition: interaction.hpp:373
LJRepulsive(const flt epsilon, const flt sigma)
Definition: interaction.hpp:119
IEpsRepsSigExpCutAtom(AtomID a, vector< flt > epsilons, flt repeps, flt sigma, flt n, uint indx, flt cut)
Definition: interaction.hpp:939
Definition: interaction.hpp:1462
flt V()
Volume. Can return NaN.
Definition: interaction.hpp:2073
AtomID atom2
Definition: interaction.hpp:1168
flt theta0
Definition: interaction.hpp:280
flt max_size()
Definition: interaction.hpp:1198
static flt energy(const Vec diff, const flt eps, const flt sig, const flt cutsig)
Definition: interaction.hpp:199
void add(WallAtom a)
Definition: interaction.hpp:1998
vector< flt > cos_coefficients
Definition: interaction.hpp:306
uint size() const
Definition: interaction.hpp:1969
flt force_mag
Definition: interaction.hpp:543
void set_force(Box &box)
Definition: interaction.hpp:362
LoisOhernPairMinCLs(LoisOhernAtom a1, LoisOhernAtom a2)
Definition: interaction.hpp:1444
SoftWall(Vec loc, Vec norm, flt expt=2.0)
Definition: interaction.hpp:1996
RandomForceType force_type
Definition: interaction.hpp:545
Charged()
Definition: interaction.hpp:756
EisMclachlanAtom(AtomID a, EisMclachlanAtom other)
Definition: interaction.hpp:1159
void add(A atm)
Definition: interaction.hpp:1584
flt resize_to_V(flt newV)
Definition: interaction.hpp:2084
flt sig
Definition: interaction.hpp:1257
flt get_last_f()
Definition: interaction.hpp:2037
flt sigma
Definition: interaction.hpp:150
flt epsilon
Definition: interaction.hpp:230
Definition: interaction.hpp:403
flt sigma
Definition: interaction.hpp:193
flt set_forces_get_pressure(Box &box)
Set forces (Atom.f) and return at the same time (see pressure()).
Definition: interaction.hpp:1733
flt get_epsilon(IEpsRepsSigCutAtom &other)
Definition: interaction.hpp:1080
Vec loc
Definition: interaction.hpp:2014
Definition: interaction.hpp:265
flt energy(Box &box)
Definition: interaction.hpp:790
Definition: interaction.hpp:639
AtomID atom2
Definition: interaction.hpp:967
LJAttractCutPair(IEpsSigCutAtom a1, IEpsSigCutAtom a2)
Definition: interaction.hpp:915
flt cutoff
Definition: interaction.hpp:334
void set_location(Vec newloc)
Definition: interaction.hpp:2004
FixedForceRegion(vector< FixedForceRegionAtom > atoms=vector< FixedForceRegionAtom >())
Definition: interaction.hpp:408
virtual Vec diff(Vec r1, Vec r2)=0
Distance between two points, given boundary conditions.
flt energy(Box &box)
Definition: interaction.hpp:361
vector< Charged > atoms
Definition: interaction.hpp:1953
Matrix stress(Box &box)
The force-moment tensor for the current simulation:
Definition: interaction.hpp:1938
Vec F
Definition: interaction.hpp:358
flt sigcut
Definition: interaction.hpp:885
PairList ignorepairs
Definition: interaction.hpp:1954
Definition: interaction.hpp:1512
vector< AngleGrouping > triples
Definition: interaction.hpp:652
Definition: interaction.hpp:1979
void apply_force(Box &box, Vec f, SpheroCylinderDiff diff, flt I)
Definition: interaction.hpp:2157
Vec forces(Box &box)
Definition: interaction.hpp:997
Vec f
forces
Definition: box.hpp:240
vector< flt > epsilons
Definition: interaction.hpp:1066
void add(vector< flt > nums, AtomID a1, AtomID a2, AtomID a3, AtomID a4)
Definition: interaction.hpp:692
LoisLinPair(LoisLinAtom a1, LoisLinAtom a2, flt eps, flt sig, flt f, flt l)
Definition: interaction.hpp:1479
flt sigcut
Definition: interaction.hpp:822
flt sigma
Definition: interaction.hpp:966
Definition: interaction.hpp:1388
flt sig
Definition: interaction.hpp:1330
flt q2
Definition: interaction.hpp:333
Definition: interaction.hpp:2145
EpsSigCutAtom(AtomID a, flt epsilon, flt sigma, flt cut)
Definition: interaction.hpp:803
flt c2
Definition: interaction.hpp:1165
flt sig
Definition: interaction.hpp:2164
virtual ~FPairXFunct()=0
Definition: interaction.hpp:745
vector< SoftWall * > walls
Definition: interaction.hpp:2044
IDPair & p2
Definition: interaction.hpp:2147
IEpsRepsSigExpCutAtom()
Definition: interaction.hpp:938
static flt energy(const flt rsig)
Definition: interaction.hpp:162
flt max_size()
Definition: interaction.hpp:1161
vector< BondGrouping > pairs
Definition: interaction.hpp:610
vector< P > & pair_iter()
Definition: interaction.hpp:1604
Vec get_axis()
Definition: interaction.hpp:2036
IEpsSigCutAtom()
Definition: interaction.hpp:886
static flt energy(const Vec diff, const flt eps, const flt sig)
Definition: interaction.hpp:154
RepulsionPair(EpsSigExpAtom a1, EpsSigExpAtom a2)
Definition: interaction.hpp:1259
LoisLinAtom(AtomID a, LoisLinAtom other)
Definition: interaction.hpp:1467
The basic Interaction class, used to represent a potential function.
Definition: interaction.hpp:58
AtomID atom2
Definition: interaction.hpp:762
void set_forces(Box &box)
Definition: interaction.hpp:422
uint size()
Definition: interaction.hpp:1548
sptr< SCBox > box
Definition: interaction.hpp:1523
flt xij
Definition: interaction.hpp:734
NListed(sptr< AtomVec > vec, sptr< NeighborList > neighbors)
Definition: interaction.hpp:1582
Vec direction
Definition: interaction.hpp:393
void ignore(AtomID a1, AtomID a2)
Definition: interaction.hpp:2220
Definition: interaction.hpp:1443
Definition: interaction.hpp:537
AtomID atom2
Definition: interaction.hpp:1258
flt energy(Box &box)
Definition: interaction.hpp:1114
flt get_sigma(IEpsISigExpAtom &other)
Definition: interaction.hpp:1235
EpsSigCutAtom()
Definition: interaction.hpp:802
void set_forces(Box &box)
Definition: interaction.hpp:1828
Vec get_normal()
Definition: interaction.hpp:2007
flt energy(Box &box)
Definition: interaction.hpp:1336
Definition: interaction.hpp:1329
Definition: interaction.hpp:676
Definition: interaction.hpp:225
Vec forces(Box &box)
Definition: interaction.hpp:1135
SCSpringList(SCAtomVec *scs, flt eps, flt sig, flt l)
Create an SCSpringList based on scs, using an epsilon of eps, a diameter of sigma, and a "length" of l, where l is cap center-to-cap center distance.
Definition: interaction.hpp:2201
AtomID first() const
Definition: box.hpp:288
flt pressure(Box &box)
Partial pressure due to this Interaction.
Definition: interaction.hpp:576
flt get_sigma(IEpsISigCutAtom &other)
Definition: interaction.hpp:838
AtomID get_id(cuint n)
Definition: box.hpp:442
virtual Atom & get(cuint n)
Definition: box.hpp:318
flt m2
Definition: interaction.hpp:485
virtual flt setForcesGetEnergy(Box &box)
Definition: interaction.hpp:1840
SCSpringList(SCAtomVec *scs, flt eps, flt sig, vector< flt > ls)
Definition: interaction.hpp:2203
IEpsISigExpAtom(AtomID a, IEpsISigCutAtom other)
Definition: interaction.hpp:1226
RandomForce(AtomGroup &agroup, flt force_mag, flt freq, RandomForceType force_type=UNIFORM)
Definition: interaction.hpp:557
EisMclachlanPair(EisMclachlanAtom a1, EisMclachlanAtom a2)
Definition: interaction.hpp:1169
vector< flt > epsilons
Definition: interaction.hpp:880
void fill(Box &box, ForcePairX &fpair)
Definition: interaction.hpp:1296
flt set_forces_get_pressure(Box &box)
Set forces (Atom.f) and return at the same time (see pressure()).
Definition: interaction.hpp:669
Definition: interaction.hpp:1201
virtual AtomID get_id(cuint n)=0
SoftWallCylinder(Vec loc, Vec axis, flt radius, flt expt=2.0)
Definition: interaction.hpp:2021
Definition: interaction.hpp:189
flt sigma
Definition: interaction.hpp:784
LJAttractCutPair(IEpsISigCutAtom a1, IEpsISigCutAtom a2)
Definition: interaction.hpp:910
flt cut_energy
Definition: interaction.hpp:232
void add(FixedSpringAtom a)
Definition: interaction.hpp:461
flt set_forces_get_pressure(Box &box)
Set forces (Atom.f) and return at the same time (see pressure()).
Definition: interaction.hpp:1669
AtomID atom2
Definition: interaction.hpp:1097
IEpsISigExpAtom(AtomID a, vector< flt > epsilons, vector< flt > sigmas, uint indx, flt exponent=2.5)
Definition: interaction.hpp:1220
virtual Matrix stress(Box &box)
The force-moment tensor for the current simulation:
Definition: interaction.hpp:93
flt sigma
Definition: interaction.hpp:1192
AtomID last() const
Definition: box.hpp:289
Definition: interaction.hpp:2191
Definition: interaction.hpp:541
AtomGroup * g2
Definition: interaction.hpp:483
AtomID atom2
Definition: interaction.hpp:904
void add(A atm)
Definition: interaction.hpp:1544
Vec forces(const Vec r)
Definition: interaction.hpp:341
A 3x3 matrix, with methods for adding, subtracting, dot product, etc.
Definition: vec.hpp:360
Definition: interaction.hpp:739
vector< A > group
Definition: interaction.hpp:1525
SimpleListed()
Definition: interaction.hpp:1543
flt sigcut
Definition: interaction.hpp:1473
IEpsRepsSigCutAtom(AtomID a, IEpsRepsSigCutAtom other)
Definition: interaction.hpp:1077
Vec diff(Vec r1, Vec r2)
Distance between two points, given boundary conditions.
Definition: interaction.hpp:2061
flt pressure(Box &box)
Partial pressure due to this Interaction.
Definition: interaction.hpp:388
Definition: interaction.hpp:588
virtual flt set_forces_get_pressure(Box &box)
Set forces (Atom.f) and return at the same time (see pressure()).
Definition: interaction.hpp:1628
flt x0
Definition: interaction.hpp:268
flt q1q2
Definition: interaction.hpp:761
uint indx
Definition: interaction.hpp:1216
Definition: interaction.hpp:732
bool same_atoms(AngleGrouping &other)
Definition: interaction.hpp:644
flt energy(Box &box, IDPair &pair)
Definition: interaction.hpp:1784
flt energy(const Vec &diff)
Definition: interaction.hpp:206
Definition: interaction.hpp:1015
Definition: interaction.hpp:902
static flt forces(const flt rsig, const flt cutsig)
Definition: interaction.hpp:255
flt exponent
Definition: interaction.hpp:932
void add(AtomID a, Vec loc, flt k, bool usex=true, bool usey=true, bool usez=true)
Definition: interaction.hpp:462
void set_forces(Box &box)
Definition: interaction.hpp:1719
flt set_forces_get_pressure(Box &box)
Set forces (Atom.f) and return at the same time (see pressure()).
Definition: interaction.hpp:580
flt energy(Box &box)
Potential energy due to this Interaction.
Definition: interaction.hpp:378
flt x0
Definition: interaction.hpp:484
uint ignore_size() const
Definition: interaction.hpp:1968
void add_pair(AtomID a1, AtomID a2)
Add a pair to the list.
Definition: trackers.hpp:61
AtomID a
Definition: interaction.hpp:359
static flt get_angle(const Vec &diff1, const Vec &diff2, const Vec &diff3)
Definition: interaction.cpp:293
Vec forces(Box &box)
Definition: interaction.hpp:1182
vector< DihedralGrouping > groups
Definition: interaction.hpp:688
flt sigma
Definition: interaction.hpp:884
vector< FixedForceRegionAtom > atoms
Definition: interaction.hpp:406
uint size() const
Definition: interaction.hpp:466
Vec forces(Box &box)
Definition: interaction.hpp:1272
flt lastf
Definition: interaction.hpp:1994
flt energy(const Vec &diff)
Definition: interaction.hpp:170
array< int, NDIM > fixed_box
Definition: interaction.hpp:597
flt sigma
Definition: interaction.hpp:1389
SCSpringPair(IDPair &p1, IDPair &p2, flt eps, flt sig, flt l)
Definition: interaction.hpp:2168
flt pressure(Box &box)
Partial pressure due to this Interaction.
Definition: interaction.hpp:720
flt set_forces_get_pressure(Box &box)
Set forces (Atom.f) and return at the same time (see pressure()).
Definition: interaction.hpp:721
Definition: interaction.hpp:965
flt get_sigma(IEpsRepsSigExpCutAtom &other)
Definition: interaction.hpp:958
uint size() const
Definition: interaction.hpp:670
flt pressure(Box &box)
Partial pressure due to this Interaction.
Definition: interaction.hpp:509
flt energy(Box &box)
Definition: interaction.hpp:1032
Vec get_location()
Definition: interaction.hpp:2034
EnergyForce get_energy_forces(Box &box)
Definition: interaction.hpp:1279
Definition: interaction.hpp:147
flt k
Definition: interaction.hpp:431
Definition: interaction.hpp:1538
vector< RandomForceAtom > group
Definition: interaction.hpp:553
Vec norm
Definition: interaction.hpp:1991
vector< flt > sigmas
Definition: interaction.hpp:1214
vector< A > & atom_list()
Definition: interaction.hpp:1553
EpsSigAtom(AtomID a, flt epsilon, flt sigma)
Definition: interaction.hpp:769
bool ywalls
Definition: interaction.hpp:2043
Definition: interaction.hpp:2113
bool same_atoms(BondGrouping &other)
Definition: interaction.hpp:602
#define NDIM
Definition: vecrand.hpp:12
IDPair pair(cuint n)
Definition: interaction.hpp:2130
Definition: box.hpp:282
AtomID a1
Definition: interaction.hpp:733
virtual flt set_forces_get_pressure(Box &box)
Set forces (Atom.f) and return at the same time (see pressure()).
Definition: interaction.hpp:69
flt get_epsilon(IEpsRepsSigExpCutAtom &other)
Definition: interaction.hpp:953
vector< FixedSpringAtom > atoms
Definition: interaction.hpp:458
AtomID atom2
Definition: interaction.hpp:785
bool usecos
Definition: interaction.hpp:281
flt max_size()
Definition: interaction.hpp:1469
BondDiffType diff_type
Definition: interaction.hpp:596
virtual ~InteractionPairsX()
Definition: interaction.hpp:751
LJRepulsePair(EpsSigAtom LJ1, EpsSigAtom LJ2)
Definition: interaction.hpp:786
EpsSigAtom()
Definition: interaction.hpp:768
flt expt
Definition: interaction.hpp:1992
Definition: interaction.hpp:747
AngleGrouping(flt k, flt x0, AtomID a1, AtomID a2, AtomID a3)
Definition: interaction.hpp:642
LoisLinPairMin(LoisLinAtom a1, LoisLinAtom a2)
Definition: interaction.hpp:1513
flt energy(Box &box)
Definition: interaction.hpp:440
Definition: interaction.hpp:1521
flt pressure(Box &box)
Partial pressure due to this Interaction.
Definition: interaction.hpp:477
AtomID atom2
Definition: interaction.hpp:1400
Definition: interaction.hpp:2162
Definition: interaction.hpp:2140
void add(Charged a)
Definition: interaction.hpp:1963
uint size() const
Definition: interaction.hpp:414
flt springk
Definition: interaction.hpp:279
Vec forces_pair(P pair, Box &box)
Definition: interaction.hpp:1611
vector< P > pairs
Definition: interaction.hpp:1578
Definition: interaction.hpp:357
Definition: interaction.hpp:551
Definition: interaction.hpp:727
flt energy(const Vec &diff)
Definition: interaction.hpp:240
NListedVirial(sptr< AtomVec > vec, sptr< NeighborList > neighbors)
Definition: interaction.hpp:1624
sptr< NeighborList > neighbor_list()
Definition: interaction.hpp:1613
flt set_forces_get_pressure(Box &box)
Set forces (Atom.f) and return at the same time (see pressure()).
Definition: interaction.hpp:1890
Repulsive LJ: .
Definition: interaction.hpp:783
virtual void set_forces(Box &box)=0
Definition: interaction.hpp:879
Vec forces(const Vec &diff)
Definition: interaction.hpp:262
SCBoxed(sptr< AtomVec > atomv, sptr< SCBox > box)
Definition: interaction.hpp:1527
flt sigma
Definition: interaction.hpp:931
RandomForceType
Definition: interaction.hpp:535
int get_fixed(uint i)
Definition: interaction.hpp:601
flt sigma
Definition: interaction.hpp:767
EpsSigExpAtom(AtomID a, EpsSigExpAtom other)
Definition: interaction.hpp:1196
flt energy(Box &box)
Potential energy due to this Interaction.
Definition: interaction.hpp:1706
vector< A > & atom_list()
Definition: interaction.hpp:1612
flt get_epsilon(IEpsISigCutAtom &other)
Definition: interaction.hpp:831
Definition: interaction.hpp:276
flt get_epsilon(IEpsSigCutAtom &other)
Definition: interaction.hpp:892
flt get_last_f()
Definition: interaction.hpp:2009
virtual ~Interaction()
Definition: interaction.hpp:98
Vec loc
Definition: interaction.hpp:1990
flt sig
Definition: interaction.hpp:1094
The virtual interface for the shape of the space and its boundaries.
Definition: box.hpp:50
bool usepow
Definition: interaction.hpp:308
Definition: interaction.hpp:754
Definition: interaction.hpp:1620
flt sigma
Definition: interaction.hpp:1981
flt costheta
Definition: interaction.hpp:300
ChargePair(Charged a1, Charged a2)
Definition: interaction.hpp:763
LJishPair(IEpsRepsSigExpCutAtom LJ1, IEpsRepsSigExpCutAtom LJ2)
Definition: interaction.hpp:968
LoisOhernAtom()
Definition: interaction.hpp:1390
Eigen::Matrix< flt, NDIM, 1 > Vec
The basic physics vector.
Definition: vecrand.hpp:53
flt cut_energy
Definition: interaction.hpp:194
FixedSpring(vector< FixedSpringAtom > atoms=vector< FixedSpringAtom >())
Definition: interaction.hpp:460
Definition: interaction.hpp:2012
Definition: interaction.hpp:429
FixedForce(vector< FixedForceAtom > atoms=vector< FixedForceAtom >())
Definition: interaction.hpp:369
static Vec forces(const Vec diff, const flt eps, const flt sig, const flt cutsig)
Definition: interaction.hpp:211
~BondAngle()
Definition: interaction.hpp:294
Vec forces(Box &box)
Definition: interaction.hpp:1424
static Vec forces(const Vec diff, const flt eps, const flt sig)
Definition: interaction.hpp:130
flt energy(Box &box, SpheroCylinderDiff diff)
Definition: interaction.hpp:2174
virtual flt energy(Box &box)
Potential energy due to this Interaction.
Definition: interaction.hpp:1630
EpsSigAtom(AtomID a, EpsSigAtom other)
Definition: interaction.hpp:771
bool add(flt k, flt x0, AtomID a1, AtomID a2, AtomID a3, bool replace=true)
Definition: interaction.hpp:660
flt lambda2
Definition: interaction.hpp:2142
Definition: interaction.hpp:104
Vec rand_vec_boxed()
Generate a random vector inside a box with sides of length 1.
Definition: vecrand.cpp:19
LJAttractFixedRepulsePair(IEpsRepsSigCutAtom a1, IEpsRepsSigCutAtom a2)
Definition: interaction.hpp:1099
void add(WallAtom a)
Definition: interaction.hpp:2023
EpsSigExpDragAtom(AtomID a, EpsSigExpDragAtom other)
Definition: interaction.hpp:1324
flt energy(Box &box)
Definition: interaction.hpp:1265
uint size() const
for iterating over neighbors
Definition: trackers.hpp:91
Definition: interaction.hpp:330
flt max_size()
Definition: interaction.hpp:1085
COMSpring(AtomGroup *g1, AtomGroup *g2, flt k, flt x0=0)
Definition: interaction.hpp:487
vector< flt > epsilons
Definition: interaction.hpp:1213
IEpsRepsSigCutAtom(AtomID a, vector< flt > epsilons, flt repeps, flt sigma, uint indx, flt cut)
Definition: interaction.hpp:1071
flt q
Definition: interaction.hpp:755
flt screen
Definition: interaction.hpp:1955
flt sigma
Definition: interaction.hpp:231
flt energy(Box &box)
Potential energy due to this Interaction.
Definition: interaction.hpp:416
flt sigmai
Definition: interaction.hpp:1155
LoisLinAtom(AtomID a, flt eps, flt sigma, flt depth, flt width)
Definition: interaction.hpp:1465
P get_pair(IDPair &pair)
Definition: interaction.hpp:1593
flt x0
Definition: interaction.hpp:594
vector< A > atoms
Definition: interaction.hpp:1540
DihedralGrouping(vector< flt > cos_coefficients, vector< flt > sincoeffs, AtomID a1, AtomID a2, AtomID a3, AtomID a4, bool usepow=true)
Definition: interaction.hpp:680
flt lastf
Definition: interaction.hpp:2018
flt pressure(Box &box)
Partial pressure due to this Interaction.
Definition: interaction.hpp:1752
EpsSigCutAtom(AtomID a, EpsSigCutAtom other)
Definition: interaction.hpp:805
flt energy(Box &box)
Potential energy due to this Interaction.
Definition: interaction.hpp:489
WalledBox2D(Vec size, bool xwalled, bool ywalled, flt expt=2.0)
Definition: interaction.hpp:2046
Definition: interaction.hpp:800
Definition: interaction.hpp:1154
const flt OVERNDIM
A constant equal to , where is the number of dimensions.
Definition: vecrand.hpp:66
AtomID a2
Definition: interaction.hpp:733
tuple p
Definition: namereplacer.py:23
Vec forces(Box &box)
Definition: interaction.hpp:792
Definition: interaction.hpp:607
A pointer to an Atom.
Definition: box.hpp:247
AtomID a1
Definition: interaction.hpp:641
uint size() const
Definition: interaction.hpp:563
flt sigcut
Definition: interaction.hpp:937
AtomGroup * g1
Definition: interaction.hpp:482
flt max_size()
Definition: interaction.hpp:1326
Definition: interaction.hpp:1951
flt energy(Box &box)
Definition: interaction.hpp:1408
Truncated and shifted Lennard-Jones, in the form .
Definition: interaction.hpp:861
flt sig
Definition: interaction.hpp:1067
unsigned long long contacts(Box &box)
number of Atom pairs with E != 0
Definition: interaction.hpp:1792
IEpsSigCutAtom(AtomID a, vector< flt > epsilons, uint indx, flt sigma, flt cut)
Definition: interaction.hpp:887
void set_forces(Box &box)
Definition: interaction.hpp:1653
Vec forces(Box &box)
Definition: interaction.hpp:1047
Spring(const flt k, const flt x0)
Definition: interaction.hpp:270
Definition: interaction.hpp:1988
AtomID a3
Definition: interaction.hpp:641
const unsigned int cuint
Definition: box.hpp:13
SCAtomVec(uint N, flt mass)
Definition: interaction.hpp:2125
EpsSigExpDragAtom(AtomID a, flt eps, flt sigma, flt gamma, flt exponent=2.5)
Definition: interaction.hpp:1322
vector< flt > Fs
Definition: interaction.hpp:397
Vec diff(Vec r1, Vec r2, array< int, NDIM > boxes)
Definition: interaction.hpp:2067
EisMclachlanAtom(AtomID a, flt dist, flt sigmai)
Definition: interaction.hpp:1157