Path Integral Quantum Monte Carlo
potential.cpp
Go to the documentation of this file.
1 
8 #include "potential.h"
9 #include "path.h"
10 #include "lookuptable.h"
11 #include "communicator.h"
12 
13 #include <boost/math/special_functions/ellint_1.hpp>
14 #include <boost/math/special_functions/ellint_2.hpp>
15 
16 // ---------------------------------------------------------------------------
17 // ---------------------------------------------------------------------------
18 // POTENTIAL BASE CLASS ------------------------------------------------------
19 // ---------------------------------------------------------------------------
20 // ---------------------------------------------------------------------------
21 
22 /**************************************************************************/
26 
27 }
28 
29 /**************************************************************************/
33 }
34 
35 /**************************************************************************/
41 Array<dVec,1> PotentialBase::initialConfig(const Container *boxPtr, MTRand &random,
42  const int numParticles) {
43 
44  /* The particle configuration */
45  Array<dVec,1> initialPos(numParticles);
46  initialPos = 0.0;
47 
48  /* Get the linear size per particle, and the number of particles */
49  double initSide = pow((1.0*numParticles/boxPtr->volume),-1.0/(1.0*NDIM));
50 
51  /* We determine the number of initial grid boxes there are in
52  * in each dimension and compute their size */
53  int totNumGridBoxes = 1;
54  iVec numNNGrid;
55  dVec sizeNNGrid;
56 
57  for (int i = 0; i < NDIM; i++) {
58  numNNGrid[i] = static_cast<int>(ceil((boxPtr->side[i] / initSide) - EPS));
59 
60  /* Make sure we have at least one grid box */
61  if (numNNGrid[i] < 1)
62  numNNGrid[i] = 1;
63 
64  /* Compute the actual size of the grid */
65  sizeNNGrid[i] = boxPtr->side[i] / (1.0 * numNNGrid[i]);
66 
67  /* Determine the total number of grid boxes */
68  totNumGridBoxes *= numNNGrid[i];
69  }
70 
71  /* Now, we place the particles at the middle of each box */
72  PIMC_ASSERT(totNumGridBoxes>=numParticles);
73  dVec pos;
74  for (int n = 0; n < totNumGridBoxes; n++) {
75 
76  iVec gridIndex;
77  for (int i = 0; i < NDIM; i++) {
78  int scale = 1;
79  for (int j = i+1; j < NDIM; j++)
80  scale *= numNNGrid[j];
81  gridIndex[i] = (n/scale) % numNNGrid[i];
82  }
83 
84  for (int i = 0; i < NDIM; i++)
85  pos[i] = (gridIndex[i]+0.5)*sizeNNGrid[i] - 0.5*boxPtr->side[i];
86 
87  boxPtr->putInside(pos);
88 
89  if (n < numParticles)
90  initialPos(n) = pos;
91  else
92  break;
93  }
94 
95  return initialPos;
96 }
97 
98 /**************************************************************************/
106 void PotentialBase::output(const double maxSep) {
107  dVec sep;
108  sep = 0.0;
109  for (double d = 0; d < maxSep; d+= (maxSep/1.0E6)) {
110  sep[0] = d;
111  communicate()->file("debug")->stream()
112  << format("%16.12E\t%16.12E\n") % d % V(sep);
113  }
114 }
115 
116 /**************************************************************************/
119 double PotentialBase::deltaSeparation(double sep1, double sep2) const {
120 
121  double delta = sep2-sep1;
122  while (delta >= 0.5*constants()->L())
123  delta -= constants()->L();
124  while (delta < -0.5*constants()->L())
125  delta += constants()->L();
126 
127  return delta;
128 }
129 
130 /**************************************************************************/
135 Array<double,1> PotentialBase::getExcLen() {
136  /* The particle configuration */
137  Array<double,1> excLens(0);
138  return excLens;
139 }
140 // ---------------------------------------------------------------------------
141 // ---------------------------------------------------------------------------
142 // TABULATED POTENTIAL CLASS -------------------------------------------------
143 // ---------------------------------------------------------------------------
144 // ---------------------------------------------------------------------------
145 
146 /**************************************************************************/
150  extV = 0.0;
151  extdVdr = 0.0;
152  extd2Vdr2 = 0.0;
153 }
154 
155 /**************************************************************************/
159  lookupV.free();
160  lookupdVdr.free();
161  lookupd2Vdr2.free();
162 }
163 
164 /**************************************************************************/
168 void TabulatedPotential::initLookupTable(const double _dr, const double maxSep) {
169 
170  /* We now calculate the lookup tables for the interaction potential and
171  * its first and second derivatives. */
172  dr = _dr;
173  tableLength = int(maxSep/dr);
174  lookupV.resize(tableLength);
175  lookupdVdr.resize(tableLength);
176  lookupd2Vdr2.resize(tableLength);
177  lookupV = 0.0;
178  lookupdVdr = 0.0;
179  lookupd2Vdr2 = 0.0;
180 
181  double r = 0;
182 
183  for (int n = 0; n < tableLength; n++) {
184  lookupV(n) = valueV(r);
185  lookupdVdr(n) = valuedVdr(r);
186  lookupd2Vdr2(n) = valued2Vdr2(r);
187  r += dr;
188  }
189 
190  /* r = 0.0; */
191  /* for (int n = 0; n < tableLength; n++) { */
192  /* communicate()->file("debug")->stream() << format("%24.16e %24.16e\n") */
193  /* % r % lookupV(n); */
194  /* r += dr; */
195  /* }; */
196 
197  /* exit(-1); */
198 
199 
200 // cout << format("%16.8E%16.8E%16.8E%16.8E%16.8E%16.8E%16.8E\n") % r % lookupV(n) % valueV(r) %
201 
202 // double rc = constants()->rc();
203 // for (int n = 0; n < tableLength; n++) {
204 // r += dr;
205 // if (r <= rc) {
206 // lookupV(n) = valueV(r) - valueV(rc) - valuedVdr(rc)*(r-rc);
207 // lookupdVdr(n) = valuedVdr(r) - valuedVdr(rc);
208 // }
209 // else {
210 // lookupV(n) = 0.0;
211 // lookupdVdr(n) = 0.0;
212 // }
213 // cout << format("%16.8E%16.8E%16.8E%16.8E%16.8E%16.8E%16.8E\n") % r % lookupV(n) % valueV(r) %
214 // lookupdVdr(n) % valuedVdr(r) % (lookupV(n) - valueV(r)) % (lookupdVdr(n) - valuedVdr(r));
215 // }
216 
217 }
218 
219 /**************************************************************************/
226 double TabulatedPotential::newtonGregory(const Array<double,1> &VTable,
227  const TinyVector<double,2> &extVal, const double r) {
228 
229  double rdr = r/dr;
230  int k = int(rdr);
231 
232  if (k <= 0)
233  return extVal[0];
234 
235  if (k >= tableLength)
236  return extVal[1];
237 
238  double xi = rdr - 1.0*k;
239  double vkm1 = VTable(k-1);
240  double vk = VTable(k);
241  double vkp1 = VTable(k+1);
242 
243  double T1 = vkm1 + (vk - vkm1) * xi;
244  double T2 = vk + (vkp1 - vk) * (xi - 1.0);
245 
246  return (T1 + 0.5 * (T2 - T1) * xi);
247 }
248 
249 /**************************************************************************/
255 double TabulatedPotential::direct(const Array<double,1> &VTable,
256  const TinyVector<double,2> &extVal, const double r) {
257 
258  int k = int(r/dr);
259  if (k <= 0)
260  return extVal[0];
261 
262  if (k >= tableLength)
263  return extVal[1];
264 
265  return VTable(k);
266 }
267 
268 // ---------------------------------------------------------------------------
269 // ---------------------------------------------------------------------------
270 // FREE POTENTIAL CLASS ------------------------------------------------------
271 // ---------------------------------------------------------------------------
272 // ---------------------------------------------------------------------------
273 
274 /**************************************************************************/
278 }
279 
280 /**************************************************************************/
284 }
285 
286 // ---------------------------------------------------------------------------
287 // ---------------------------------------------------------------------------
288 // SINGLE WELL POTENTIAL CLASS -----------------------------------------------
289 // ---------------------------------------------------------------------------
290 // ---------------------------------------------------------------------------
291 
292 /**************************************************************************/
296 }
297 
298 /**************************************************************************/
302 }
303 
304 // ---------------------------------------------------------------------------
305 // ---------------------------------------------------------------------------
306 // HARMONIC POTENTIAL CLASS --------------------------------------------------
307 // ---------------------------------------------------------------------------
308 // ---------------------------------------------------------------------------
309 
310 /**************************************************************************/
313 HarmonicPotential::HarmonicPotential(double _omega) : PotentialBase() {
314  omega2 = _omega*_omega;
315 }
316 
317 HarmonicPotential::HarmonicPotential() : PotentialBase() {
318  omega2 = 1.0;
319 }
320 
321 /**************************************************************************/
325 }
326 
327 /**************************************************************************/
332 Array<dVec,1> HarmonicPotential::initialConfig(const Container *boxPtr, MTRand &random,
333  const int numParticles) {
334 
335  /* The particle configuration */
336  Array<dVec,1> initialPos(numParticles);
337  initialPos = 0.0;
338 
339  for (int n = 0; n < numParticles; n++) {
340  for (int i = 0; i < NDIM; i++)
341  initialPos(n)[i] = 0.1*boxPtr->side[i]*(-1.0 + 2.0*random.rand());
342  }
343 
344  return initialPos;
345 }
346 
347 // ---------------------------------------------------------------------------
348 // ---------------------------------------------------------------------------
349 // HARMONIC TUBE POTENTIAL CLASS ---------------------------------------------
350 // ---------------------------------------------------------------------------
351 // ---------------------------------------------------------------------------
352 
353 /**************************************************************************/
359 {
360  /* c is a dimensionless constant */
361  c = 1.20272;
362 
363  /* We have to determine the frequency of the oscillator from it's length.
364  * w = \hbar / (m R^2). It is measured in THz */
365  w = 6.35077 / (radius*radius*constants()->m());
366 }
367 
368 /**************************************************************************/
372 }
373 
374 // ---------------------------------------------------------------------------
375 // ---------------------------------------------------------------------------
376 // DELTA POTENTIAL CLASS -----------------------------------------------------
377 // ---------------------------------------------------------------------------
378 // ---------------------------------------------------------------------------
379 
380 /**************************************************************************/
387 DeltaPotential::DeltaPotential(double _sigma, double _g) : PotentialBase()
388 {
389  /* Define the parameters of the delta function potential. */
390  norm = _g/sqrt(2.0*_sigma*_sigma*M_PI);
391  inv2sigma2 = 1.0/(2.0*_sigma*_sigma);
392 }
393 
394 /**************************************************************************/
398 }
399 
400 // ---------------------------------------------------------------------------
401 // ---------------------------------------------------------------------------
402 // LORENTZIAN POTENTIAL CLASS ------------------------------------------------
403 // ---------------------------------------------------------------------------
404 // ---------------------------------------------------------------------------
405 
406 /**************************************************************************/
414 {
415  /* Define the parameters of the Lorentzian delta function potential. */
416  a = _a;
417  c = _c;
418  norm = 2.0 * c * a / M_PI;
419 }
420 
421 /**************************************************************************/
425 }
426 
427 // ---------------------------------------------------------------------------
428 // ---------------------------------------------------------------------------
429 // SUTHERLAND POTENTIAL CLASS ------------------------------------------------
430 // ---------------------------------------------------------------------------
431 // ---------------------------------------------------------------------------
432 
433 /**************************************************************************/
440 {
441  g = 2.0*constants()->lambda() * _g * (_g - 1.0);
442  pioL = M_PI / constants()->L();
443 }
444 
445 /**************************************************************************/
449  // empty
450 }
451 
452 // ---------------------------------------------------------------------------
453 // ---------------------------------------------------------------------------
454 // DIPOLE POTENTIAL CLASS ----------------------------------------------------
455 // ---------------------------------------------------------------------------
456 // ---------------------------------------------------------------------------
457 
458 /**************************************************************************/
463 {
464  /* if (NDIM==2) { */
465  /* int N = constants()->initialNumParticles(); */
466  /* tailV = M_PI * N * N / (constants()->L() * constants()->L() * constants()->rc()); */
467  /* } */
468 }
469 
470 /**************************************************************************/
474  // empty
475 }
476 
477 // ---------------------------------------------------------------------------
478 // ---------------------------------------------------------------------------
479 // FIXED AZIZ POTENTIAL CLASS ------------------------------------------------
480 // ---------------------------------------------------------------------------
481 // ---------------------------------------------------------------------------
482 
483 /**************************************************************************/
491  aziz(_boxPtr) {
492 
493  char state; // Fixed or updateable?
494  dVec pos; // The loaded position
495 
496  /* Initialize the cutoff^2 */
497  rc2 = constants()->rc2();
498 
499  /* We start with an array of size 500 */
500  fixedParticles.resize(500);
501 
502  /* Here we load both the number and location of fixed helium atoms from disk. */
503  numFixedParticles = 0;
504  int n = 0;
505  while (!communicate()->file("fixed")->stream().eof()) {
506  if (communicate()->file("fixed")->stream().peek() == '#') {
507  communicate()->file("fixed")->stream().ignore(512,'\n');
508  }
509  else {
510  communicate()->file("fixed")->stream() >> state;
511  for (int i = 0; i < NDIM; i++)
512  communicate()->file("fixed")->stream() >> pos[i];
513 
514  /* If the particle is labelled with an 'F' it is fixed and should
515  * be included here */
516  if (state == 'F') {
517  numFixedParticles++;
518  if (numFixedParticles >= int(fixedParticles.size()))
519  fixedParticles.resizeAndPreserve(numFixedParticles);
520 
521  /* Put the initial position in the container */
522  _boxPtr->putInside(pos);
523  fixedParticles(n) = pos;
524  n++;
525  }
526  communicate()->file("fixed")->stream().ignore();
527  }
528  }
529 
530  fixedParticles.resizeAndPreserve(numFixedParticles);
531 
532  /* Now that we have the particle positions, create a new lookup table pointer
533  * and initialize it */
534  lookupPtr = new LookupTable(_boxPtr,1,numFixedParticles);
535  lookupPtr->updateGrid(fixedParticles);
536 
537  /* Resize and initialize our local grid box arrays */
538  fixedBeadsInGrid.resize(lookupPtr->getTotNumGridBoxes(),numFixedParticles);
539  numFixedBeadsInGrid.resize(lookupPtr->getTotNumGridBoxes());
540  fixedBeadsInGrid = XXX;
541  numFixedBeadsInGrid = 0;
542 
543  /* Create a local copy of all beads in each grid box plus nearest neighbors.
544  * This will drastically speed up the computing of potential energies. */
545  for (n = 0; n < lookupPtr->getTotNumGridBoxes(); n++) {
546  lookupPtr->updateFullInteractionList(n,0);
547  numFixedBeadsInGrid(n) = lookupPtr->fullNumBeads;
548  for (int m = 0; m < lookupPtr->fullNumBeads; m++)
549  fixedBeadsInGrid(n,m) = lookupPtr->fullBeadList(m)[1];
550  }
551 
552 }
553 
554 /**************************************************************************/
558  delete lookupPtr;
559  fixedParticles.free();
560  fixedBeadsInGrid.free();
561  numFixedBeadsInGrid.free();
562 }
563 
564 /**************************************************************************/
568 double FixedAzizPotential::V(const dVec &pos) {
569 
570  double totV = 0.0;
571 
572  /* We first find the grid box number that the particle resides in */
573  int gridNumber = lookupPtr->gridNumber(pos);
574 
575  /* We now loop over all fixed particles in this grid box, only computing
576  * interactions when the separation is less than the cutoff */
577  dVec sep;
578  for (int n = 0; n < numFixedBeadsInGrid(gridNumber); n++) {
579  sep = fixedParticles(fixedBeadsInGrid(gridNumber,n)) - pos;
580  lookupPtr->boxPtr->putInBC(sep);
581  if (dot(sep,sep) < rc2)
582  totV += aziz.V(sep);
583  }
584 
585  return totV;
586 }
587 
588 /**************************************************************************/
593 
594  dVec totGradV;
595  totGradV = 0.0;
596 
597  /* We first find the grid box number that the particle resides in */
598  int gridNumber = lookupPtr->gridNumber(pos);
599 
600  /* We now loop over all fixed particles in this grid box, only computing
601  * the gradient of interactions when the separation is less than the cutoff */
602  dVec sep;
603  for (int n = 0; n < numFixedBeadsInGrid(gridNumber); n++) {
604  sep = fixedParticles(fixedBeadsInGrid(gridNumber,n)) - pos;
605  lookupPtr->boxPtr->putInBC(sep);
606  if (dot(sep,sep) < rc2)
607  totGradV += aziz.gradV(sep);
608  }
609 
610  return totGradV;
611 }
612 
613 /**************************************************************************/
620 Array<dVec,1> FixedAzizPotential::initialConfig(const Container *boxPtr, MTRand &random,
621  const int numParticles) {
622 
623  /* The particle configuration */
624  Array<dVec,1> initialPos(1);
625  initialPos = 0.0;
626 
627  int locNumParticles = 0; // Number of updateable particles
628  char state; // Update or Fix
629  dVec pos; // The current position
630 
631  /* We go through all lines in the fixed input file, discarding any comments
632  * and assign the initial positions of the particles */
633  int n = 0;
634  while (!communicate()->file("fixed")->stream().eof()) {
635  if (communicate()->file("fixed")->stream().peek() == '#') {
636  communicate()->file("fixed")->stream().ignore(512,'\n');
637  }
638  else {
639  communicate()->file("fixed")->stream() >> state;
640  for (int i = 0; i < NDIM; i++)
641  communicate()->file("fixed")->stream() >> pos[i];
642 
643  /* If the particle is labelled with an 'U' it is updatable and should
644  * be included */
645  if (state == 'U') {
646  locNumParticles++;
647  initialPos.resizeAndPreserve(locNumParticles);
648 
649  /* Put the initial position in the box */
650  boxPtr->putInside(pos);
651 
652  /* Assign the position to all time slices*/
653  initialPos(n) = pos;
654  n++;
655  }
656  communicate()->file("fixed")->stream().ignore();
657  }
658  }
659 
660  /* Reset the file pointer */
661  communicate()->file("fixed")->stream().seekg(0, ios::beg);
662 
663  /* Return the initial Positions */
664  return initialPos;
665 }
666 
667 // ---------------------------------------------------------------------------
668 // ---------------------------------------------------------------------------
669 // FixedPositionLJPotential Class---------------------------------------------
670 // ---------------------------------------------------------------------------
671 // ---------------------------------------------------------------------------
672 
673 /**************************************************************************/
676 FixedPositionLJPotential::FixedPositionLJPotential (double _sigma, double _epsilon,
677  const Container *_boxPtr) : PotentialBase() {
678 
679  boxPtr = _boxPtr;
680  sigma = _sigma;
681  epsilon = _epsilon;
682 
683  Lz = boxPtr->side[NDIM-1];
684 
685  /* Fixed positions of FILENAME */
686  dVec pos; // The loaded position
687 
688  /* We start with an array of size 500 */
689  fixedParticles.resize(500);
690 
691  /* Here we load both the number and location of fixed positions from disk. */
692  numFixedParticles = 0;
693  int n = 0;
694  while (!communicate()->file("fixed")->stream().eof()) {
695  if (communicate()->file("fixed")->stream().peek() == '#') {
696  communicate()->file("fixed")->stream().ignore(512,'\n');
697  }
698  else {
699  for (int i = 0; i < NDIM; i++)
700  communicate()->file("fixed")->stream() >> pos[i];
701  numFixedParticles++;
702  if (numFixedParticles >= int(fixedParticles.size()))
703  fixedParticles.resizeAndPreserve(numFixedParticles);
704  fixedParticles(n) = pos;
705  n++;
706  communicate()->file("fixed")->stream().ignore();
707  }
708  }
709  fixedParticles.resizeAndPreserve(numFixedParticles);
710 
711  /* print out the potential to disk */
712  /* int numPoints = 200; */
713  /* double dx = _boxPtr->side[0]/numPoints; */
714  /* double dy = _boxPtr->side[1]/numPoints; */
715  /* pos[2] = -0.5*_boxPtr->side[2] + 2.635; */
716  /* for (int i = 0; i < numPoints; i++) { */
717  /* pos[0] = -0.5*_boxPtr->side[0] + i*dx; */
718  /* for (int j = 0; j < numPoints; j++) { */
719  /* pos[1] = -0.5*_boxPtr->side[1] + j*dy; */
720  /* communicate()->file("debug")->stream() << format("%24.16e %24.16e %24.16e\n") % pos[0] % pos[1] % V(pos); */
721  /* } */
722  /* } */
723 
724  /* exit(-1); */
725 }
726 
727 
728 /**************************************************************************/
732  fixedParticles.free();
733 }
734 
735 /**************************************************************************/
743 
744  /* Notes: for now I hard-code the potential at 1.5 \AA and a LJ-cutoff of
745  * 20 \AA */
746 
747  if (r[NDIM-1] < (-0.5*Lz + 1.5) )
748  return 87292.0;
749 
750  else if (r[NDIM-1] > 0.0)
751  return 0.0;
752 
753  double v = 0.0;
754  double sor = 0.0;
755  double x = 0.0;
756  dVec sep;
757  for (int i = 0; i < numFixedParticles; i++) {
758  sep[0] = fixedParticles(i)[0] - r[0];
759  sep[1] = fixedParticles(i)[1] - r[1];
760  boxPtr->putInBC(sep);
761  sep[2] = fixedParticles(i)[2] - r[2];
762  x = sqrt(dot(sep,sep));
763  if (x < 20.0) {
764  sor = sigma/x;
765  v += pow(sor,12)-pow(sor,6);
766  }
767  }
768  return 4.0*epsilon*v;
769 }
770 
771 // ---------------------------------------------------------------------------
772 // ---------------------------------------------------------------------------
773 // HARD CYLINDER POTENTIAL CLASS ---------------------------------------------
774 // ---------------------------------------------------------------------------
775 // ---------------------------------------------------------------------------
776 
777 /**************************************************************************/
782  PotentialBase(),
783  R(radius) {
784 }
785 
786 /**************************************************************************/
790 }
791 
792 /**************************************************************************/
800 PlatedLJCylinderPotential::PlatedLJCylinderPotential(const double Ro_, const double Rw,
801  const double sigmaPlated_, const double epsilonPlated_, const double densityPlated_) :
802  PotentialBase(),
804 {
805  /* The outer radius of the tube */
806  Ro = Ro_;
807 
808  /* Plating substrates */
809  densityPlated = densityPlated_; //0.207; // atoms / angstrom^3
810  epsilonPlated = epsilonPlated_; //36.13613; // Kelvin
811  sigmaPlated = sigmaPlated_; //3.0225; // angstroms
812  Ri = Ro-Rw; //3.5;
813 
814  // Neon Values
815  /* densityPlated = 0.207; // atoms / angstrom^3 FIXME */
816  /* epsilonPlated = 20.161242; // Kelvin */
817  /* sigmaPlated = 2.711; // angstroms */
818  /* Ri = Ro-(1.52*2); //FIXME Fix hardcoded R2 */
819 
820  /* Substrate parameters for MCM-41 obtained by fitting
821  * @see https://nano.delmaestro.org/index.php?title=Effective_external_potential_(nanopores)
822  */
823  density = 1.000; // atoms / angstrom^3
824  epsilon = 1.59; // Kelvin
825  sigma = 3.44; // angstroms
826 
827  /* We choose a mesh consisting of 10^6 points, and create the lookup table */
828  dR = (1.0E-6)*Ri;
829  initLookupTable(dR,Ri);
830 
831  /* Find the minimun of the potential */
832  minV = 1.0E5;
833  for (int n = 0; n < tableLength; n++) {
834  if (lookupV(n) < minV)
835  minV = lookupV(n);
836  }
837 
838  /* The extremal values for the lookup table */
839  extV = valueV(0.0),valueV(Ri);
840  extdVdr = valuedVdr(0.0),valuedVdr(Ri);
841 }
842 
843 /**************************************************************************/
847 }
848 
849 /**************************************************************************/
855 double PlatedLJCylinderPotential::V_(const double r, const double R,
856  const double sig, const double eps, const double dens)
857 {
858  double x = (r - (r>=R)*EPS) / R;
859 
860  double x2 = x*x;
861  double x4 = x2*x2;
862  double x6 = x2*x4;
863  double x8 = x4*x4;
864  double f1 = 1.0 / (1.0 - x2);
865  double sigoR3 = pow(sig/R,3.0);
866  double sigoR9 = sigoR3*sigoR3*sigoR3;
867 
868  double Kx2 = boost::math::ellint_1(x);
869  double Ex2 = boost::math::ellint_2(x);
870 
871  double v9 = (1.0*pow(f1,9.0)/(240.0)) * (
872  (1091.0 + 11156*x2 + 16434*x4 + 4052*x6 + 35*x8)*Ex2 -
873  8.0*(1.0 - x2)*(1.0 + 7*x2)*(97.0 + 134*x2 + 25*x4)*Kx2);
874  double v3 = 2.0*pow(f1,3.0) * ((7.0 + x2)*Ex2 - 4.0*(1.0-x2)*Kx2);
875 
876  return ((M_PI*eps*sig*sig*sig*dens/3.0)*(sigoR9*v9 - sigoR3*v3));
877 }
878 
879 /**************************************************************************/
884 double PlatedLJCylinderPotential::dVdr_(const double r, const double R,
885  const double sig, const double eps, const double dens)
886 {
887  double x = (r - (r>=R)*EPS) / R;
888 
889  /* dV/dr */
890  if (x < EPS)
891  return (1.28121E8/pow(R,11.0) - 102245.0/pow(R,5.0))*x;
892  else {
893  double x2 = x*x;
894  double x4 = x2*x2;
895  double x6 = x2*x4;
896  double x8 = x4*x4;
897  double f1 = 1.0 / (1.0 - x2);
898  double sigoR3 = pow(sig/R,3.0);
899  double sigoR9 = sigoR3*sigoR3*sigoR3;
900 
901  double Kx2 = boost::math::ellint_1(x);
902  double Ex2 = boost::math::ellint_2(x);
903 
904  double dv9dx =(3.0*pow(f1,10.0)/(80.0*x)) *
905  ( (1.0 + x2)*(35.0 + 5108*x2 + 22482*x4 + 5108*x6 + 35*x8)*Ex2 -
906  (1.0 - x2)*(35.0 + 3428*x2 + 15234*x4 + 12356*x6 +1715*x8)*Kx2 );
907  double dv3dx = (6.0*pow(f1,4.0)/x) *
908  ( (1.0 + 14*x2 + x4)*Ex2 - (1.0 + 6*x2 - 7*x4)*Kx2 );
909  return ((M_PI*eps*sig*sig*sig*dens/(3.0*R))*(sigoR9*dv9dx - sigoR3*dv3dx));
910  }
911 }
912 
913 /**************************************************************************/
917 double PlatedLJCylinderPotential::valueV(const double r) {
918 
919  return V_(r,Ri,sigmaPlated,epsilonPlated,densityPlated) -
920  V_(r,Ro,sigmaPlated,epsilonPlated,densityPlated) +
921  V_(r,Ro,sigma,epsilon,density);
922 }
923 
924 /**************************************************************************/
929 double PlatedLJCylinderPotential::valuedVdr(const double r) {
930  return dVdr_(r,Ri,sigmaPlated,epsilonPlated,densityPlated) -
931  dVdr_(r,Ro,sigmaPlated,epsilonPlated,densityPlated) +
932  dVdr_(r,Ro,sigma,epsilon,density);
933 }
934 
935 /**************************************************************************/
939 double PlatedLJCylinderPotential::valued2Vdr2(const double r) {
940  return 0.0;
941 }
942 
943 /**************************************************************************/
948 Array<dVec,1> PlatedLJCylinderPotential::initialConfig(const Container *boxPtr, MTRand &random,
949  const int numParticles) {
950 
951  /* The particle configuration */
952  Array<dVec,1> initialPos(numParticles);
953  initialPos = 0.0;
954 
955  /* We shift the radius inward to account for the excluded volume from the
956  * hard wall. This represents the largest prism that can be put
957  * inside a cylinder. */
958  dVec lside;
959  lside[0] = lside[1] = sqrt(2.0)*(Ri-sigmaPlated);
960  lside[2] = boxPtr->side[NDIM-1];
961 
962  /* Get the linear size per particle */
963  double initSide = pow((1.0*numParticles/product(lside)),-1.0/(1.0*NDIM));
964 
965  /* We determine the number of initial grid boxes there are in
966  * in each dimension and compute their size */
967  int totNumGridBoxes = 1;
968  iVec numNNGrid;
969  dVec sizeNNGrid;
970 
971  for (int i = 0; i < NDIM; i++) {
972  numNNGrid[i] = static_cast<int>(ceil((lside[i] / initSide) - EPS));
973 
974  /* Make sure we have at least one grid box */
975  if (numNNGrid[i] < 1)
976  numNNGrid[i] = 1;
977 
978  /* Compute the actual size of the grid */
979  sizeNNGrid[i] = lside[i] / (1.0 * numNNGrid[i]);
980 
981  /* Determine the total number of grid boxes */
982  totNumGridBoxes *= numNNGrid[i];
983  }
984 
985  /* Now, we place the particles at the middle of each box */
986  PIMC_ASSERT(totNumGridBoxes>=numParticles);
987  dVec pos;
988  for (int n = 0; n < totNumGridBoxes; n++) {
989 
990  iVec gridIndex;
991  for (int i = 0; i < NDIM; i++) {
992  int scale = 1;
993  for (int j = i+1; j < NDIM; j++)
994  scale *= numNNGrid[j];
995  gridIndex[i] = (n/scale) % numNNGrid[i];
996  }
997 
998  for (int i = 0; i < NDIM; i++)
999  pos[i] = (gridIndex[i]+0.5)*sizeNNGrid[i] - 0.5*lside[i];
1000 
1001  boxPtr->putInside(pos);
1002 
1003  if (n < numParticles)
1004  initialPos(n) = pos;
1005  else
1006  break;
1007  }
1008 
1009  return initialPos;
1010 }
1011 
1012 // ---------------------------------------------------------------------------
1013 // ---------------------------------------------------------------------------
1014 // LJ CYLINDER POTENTIAL CLASS -----------------------------------------------
1015 // ---------------------------------------------------------------------------
1016 // ---------------------------------------------------------------------------
1017 
1018 /**************************************************************************/
1028  PotentialBase(),
1030 {
1031 
1032  /* The radius of the tube */
1033  R = radius;
1034 
1035  /* The density of nitrogen in silicon nitride */
1036  density = 0.078; // atoms / angstrom^3
1037 // density = 0.008; // atoms / angstrom^3
1038 
1039  /* We define the values of epsilon and sigma for N and He */
1040 // double epsilonHe = 10.216; // Kelvin
1041 // double sigmaHe = 2.556; // angstroms
1042 // double sigmaN = 3.299; // angstroms
1043 // double epsilonN = 36.2; // Kelvin
1044 // epsilon = sqrt(epsilonHe*epsilonN);
1045 // sigma = 0.5*(sigmaHe + sigmaN);
1046 
1047  /* We operate under the assumption that the silicon can be neglected in the
1048  * silicon-nitride, and thus only consider the Nitrogen. We use a
1049  * Kiselov type model to extract the actual parameters. We assume that
1050  * silicate and silicon-nitride are roughly equivalent. */
1051  epsilon = 10.22; // Kelvin
1052  sigma = 2.628; // angstroms
1053 
1054 // epsilon = 32; // Kelvin
1055 // sigma = 3.08; // angstroms
1056 
1057  /* We choose a mesh consisting of 10^6 points, and create the lookup table */
1058  dR = (1.0E-6)*R;
1059 
1060  initLookupTable(dR,R);
1061 
1062  /* Find the minimun of the potential */
1063  minV = 1.0E5;
1064  for (int n = 0; n < tableLength; n++) {
1065  if (lookupV(n) < minV)
1066  minV = lookupV(n);
1067  }
1068 
1069  /* The extremal values for the lookup table */
1070  extV = valueV(0.0),valueV(R);
1071  extdVdr = valuedVdr(0.0),valuedVdr(R);
1072 }
1073 
1074 /**************************************************************************/
1078 }
1079 
1080 /**************************************************************************/
1086 double LJCylinderPotential::valueV(const double r) {
1087  double x = r / R;
1088 
1089  if (x >= 1.0)
1090  x = 1.0 - EPS;
1091 
1092  double x2 = x*x;
1093  double x4 = x2*x2;
1094  double x6 = x2*x4;
1095  double x8 = x4*x4;
1096  double f1 = 1.0 / (1.0 - x2);
1097  double sigoR3 = pow(sigma/R,3.0);
1098  double sigoR9 = sigoR3*sigoR3*sigoR3;
1099 
1100  double Kx2 = boost::math::ellint_1(x);
1101  double Ex2 = boost::math::ellint_2(x);
1102 
1103  double v9 = (1.0*pow(f1,9.0)/(240.0)) * (
1104  (1091.0 + 11156*x2 + 16434*x4 + 4052*x6 + 35*x8)*Ex2 -
1105  8.0*(1.0 - x2)*(1.0 + 7*x2)*(97.0 + 134*x2 + 25*x4)*Kx2);
1106  double v3 = 2.0*pow(f1,3.0) * ((7.0 + x2)*Ex2 - 4.0*(1.0-x2)*Kx2);
1107 
1108  return ((M_PI*epsilon*sigma*sigma*sigma*density/3.0)*(sigoR9*v9 - sigoR3*v3));
1109 }
1110 
1111 /**************************************************************************/
1116 double LJCylinderPotential::valuedVdr(const double r) {
1117 
1118  double x = r / R;
1119 
1120  if (x >= 1.0)
1121  x = 1.0 - EPS;
1122 
1123  /* dV/dr */
1124  if (x < EPS)
1125  return (1.28121E8/pow(R,11.0) - 102245.0/pow(R,5.0))*x;
1126  else {
1127  double x2 = x*x;
1128  double x4 = x2*x2;
1129  double x6 = x2*x4;
1130  double x8 = x4*x4;
1131  double f1 = 1.0 / (1.0 - x2);
1132  double sigoR3 = pow(sigma/R,3.0);
1133  double sigoR9 = sigoR3*sigoR3*sigoR3;
1134 
1135  double Kx2 = boost::math::ellint_1(x);
1136  double Ex2 = boost::math::ellint_2(x);
1137 
1138  double dv9dx =(3.0*pow(f1,10.0)/(80.0*x)) *
1139  ( (1.0 + x2)*(35.0 + 5108*x2 + 22482*x4 + 5108*x6 + 35*x8)*Ex2 -
1140  (1.0 - x2)*(35.0 + 3428*x2 + 15234*x4 + 12356*x6 +1715*x8)*Kx2 );
1141  double dv3dx = (6.0*pow(f1,4.0)/x) *
1142  ( (1.0 + 14*x2 + x4)*Ex2 - (1.0 + 6*x2 - 7*x4)*Kx2 );
1143  return ((M_PI*epsilon*sigma*sigma*sigma*density/(3.0*R))*(sigoR9*dv9dx - sigoR3*dv3dx));
1144  }
1145 }
1146 
1147 /**************************************************************************/
1152 double LJCylinderPotential::valued2Vdr2(const double r) {
1153 
1154  double x = r / R;
1155 
1156  if (x >= 1.0)
1157  x = 1.0 - EPS;
1158 
1159  /* d2V/dr2 */
1160  /*if (x < EPS){
1161  // related to hard core limit, this will likely need to be implemented.
1162  return (1.28121E8/pow(R,11.0) - 102245.0/pow(R,5.0))*x;
1163  }
1164  else {*/
1165  double x2 = x*x;
1166  double x4 = x2*x2;
1167  double x6 = x2*x4;
1168  double x8 = x4*x4;
1169  double x10 = x8*x2;
1170  double f1 = 1.0 / (1.0 - x2);
1171  double sigoR3 = pow(sigma/R,3.0);
1172  double sigoR9 = sigoR3*sigoR3*sigoR3;
1173 
1174  double Kx2 = boost::math::ellint_1(x);
1175  double Ex2 = boost::math::ellint_2(x);
1176 
1177  double d2v9dx2 = (2.0/240)*pow(f1,11)*((1925.0*x10 + 319451.0*x8 + 2079074.0*x6 + 2711942.0*x4
1178  + 764873.0*x2 + 20975.0)*Ex2
1179  - (729400.0*x10 + 430024.0*x8 + 344752.0*x6 + 767248.0*x4 + 386200.0*x2 + 12712.0)*Kx2); // checked --MTG
1180  double d2v3dx2 = 8.0*pow(f1,5)*((11.0 + 80.0*x2 + 5.0*x4)*Ex2 + 4.0*(5.0*x4 - 4.0*x2 - 1.0)*Kx2); // checked --MTG
1181 
1182  return ((M_PI*epsilon*sigma*sigma*sigma*density/(3.0*R*R))*(sigoR9*d2v9dx2 - sigoR3*d2v3dx2));
1183  //}
1184 }
1185 
1186 /**************************************************************************/
1191 Array<dVec,1> LJCylinderPotential::initialConfig(const Container *boxPtr, MTRand &random,
1192  const int numParticles) {
1193 
1194  /* The particle configuration */
1195  Array<dVec,1> initialPos(numParticles);
1196  initialPos = 0.0;
1197 
1198  /* We shift the radius inward to account for the excluded volume from the
1199  * hard wall. This represents the largest prism that can be put
1200  * inside a cylinder. */
1201  dVec lside;
1202  lside[0] = lside[1] = sqrt(2.0)*(R-sigma);
1203  lside[2] = boxPtr->side[NDIM-1];
1204 
1205  /* Get the linear size per particle */
1206  double initSide = pow((1.0*numParticles/product(lside)),-1.0/(1.0*NDIM));
1207 
1208  /* We determine the number of initial grid boxes there are in
1209  * in each dimension and compute their size */
1210  int totNumGridBoxes = 1;
1211  iVec numNNGrid;
1212  dVec sizeNNGrid;
1213 
1214  for (int i = 0; i < NDIM; i++) {
1215  numNNGrid[i] = static_cast<int>(ceil((lside[i] / initSide) - EPS));
1216 
1217  /* Make sure we have at least one grid box */
1218  if (numNNGrid[i] < 1)
1219  numNNGrid[i] = 1;
1220 
1221  /* Compute the actual size of the grid */
1222  sizeNNGrid[i] = lside[i] / (1.0 * numNNGrid[i]);
1223 
1224  /* Determine the total number of grid boxes */
1225  totNumGridBoxes *= numNNGrid[i];
1226  }
1227 
1228  /* Now, we place the particles at the middle of each box */
1229  PIMC_ASSERT(totNumGridBoxes>=numParticles);
1230  dVec pos;
1231  for (int n = 0; n < totNumGridBoxes; n++) {
1232 
1233  iVec gridIndex;
1234  for (int i = 0; i < NDIM; i++) {
1235  int scale = 1;
1236  for (int j = i+1; j < NDIM; j++)
1237  scale *= numNNGrid[j];
1238  gridIndex[i] = (n/scale) % numNNGrid[i];
1239  }
1240 
1241  for (int i = 0; i < NDIM; i++)
1242  pos[i] = (gridIndex[i]+0.5)*sizeNNGrid[i] - 0.5*lside[i];
1243 
1244  boxPtr->putInside(pos);
1245 
1246  if (n < numParticles)
1247  initialPos(n) = pos;
1248  else
1249  break;
1250  }
1251 
1252  return initialPos;
1253 }
1254 
1255 // ---------------------------------------------------------------------------
1256 // ---------------------------------------------------------------------------
1257 // LJ HOUR GLASS POTENTIAL CLASS ---------------------------------------------
1258 // ---------------------------------------------------------------------------
1259 // ---------------------------------------------------------------------------
1260 
1261 /**************************************************************************/
1274  const double deltaRadius, const double deltaWidth) : PotentialBase()
1275 {
1276  /* The radius of the tube */
1277  R = radius;
1278 
1279  /* The modification in radius: R(z=0) = R - dR */
1280  dR = deltaRadius;
1281 
1282  /* the length of the pore */
1283  L = constants()->L();
1284 
1285  /* The inverse length scale over which the radius changes */
1286  invd = 1.0/deltaWidth;
1287  R0 = 1.0/tanh(0.5*L*invd);
1288 
1289  /* The density of nitrogen in silicon nitride */
1290  density = 0.078; // atoms / angstrom^3
1291 
1292  /* These are the values we have historically used for amorphous silicon
1293  * nitride. */
1294  epsilon = 10.22; // Kelvin
1295  sigma = 2.628; // angstroms
1296 }
1297 
1298 /**************************************************************************/
1302 }
1303 
1304 /**************************************************************************/
1314 double LJHourGlassPotential::V(const dVec &r) {
1315 
1316  double x = 0.0;
1317  for (int i = 0; i < NDIM-1; i++)
1318  x += r[i]*r[i];
1319 
1320  double Rz = Rtanh(r[NDIM-1]);
1321 
1322  x = sqrt(x)/Rz;
1323 
1324  if (x >= 1.0)
1325  x = 1.0 - EPS;
1326 
1327  double x2 = x*x;
1328  double x4 = x2*x2;
1329  double x6 = x2*x4;
1330  double x8 = x4*x4;
1331  double f1 = 1.0 / (1.0 - x2);
1332  double sigoR3 = pow(sigma/Rz,3.0);
1333  double sigoR9 = sigoR3*sigoR3*sigoR3;
1334 
1335  double Kx2 = boost::math::ellint_1(x);
1336  double Ex2 = boost::math::ellint_2(x);
1337 
1338  double v9 = (1.0*pow(f1,9.0)/(240.0)) * (
1339  (1091.0 + 11156*x2 + 16434*x4 + 4052*x6 + 35*x8)*Ex2 -
1340  8.0*(1.0 - x2)*(1.0 + 7*x2)*(97.0 + 134*x2 + 25*x4)*Kx2);
1341  double v3 = 2.0*pow(f1,3.0) * ((7.0 + x2)*Ex2 - 4.0*(1.0-x2)*Kx2);
1342 
1343  return ((M_PI*epsilon*sigma*sigma*sigma*density/3.0)*(sigoR9*v9 - sigoR3*v3));
1344 }
1345 
1346 /**************************************************************************/
1355  MTRand &random, const int numParticles) {
1356 
1357  /* The particle configuration */
1358  Array<dVec,1> initialPos(numParticles);
1359  initialPos = 0.0;
1360 
1361  /* We shift the radius inward to account for the excluded volume from the
1362  * hard wall. This represents the largest prism that can be put
1363  * inside a cylinder. */
1364  dVec lside;
1365  lside[0] = lside[1] = sqrt(2.0)*(R-dR-sigma);
1366  lside[2] = boxPtr->side[NDIM-1];
1367 
1368  /* Make sure our radius isn't too small */
1369  PIMC_ASSERT(lside[0]>0.0);
1370 
1371  /* Get the linear size per particle */
1372  double initSide = pow((1.0*numParticles/product(lside)),-1.0/(1.0*NDIM));
1373 
1374  /* We determine the number of initial grid boxes there are in
1375  * in each dimension and compute their size */
1376  int totNumGridBoxes = 1;
1377  iVec numNNGrid;
1378  dVec sizeNNGrid;
1379 
1380  for (int i = 0; i < NDIM; i++) {
1381  numNNGrid[i] = static_cast<int>(ceil((lside[i] / initSide) - EPS));
1382 
1383  /* Make sure we have at least one grid box */
1384  if (numNNGrid[i] < 1)
1385  numNNGrid[i] = 1;
1386 
1387  /* Compute the actual size of the grid */
1388  sizeNNGrid[i] = lside[i] / (1.0 * numNNGrid[i]);
1389 
1390  /* Determine the total number of grid boxes */
1391  totNumGridBoxes *= numNNGrid[i];
1392  }
1393 
1394  /* Now, we place the particles at the middle of each box */
1395  PIMC_ASSERT(totNumGridBoxes>=numParticles);
1396  dVec pos;
1397  for (int n = 0; n < totNumGridBoxes; n++) {
1398 
1399  iVec gridIndex;
1400  for (int i = 0; i < NDIM; i++) {
1401  int scale = 1;
1402  for (int j = i+1; j < NDIM; j++)
1403  scale *= numNNGrid[j];
1404  gridIndex[i] = (n/scale) % numNNGrid[i];
1405  }
1406 
1407  for (int i = 0; i < NDIM; i++)
1408  pos[i] = (gridIndex[i]+0.5)*sizeNNGrid[i] - 0.5*lside[i];
1409 
1410  boxPtr->putInside(pos);
1411 
1412  if (n < numParticles)
1413  initialPos(n) = pos;
1414  else
1415  break;
1416  }
1417 
1418  return initialPos;
1419 }
1420 
1421 /**************************************************************************/
1429 Array<dVec,1> LJHourGlassPotential::initialConfig(const Container *boxPtr,
1430  MTRand &random, const int numParticles) {
1431 
1432  /* The particle configuration */
1433  Array<dVec,1> initialPos(numParticles);
1434  initialPos = 0.0;
1435 
1436  dVec pos;
1437  pos = 0.0;
1438 
1439  /* We randomly place the particles inside the cylinder taking acount of the
1440  * pinched radius.
1441  * @see http://mathworld.wolfram.com/DiskPointPicking.html
1442  */
1443  for (int n = 0; n < numParticles; n++) {
1444 
1445  /* Uniform position along the pore */
1446  pos[NDIM-1] = L*(-0.5 + random.rand());
1447 
1448  /* Uniform position in a disk of z-dependent radius*/
1449  double theta = 2.0*M_PI*random.rand();
1450  double r = (Rtanh(pos[NDIM-1]) - sigma)*sqrt(random.rand());
1451 
1452  pos[0] = r*cos(theta);
1453  pos[1] = r*sin(theta);
1454 
1455  boxPtr->putInside(pos);
1456 
1457  initialPos(n) = pos;
1458  }
1459 
1460  return initialPos;
1461 }
1462 
1463 // ---------------------------------------------------------------------------
1464 // ---------------------------------------------------------------------------
1465 // AZIZ POTENTIAL CLASS ------------------------------------------------------
1466 // ---------------------------------------------------------------------------
1467 // ---------------------------------------------------------------------------
1468 
1469 /**************************************************************************/
1476 {
1477  /* Define all variables for the Aziz potential */
1478  /* R.A. Aziz et al. J. Chem. Phys. 70, 4330 (1979) */
1479  rm = 2.9673; // A
1480  A = 0.5448504E6;
1481  epsilon = 10.8; // K
1482  alpha = 13.353384;
1483  D = 1.241314;
1484  C6 = 1.3732412;
1485  C8 = 0.4253785;
1486  C10 = 0.1781;
1487 
1488  /* The extremal values are all zero here */
1489  extV = 0.0;
1490  extdVdr = 0.0;
1491  extd2Vdr2 = 0.0;
1492 
1493  /* We take the maximum possible separation */
1494  double L = _boxPtr->maxSep;
1495 
1496  /* Create the potential lookup tables */
1497  // initLookupTable(0.00005*rm,L);
1498  initLookupTable((1.0E-6)*rm,L);
1499 
1500  /* Now we compute the tail correction */
1501  double rmoL = rm / L;
1502  double rm3 = rm*rm*rm;
1503  double t1 = A*exp(-alpha*L/(2.0*rm))*rm*(8.0*rm*rm + 4.0*L*rm * alpha + L*L*alpha*alpha)
1504  / (4.0*alpha*alpha*alpha);
1505  double t2 = 8.0*C6*pow(rmoL,3.0)/3.0;
1506  double t3 = 32.0*C8*pow(rmoL,5.0)/5.0;
1507  double t4 = 128.0*C10*pow(rmoL,7.0)/7.0;
1508 
1509  tailV = 2.0*M_PI*epsilon*(t1 - rm3*(t2+t3+t4));
1510 }
1511 
1512 /**************************************************************************/
1516 }
1517 
1518 /**************************************************************************/
1524 double AzizPotential::valueV(const double r) {
1525  double x = r / rm;
1526 
1527  double Urep = A * exp(-alpha*x);
1528 
1529  /* No self interactions */
1530  if (x < EPS)
1531  return 0.0;
1532  /* Hard core limit */
1533  else if (x < 0.01)
1534  return (epsilon * Urep);
1535  else {
1536  double ix2 = 1.0 / (x * x);
1537  double ix6 = ix2 * ix2 * ix2;
1538  double ix8 = ix6 * ix2;
1539  double ix10 = ix8 * ix2;
1540  double Uatt = -( C6*ix6 + C8*ix8 + C10*ix10 ) * F(x);
1541  return ( epsilon * (Urep + Uatt) );
1542  }
1543 }
1544 
1545 /**************************************************************************/
1550 double AzizPotential::valuedVdr(const double r) {
1551  double x = r / rm;
1552 
1553  double T1 = -A * alpha * exp(-alpha*x);
1554 
1555  /* dV/dR */
1556  /* No self interactions */
1557  if (x < EPS)
1558  return 0.0;
1559  /* Hard core limit */
1560  else if (x < 0.01)
1561  return ( ( epsilon / rm ) * T1 );
1562  else {
1563  /* The various inverse powers of x */
1564  double ix = 1.0 / x;
1565  double ix2 = ix*ix;
1566  double ix6 = ix2 * ix2 * ix2;
1567  double ix7 = ix6 * ix;
1568  double ix8 = ix6 * ix2;
1569  double ix9 = ix8 * ix;
1570  double ix10 = ix8 * ix2;
1571  double ix11 = ix10 * ix;
1572  double T2 = ( 6.0*C6*ix7 + 8.0*C8*ix9 + 10.0*C10*ix11 ) * F(x);
1573  double T3 = -( C6*ix6 + C8*ix8 + C10*ix10 ) * dF(x);
1574  return ( ( epsilon / rm ) * (T1 + T2 + T3) );
1575  }
1576 }
1577 
1578 /**************************************************************************/
1583 double AzizPotential::valued2Vdr2(const double r) {
1584  double x = r / rm;
1585 
1586  double T1 = A * alpha * alpha * exp(-alpha*x);
1587 
1588  /* d^2V/dR^2 */
1589  /* No self interactions */
1590  if (x < EPS)
1591  return 0.0;
1592  /* Hard core limit */
1593  else if (x < 0.01)
1594  return ( ( epsilon / rm ) * T1 );
1595  else {
1596  /* The various inverse powers of x */
1597  double ix = 1.0 / x;
1598  double ix2 = ix*ix;
1599  double ix6 = ix2 * ix2 * ix2;
1600  double ix7 = ix6 * ix;
1601  double ix8 = ix6 * ix2;
1602  double ix9 = ix8 * ix;
1603  double ix10 = ix8 * ix2;
1604  double ix11 = ix10 * ix;
1605  double ix12 = ix11 * ix;
1606  double T2 = - ( 42.0*C6*ix8 + 72.0*C8*ix10 + 110.0*C10*ix12 ) * F(x);
1607  double T3 = 2.0*( 6.0*C6*ix7 + 8.0*C8*ix9 + 10.0*C10*ix11 ) * dF(x);
1608  double T4 = - ( C6*ix6 + C8*ix8 + C10*ix10 ) * d2F(x);
1609  return ( ( epsilon / (rm*rm) ) * (T1 + T2 + T3 + T4) );
1610  }
1611 }
1612 
1613 // ---------------------------------------------------------------------------
1614 // ---------------------------------------------------------------------------
1615 // SZALEWICZ POTENTIAL CLASS ------------------------------------------------------
1616 // ---------------------------------------------------------------------------
1617 // ---------------------------------------------------------------------------
1618 
1619 /**************************************************************************/
1626 {
1627  /* Define all variables for the Szalewicz potential */
1628  // FIXME Fix comments
1629  rm = 2.9262186279335958; // Angstrom (scipy.optimize.minimize())
1630  /* The extremal values are all zero here */
1631  extV = 0.0;
1632  extdVdr = 0.0;
1633  extd2Vdr2 = 0.0;
1634 
1635  /* We take the maximum possible separation */
1636  double L = _boxPtr->maxSep;
1637 
1638  // FIXME Get rid of this stuff or put it somwhere else
1639  /*
1640  double zz = L;
1641  double dr = 1.0e-5;
1642  double zmax = 1000.0;
1643  double tot = int((zmax-L)/dr);
1644  std::vector<double> outv2;
1645  for (int iii = 0; iii < tot; iii++) {
1646  double vv = valueV(zz);
1647  outv2.push_back(vv);
1648  zz += dr;
1649  }
1650 
1651 
1652  std::cout << L << std::endl;
1653  std::cout << outv2.size() << std::endl;
1654  std::ofstream output_file2("./SzalewiczTail2.txt");
1655  //std::ostream_iterator<double> output_iterator2(output_file2, "\n");
1656  //std::copy(outv2.begin(), outv2.end(), output_iterator2);
1657  write_container(outv2,output_file2);
1658  output_file2.close();
1659  exit(0);
1660  */
1661 
1662  /* Create the potential lookup tables */
1663  // initLookupTable(0.00005*rm,L);
1664  initLookupTable((1.0E-6)*rm,L);
1665 
1666  // FIXME fix the tail correction
1667  /* Now we compute the tail correction */
1668  /*
1669  L *= 1.8897261254578281; // Convert L to Bohr
1670  double t1 = 0.0;
1671 
1672  t1 -= exp(-a * L)*((pow(a,2)*Pn[0])+(a*(1+(a*L))*Pn[1])+((2+((a*L)*(2+(a*L))))*Pn[2]))/(pow(a,3));
1673  t1 += exp(-b*L)*L*((2*Qn[0])+(L*Qn[1]))/2;
1674 
1675  double t2 = 0.0;
1676  double eL = eta*L;
1677  for (int n = 3; n < 17; n++){
1678  t2 += exp(- eL) * pow(L,-n) * ((pow(eL,n)*(eL+1)) + (exp(eL)*eL*fn(eL,n)*factorials[n])) * Cn[n]/((n-1)*eta*factorials[n]);
1679  }
1680 
1681  tailV = -t1 - t2;
1682  tailV *= 315774.65;
1683  cout << tailV << endl;
1684  cout << t1*315774.65 << endl;
1685  cout << t2*315774.65 << endl;
1686 
1687  exit(0);
1688  */
1689 }
1690 
1691 /**************************************************************************/
1695 }
1696 
1697 /**************************************************************************/
1702 double SzalewiczPotential::valueV(const double _r) {
1703  double r = _r * 1.8897261254578281; // convert from angtrom to bohr
1704  double x = _r / rm;
1705 
1706  /* No self interactions */
1707  if (x < EPS)
1708  return 0.0;
1709  else {
1710  double t1 = exp(-a*r);
1711  double t1m = 0.00;
1712  for (int i = 0; i < 3; i++){
1713  t1m += Pn[i]*pow(r,i);
1714  }
1715  t1 *= t1m;
1716 
1717  double t2 = exp(-b*r);
1718  double t2m = 0.00;
1719  for (int i = 0; i < 2; i++){
1720  t2m += Qn[i]*pow(r,i);
1721  }
1722  t2 *= t2m;
1723  double t3 = 0.0;
1724 
1725  /* Hard core limit not reached */
1726  if (x > 0.10){
1727  for (int i = 0; i < 17; i++){
1728  t3 += fn(eta*r,i)*Cn[i]/pow(r,i);
1729  }
1730  }
1731  else {
1732  for (int i = 0; i < 17; i++){
1733  t3 += fn2(eta*r,i)*Cn[i]/pow(r,i);
1734  }
1735  }
1736  return (t1 + t2 - t3)*315774.65;
1737  }
1738 }
1739 
1740 /**************************************************************************/
1745 double SzalewiczPotential::valuedVdr(const double _r) {
1746  double r = _r * 1.8897261254578281; // convert from angtrom to bohr
1747  double x = _r / rm;
1748 
1749  /* No self interactions */
1750  if (x < EPS)
1751  return 0.0;
1752  else {
1753  double t1 = exp(-a*r);
1754  double t1m = 0.00;
1755  for (int i = 0; i < 3; i++){
1756  t1m += Pn[i]*pow(r,i);
1757  }
1758  t1 *= ((Pn[1]+(2*r*Pn[2]))-(a*t1m));
1759 
1760  double t2 = exp(-b*r);
1761  double t2m = 0.00;
1762  for (int i = 0; i < 2; i++){
1763  t2m += Qn[i]*pow(r,i);
1764  }
1765  t2 *= Qn[1]-b*t2m;
1766 
1767  double t3 = 0.0;
1768  double t4 = 0.0;
1769 
1770  /* Hard core limit not reached */
1771  if (x > 0.10){
1772 
1773  for (int i = 0; i < 16; i++){
1774  t3 += (i+1)*fn(eta*r,i+1)*Cn[i+1]/pow(r,(i+2));
1775  }
1776 
1777  for (int i = 0; i < 17; i++){
1778  t4 += eta*dfn(eta*r,i)*Cn[i]/pow(r,i);
1779  }
1780 
1781  }
1782  else {
1783 
1784  for (int i = 0; i < 16; i++){
1785  t3 += (i+1)*fn2(eta*r,i+1)*Cn[i+1]/pow(r,(i+2));
1786  }
1787 
1788  for (int i = 0; i < 17; i++){
1789  t4 += eta*dfn(eta*r,i)*Cn[i]/pow(r,i);
1790  }
1791  }
1792  return (t1 + t2 + t3 - t4)*315774.65;
1793  }
1794 }
1795 
1796 /**************************************************************************/
1801 double SzalewiczPotential::valued2Vdr2(const double _r) {
1802  double r = _r * 1.8897261254578281; // convert from angtrom to bohr
1803  double x = _r / rm;
1804 
1805  /* No self interactions */
1806  if (x < EPS)
1807  return 0.0;
1808  else {
1809  double t1 = exp(-a*r);
1810  double t1m = 0.00;
1811  for (int i = 0; i < 3; i++){
1812  t1m += Pn[i]*pow(r,i);
1813  }
1814  t1 *= (2*Pn[2]-(2*a*(Pn[1]+(2*r*Pn[2])))+(a*a*t1m));
1815 
1816  double t2 = exp(-b*r);
1817  double t2m = 0.00;
1818  for (int i = 0; i < 2; i++){
1819  t2m += Qn[i]*pow(r,i);
1820  }
1821  t2 *= (b*b*t2m) - (2*b*Qn[1]);
1822 
1823  double t3 = 0.0;
1824  double t4 = 0.0;
1825  double t5 = 0.0;
1826 
1827  for (int i = 0; i < 15; i++){
1828  t3 += (i+1)*(i+2)*fn2(eta*r,i+1)*Cn[i+1]/pow(r,(i+3));
1829  }
1830 
1831  for (int i = 0; i < 16; i++){
1832  t4 += 2*(i+1)*eta*dfn(eta*r,i+1)*Cn[i+1]/pow(r,(i+2));
1833  }
1834 
1835  for (int i = 0; i < 17; i++){
1836  t5 += eta*eta*d2fn(eta*r,i)*Cn[i]/pow(r,i);
1837  }
1838 
1839  return (t1 + t2 - t3 + t4 - t5)*315774.65;
1840  }
1841 }
1842 
1843 
1844 
1845 // ---------------------------------------------------------------------------
1846 // ---------------------------------------------------------------------------
1847 // EXCLUDED VOLUME CLASS -----------------------------------------------
1848 // ---------------------------------------------------------------------------
1849 // ---------------------------------------------------------------------------
1850 /**************************************************************************/
1853 Gasparini_1_Potential::Gasparini_1_Potential(double _az, double _ay, const Container *_boxPtr) :
1854  PotentialBase(),
1855  excZ(0.5*_boxPtr->side[2]-_az),
1856  excY(0.5*_boxPtr->side[1]-_ay),
1857  V0(1.0E6)
1858 {
1859  // Empty Constructor
1860 }
1861 
1862 /**************************************************************************/
1866  // Empty Destructor
1867 }
1868 
1869 /**************************************************************************/
1875  MTRand &random, const int numParticles) {
1876 
1877  /* label the lengths of the sides of the simulation cell */
1878  dVec lside;
1879  lside[0] = boxPtr->side[0];
1880  lside[1] = boxPtr->side[1];
1881  lside[2] = boxPtr->side[NDIM-1];
1882 
1883  /* calculate actual volume */
1884  double volTot = product(lside);
1885 
1886  /* calculate density */
1887  double density = (1.0*numParticles/volTot);
1888 
1889  /* calculate excluded volume */
1890  double volExc = lside[0]*(2.0*excY)*(2.0*excZ);
1891 
1892  /* calculate actual number of particles */
1893  int correctNum = int(numParticles-density*volExc);
1894 
1895  /* The particle configuration */
1896  Array<dVec,1> initialPos(correctNum);
1897 
1898  /* get linear size per particle */
1899  double initSide = pow((1.0*correctNum/(volTot-volExc)),-1.0/(1.0*NDIM));
1900 
1901  /* For accessible volume, determine the number of
1902  * initial grid boxes there are in each dimension and compute
1903  * their size. */
1904  int totNumGridBoxes1 = 1;
1905  int totNumGridBoxes2 = 1;
1906  iVec numNNGrid1;
1907  iVec numNNGrid2;
1908  dVec sizeNNGrid1;
1909  dVec sizeNNGrid2;
1910 
1911  /* divide space into two regions, insert particles appropriately */
1912  double V1 = (lside[1]-2.0*excY)*(2.0*excZ)*lside[0];
1913  double V2 = (lside[2]-2.0*excZ)*lside[1]*lside[0];
1914 
1915  double fracV1 = V1/(V1+V2);
1916 
1917  int numIn1 = int(correctNum*fracV1);
1918  int numIn2 = (correctNum-numIn1);
1919 
1920  /* grid space in volume 1 */
1921  /* x */
1922  numNNGrid1[0] = static_cast<int>(ceil((1.0*lside[0]/initSide)-EPS));
1923  if (numNNGrid1[0] < 1)
1924  numNNGrid1[0] = 1;
1925  sizeNNGrid1[0] = 1.0*lside[0]/(1.0*numNNGrid1[0]);
1926  totNumGridBoxes1 *= numNNGrid1[0];
1927  /* y */
1928  numNNGrid1[1] = static_cast<int>(ceil(((lside[1]-2.0*excY)/initSide)-EPS));
1929  if (numNNGrid1[1] < 1)
1930  numNNGrid1[1] = 1;
1931  sizeNNGrid1[1] = (lside[1]-2.0*excY)/(1.0*numNNGrid1[1]);
1932  totNumGridBoxes1 *= numNNGrid1[1];
1933  /* z */
1934  numNNGrid1[2] = static_cast<int>(ceil(((2.0*excZ)/initSide)-EPS));
1935  if (numNNGrid1[2] < 1)
1936  numNNGrid1[2] = 1;
1937  sizeNNGrid1[2] = (2.0*excZ)/(1.0*numNNGrid1[2]);
1938  totNumGridBoxes1 *= numNNGrid1[2];
1939 
1940  /* grid space in volume 2 */
1941  /* x */
1942  numNNGrid2[0] = static_cast<int>(ceil((1.0*lside[0]/initSide)-EPS));
1943  if (numNNGrid2[0] < 1)
1944  numNNGrid2[0] = 1;
1945  sizeNNGrid2[0] = 1.0*lside[0]/(1.0*numNNGrid2[0]);
1946  totNumGridBoxes2 *= numNNGrid2[0];
1947  /* y */
1948  numNNGrid2[1] = static_cast<int>(ceil((1.0*lside[1]/initSide)-EPS));
1949  if (numNNGrid2[1] < 1)
1950  numNNGrid2[1] = 1;
1951  sizeNNGrid2[1] = 1.0*lside[1]/(1.0*numNNGrid2[1]);
1952  totNumGridBoxes2 *= numNNGrid2[1];
1953  /* z */
1954  numNNGrid2[2] = static_cast<int>(ceil(((lside[2]-2.0*excZ)/initSide)-EPS));
1955  if (numNNGrid2[2] < 1)
1956  numNNGrid2[2] = 1;
1957  sizeNNGrid2[2] = (lside[2]-2.0*excZ)/(1.0*numNNGrid2[2]);
1958  totNumGridBoxes2 *= numNNGrid2[2];
1959 
1960  /* Place particles in the middle of the boxes -- volume 1 */
1961  PIMC_ASSERT(totNumGridBoxes1>=numIn1);
1962  dVec pos1;
1963 
1964  for (int n = 0; n < totNumGridBoxes1; n++) {
1965  iVec gridIndex1;
1966  /* update grid index */
1967  for (int i = 0; i < NDIM; i++) {
1968  int scale = 1;
1969  for (int j = i+1; j < NDIM; j++)
1970  scale *= numNNGrid1[j];
1971  gridIndex1[i] = (n/scale) % numNNGrid1[i];
1972  }
1973  /* place particle in position vector, skipping over excluded volume */
1974  pos1[0] = (gridIndex1[0]+0.5)*sizeNNGrid1[0] + +0.5*lside[0] + 2.0*EPS;
1975  pos1[1] = (gridIndex1[1]+0.5)*sizeNNGrid1[1] + excY + 2.0*EPS;
1976  pos1[2] = (gridIndex1[2]+0.5)*sizeNNGrid1[2] - excZ + 2.0*EPS;
1977 
1978  if ((pos1[1]<-excY || pos1[1]>excY) || (pos1[2]<-excZ || pos1[2]>excZ))
1979  boxPtr->putInside(pos1);
1980 
1981  if (n < numIn1){
1982  initialPos(n) = pos1;
1983  }
1984  else
1985  break;
1986  }
1987 
1988  /* Place particles in the middle of the boxes -- volume 2 */
1989  PIMC_ASSERT(totNumGridBoxes2>=numIn2);
1990  dVec pos2;
1991 
1992  for (int n = 0; n < totNumGridBoxes2; n++) {
1993  iVec gridIndex2;
1994  /* update grid index */
1995  for (int i = 0; i < NDIM; i++) {
1996  int scale = 1;
1997  for (int j = i+1; j < NDIM; j++)
1998  scale *= numNNGrid2[j];
1999  gridIndex2[i] = (n/scale) % numNNGrid2[i];
2000  }
2001  /* place particles in position vectors */
2002  pos2[0] = (gridIndex2[0]+0.5)*sizeNNGrid2[0] + 0.5*lside[0] + 2.0*EPS;
2003  pos2[1] = (gridIndex2[1]+0.5)*sizeNNGrid2[1] + 0.5*lside[1] + 2.0*EPS;
2004  pos2[2] = (gridIndex2[2]+0.5)*sizeNNGrid2[2] + excZ + 2.0*EPS;
2005 
2006  if ((pos2[1]<-excY || pos2[1]>excY) || (pos2[2]<-excZ || pos2[2]>excZ))
2007  boxPtr->putInside(pos2);
2008 
2009  if (n < numIn2){
2010  initialPos(n+numIn1) = pos2;
2011  }
2012  else
2013  break;
2014  }
2015  /* do we want to output the initial config to disk? */
2016  bool outToDisk = 1;
2017  ofstream OF;
2018  if (outToDisk){
2019  OF.open("./OUTPUT/initialConfig.dat");
2020  OF<<"# Cartesian Coordinates of initial Positions (X-Y-Z)"<<endl;
2021  OF<<"# "<<lside[0]<<"\t"<<lside[1]<<"\t"<<lside[2]<<endl;
2022  OF<<"# "<< excY <<"\t"<< excZ <<endl;
2023  for (int i=0; i< int(initialPos.size()); i++)
2024  OF<<initialPos(i)(0)<< "\t"<<initialPos(i)(1)<<"\t"<<initialPos(i)(2)<<endl;
2025  OF.close();
2026  }
2027 
2028  return initialPos;
2029 }
2030 /*************************************************************************/
2034 
2035  Array<double, 1> excLens(2);
2036  excLens(0) = excY; // was ay
2037  excLens(1) = excZ; // was az
2038 
2039  return (excLens);
2040 }
2041 
2042 
2043 // ---------------------------------------------------------------------------
2044 // ---------------------------------------------------------------------------
2045 // HARD SPHERE POTENTIAL CLASS -----------------------------------------------
2046 // ---------------------------------------------------------------------------
2047 // ---------------------------------------------------------------------------
2048 
2049 /**************************************************************************/
2054  PotentialBase(),
2055  a(_a) {
2056 
2057 }
2058 
2059 /**************************************************************************/
2063 // empty deconstructor
2064 }
2065 
2066 /**************************************************************************/
2077 double HardSpherePotential::V(const dVec &sep1, const dVec &sep2)
2078 {
2079 
2080  double r1 = sqrt(dot(sep1,sep1));
2081  double r2 = sqrt(dot(sep2,sep2));
2082 
2083  if ((r1 <= a ) || (r2 <= a))
2084  return LBIG;
2085 
2086  double cosTheta = dot(sep1,sep2)/(r1*r2);
2087 
2088  double t1 = -(r1*r2 + a*a - a*(r1 + r2)) * (1.0 + cosTheta);
2089  t1 /= (4.0*constants()->lambda()*constants()->tau());
2090 
2091  double t2 = (a*(r1+r2) - a*a)/(r1*r2);
2092  double t3 = 1.0 - t2*exp(t1);
2093 
2094  return -log(t3);
2095 }
2096 
2097 /**************************************************************************/
2108 double HardSpherePotential::dVdlambda(const dVec &sep1, const dVec &sep2)
2109 {
2110 
2111  double r1 = sqrt(dot(sep1,sep1));
2112  double r2 = sqrt(dot(sep2,sep2));
2113 
2114  double cosTheta = dot(sep1,sep2)/(r1*r2);
2115 
2116  double t1 = -(r1*r2 + a*a - a*(r1 + r2)) * (1.0 + cosTheta);
2117  t1 /= (4.0*constants()->lambda()*constants()->tau());
2118 
2119  double t2 = (a*(r1+r2) - a*a)/(r1*r2);
2120  double t3 = 1.0 - t2*exp(t1);
2121 
2122  double t4 = (1.0/t3)*(t2*exp(t1))*(t1/constants()->lambda());
2123 
2124  return -t4;
2125 }
2126 
2127 /**************************************************************************/
2140 double HardSpherePotential::dVdtau(const dVec &sep1, const dVec &sep2)
2141 {
2142 
2143  double r1 = sqrt(dot(sep1,sep1));
2144  double r2 = sqrt(dot(sep2,sep2));
2145 
2146  double cosTheta = dot(sep1,sep2)/(r1*r2);
2147 
2148  double t1 = -(r1*r2 + a*a - a*(r1 + r2)) * (1.0 + cosTheta);
2149  t1 /= (4.0*constants()->lambda()*constants()->tau());
2150 
2151  double t2 = (a*(r1+r2) - a*a)/(r1*r2);
2152  double t3 = 1.0 - t2*exp(t1);
2153 
2154  double t4 = (1.0/t3)*(t2*exp(t1))*(t1/constants()->tau());
2155 
2156  return -t4;
2157 }
2158 
2159 
2160 // ---------------------------------------------------------------------------
2161 // ---------------------------------------------------------------------------
2162 // 1D Delta Potential Class -----------------------------------------------
2163 // ---------------------------------------------------------------------------
2164 // ---------------------------------------------------------------------------
2165 
2166 /**************************************************************************/
2171  PotentialBase(),
2172  g(_g)
2173 {
2174 
2175  erfCO = 7.0;
2176  l0 = 2.0*sqrt(constants()->lambda()*constants()->tau());
2177  li = 4.0*constants()->lambda()/g;
2178  xi = l0/li;
2179 
2180  xiSqOver2 = (0.5)*xi*xi;
2181  xiSqrtPIOver2 = sqrt(M_PI/2.0)*xi;
2182 }
2183 
2184 /**************************************************************************/
2188  // empty deconstructor
2189 }
2190 
2191 
2192 /**************************************************************************/
2198 double Delta1DPotential::Wint(double yt, double dxt) {
2199 
2200  double W,erfVal,expVal;
2201 
2202  if(xi + yt < erfCO)
2203  {
2204  erfVal = erfc( (xi+yt)/sqrt(2.0) );
2205  expVal = exp( xiSqOver2 + (0.5)*dxt*dxt + xi*yt );
2206  W = 1.0 - xiSqrtPIOver2*expVal*erfVal;
2207  }
2208  else
2209  {
2210  expVal = exp((-0.5)*(yt*yt-dxt*dxt));
2211  W = 1.0 - (xi/(xi+yt))*expVal;
2212  }
2213 
2214  return W;
2215 }
2216 
2217 /**************************************************************************/
2223 double Delta1DPotential::dWdxi(double yt, double dxt) {
2224 
2225  double dW,erfVal,expVal,expVal2;
2226 
2227  expVal = exp( (0.5)*(dxt*dxt-yt*yt) );
2228 
2229  if(xi + yt < erfCO)
2230  {
2231  erfVal = erfc( (xi+yt)/sqrt(2.0) );
2232  expVal2 = exp( (0.5)*(xi+yt)*(xi+yt) );
2233  dW = (-1.0)*xi*expVal*( sqrt(0.5*M_PI)*(xi +yt + (1.0/xi) )*expVal2*erfVal - 1.0 );
2234  }
2235  else
2236  dW = (-1.0)*xi*expVal*( (xi + yt + (1.0/xi) )/( xi + yt ) - 1.0 );
2237 
2238  return dW;
2239 }
2240 
2241 /**************************************************************************/
2247 double Delta1DPotential::dWdyt(double yt, double dxt) {
2248 
2249  double dW,erfVal,expVal,expVal2;
2250 
2251  expVal = exp( (0.5)*(dxt*dxt-yt*yt) );
2252 
2253  if(xi + yt < erfCO)
2254  {
2255  erfVal = erfc( (xi+yt)/sqrt(2.0) );
2256  expVal2 = exp( (0.5)*(xi+yt)*(xi+yt) );
2257  dW = xi*expVal*( 1 - sqrt(0.5*M_PI)*xi*expVal2*erfVal );
2258  }
2259  else
2260  dW = xi*expVal*( 1 - xi/(xi+yt) );
2261 
2262  return dW;
2263 }
2264 
2265 /**************************************************************************/
2271 double Delta1DPotential::dWddxt(double yt, double dxt) {
2272 
2273  double dW,erfVal,expVal;
2274 
2275  if (xi + yt < erfCO)
2276  {
2277  erfVal = erfc( (xi+yt)/sqrt(2.0) );
2278  expVal = exp( (0.5)*( dxt*dxt + xi*( xi + 2.0*yt)) );
2279  dW = (-1.0)*sqrt(0.5*M_PI)*xi*dxt*expVal*erfVal;
2280  }
2281  else
2282  {
2283  expVal = exp( (0.5)*(dxt*dxt-yt*yt) );
2284  dW = (-1.0)*dxt*expVal*( xi/( xi + yt ) );
2285  }
2286 
2287  return dW;
2288 }
2289 
2290 /**************************************************************************/
2299 double Delta1DPotential::V(const dVec &sep1, const dVec &sep2)
2300 {
2301 
2302  double dxt = deltaSeparation(sep1[0], sep2[0])/l0;
2303  double yt = (abs(sep1[0])+abs(sep2[0]))/l0;
2304 
2305  double W = Wint(yt,dxt);
2306 
2307  return (-1.0)*log(W);
2308 }
2309 
2310 /**************************************************************************/
2319 double Delta1DPotential::dVdlambda(const dVec &sep1, const dVec &sep2)
2320 {
2321 
2322  double dxt = deltaSeparation(sep1[0], sep2[0])/l0;
2323  double yt = (abs(sep1[0])+abs(sep2[0]))/l0;
2324 
2325  double W = Wint(yt,dxt);
2326  double dWdy = dWdyt(yt,dxt);
2327  double dWddx = dWddxt(yt,dxt);
2328  double dWdx = dWdxi(yt,dxt);
2329 
2330  double dWdl = ((-1.0)/(2.0*constants()->lambda()))*(yt*dWdy + dxt*dWddx+ xi*dWdx );
2331 
2332  return ((-1.0)/W)*dWdl;
2333 }
2334 
2335 /**************************************************************************/
2344 double Delta1DPotential::dVdtau(const dVec &sep1, const dVec &sep2)
2345 {
2346 
2347  double dxt = deltaSeparation(sep1[0], sep2[0])/l0;
2348  double yt = (abs(sep1[0])+abs(sep2[0]))/l0;
2349 
2350  double W = Wint(yt,dxt);
2351  double dWdt = ((1.0)/(2.0*constants()->tau()))
2352  * ( (-1.0)*yt*dWdyt(yt,dxt) - dxt*dWddxt(yt,dxt) + xi*dWdxi(yt,dxt) );
2353 
2354  return ((-1.0)/W)*dWdt;
2355 }
2356 
2357 // ---------------------------------------------------------------------------
2358 // ---------------------------------------------------------------------------
2359 // HARD ROD POTENTIAL CLASS --------------------------------------------------
2360 // ---------------------------------------------------------------------------
2361 // ---------------------------------------------------------------------------
2362 
2363 /**************************************************************************/
2368  PotentialBase(),
2369  a(_a) {
2370 
2371 }
2372 
2373 /**************************************************************************/
2377 // empty deconstructor
2378 }
2379 
2380 /**************************************************************************/
2391 double HardRodPotential::V(const dVec &sep1, const dVec &sep2)
2392 {
2393  double r1 = sqrt(dot(sep1,sep1));
2394  double r2 = sqrt(dot(sep2,sep2));
2395 
2396  /* We need to enforce the distinguishable particle constraint at short
2397  * imaginary times */
2398  /* if ( (sep1[0]*sep2[0] < 0.0) || (r1 <= a ) || (r2 <= a) ) */
2399  /* return LBIG; */
2400 
2401  if ( (r1 <= a ) || (r2 <= a) )
2402  return LBIG;
2403 
2404  double d1 = deltaSeparation(r1,a);
2405  double d2 = deltaSeparation(r2,a);
2406 
2407  double t1 = -d1*d2/(2.0*constants()->lambda()*constants()->tau());
2408 
2409  /* communicate()->file("debug")->stream() << sep1[0] << "\t" << sep2[0] << "\t" */
2410  /* << -log(1.0-exp(t1)) << endl; */
2411 
2412  return (-log(1.0 - exp(t1)));
2413 }
2414 
2415 /**************************************************************************/
2426 double HardRodPotential::dVdlambda(const dVec &sep1, const dVec &sep2)
2427 {
2428 
2429  double r1 = sqrt(dot(sep1,sep1));
2430  double r2 = sqrt(dot(sep2,sep2));
2431  double d1 = deltaSeparation(r1,a);
2432  double d2 = deltaSeparation(r2,a);
2433 
2434  double t1 = d1*d2;
2435  double t2 = t1/(2.0*constants()->lambda()*constants()->lambda()*constants()->tau());
2436 
2437  return ((0.5*t1/(exp(t2)-1.0))/(constants()->lambda()*constants()->lambda()*constants()->tau()));
2438 }
2439 
2440 /**************************************************************************/
2451 double HardRodPotential::dVdtau(const dVec &sep1, const dVec &sep2)
2452 {
2453 
2454  double r1 = sqrt(dot(sep1,sep1));
2455  double r2 = sqrt(dot(sep2,sep2));
2456  double d1 = deltaSeparation(r1,a);
2457  double d2 = deltaSeparation(r2,a);
2458 
2459  double t1 = d1*d2;
2460  double t2 = t1/(2.0*constants()->lambda()*constants()->tau());
2461 
2462  return ((0.5*t1/(exp(t2)-1.0))/(constants()->lambda()*constants()->tau()*constants()->tau()));
2463 }
2464 
2465 // ---------------------------------------------------------------------------
2466 // ---------------------------------------------------------------------------
2467 // GraphenePotential Class----------------------------------------------------
2468 // ---------------------------------------------------------------------------
2469 // ---------------------------------------------------------------------------
2470 
2471 /**************************************************************************/
2474 #include <boost/math/special_functions/bessel.hpp>
2475 GraphenePotential::GraphenePotential (double _strain, double _poisson, double _a0,
2476  double _sigma, double _epsilon) : PotentialBase() {
2477  double strain = _strain;
2478  double poisson = _poisson;
2479  double a0 = _a0;
2480 
2481  sigma = _sigma;
2482  epsilon = _epsilon;
2483 
2484  Lz = constants()->L();
2485  Lzo2 = Lz/2.0;
2486 
2487  /* Lattice vectors */
2488  /* @see: https://wiki.cmt.w3.uvm.edu/index.php?title=Bose_Hubbard_model_treatment_for_Helium_absorbed_on_graphene#strain */
2489  a1x = sqrt(3.0)*a0*0.5*(1.0-strain*poisson);
2490  a1y = 3.0*a0*0.5*(1+strain);
2491  a2x = -a1x;
2492  a2y = a1y;
2493 
2494  /* reciprocal lattice vectors */
2495  g1x = 2*M_PI/(sqrt(3.0)*a0*(1.0-strain*poisson));
2496  g1y = 2*M_PI/(3.0*a0*(1+strain));
2497  g2x = -g1x;
2498  g2y = g1y;
2499 
2500  /* basis vectors */
2501  b1x = sqrt(3.0)*a0*0.5*(1.0-strain*poisson);
2502  b1y = 0.5*a0*(1.0 + strain);
2503  b2x = 0.0;
2504  b2y = a0*(1+strain);
2505 
2506  /* area of unit cell */
2507  A = fabs((a1x*a2y) - (a1y*a2x));
2508 
2509 }
2510 
2511 
2512 /**************************************************************************/
2516 
2517 }
2518 
2519 /**************************************************************************/
2526 double GraphenePotential::V(const dVec &r) {
2527 
2528  double z = r[2] + Lzo2;
2529 
2530  /* We take care of the particle being in a forbidden region */
2531  if (z < 1.5)
2532  return 40000.;
2533  if (z > Lz-0.05)
2534  return 40000.;
2535 
2536  double x = r[0];
2537  double y = r[1];
2538  double g = 0.;
2539  double gdotb1 = 0.;
2540  double gdotb2 = 0.;
2541  double k5term = 0.;
2542  double k2term = 0.;
2543  double prefactor = 0.;
2544  double v_g = 0.;
2545  double v = (4.*M_PI/A)*epsilon*sigma*sigma*( ((2./5.)*pow((sigma/z),10)) - pow((sigma/z),4) );
2546 
2547  for (double m = -1; m < 1+1; m++) {
2548  for (double n = -1; n < 1+1; n++) {
2549  if ((m != 0) || (n != 0)){
2550  g = sqrt(pow((m*g1x + n*g2x),2.) + pow((m*g1y + n*g2y),2.));
2551  gdotb1 = ((m*g1x + n*g2x)*(b1x+x)) + ((m*g1y + n*g2y)*(b1y+y));
2552  gdotb2 = ((m*g1x + n*g2x)*(b2x+x)) + ((m*g1y + n*g2y)*(b2y+y));
2553  k5term = pow((g*sigma*sigma/2./z),5.)*boost::math::cyl_bessel_k(5., g*z)/30.;
2554  k2term = 2.*pow((g*sigma*sigma/2./z),2.)*boost::math::cyl_bessel_k(2., g*z);
2555  prefactor = epsilon*sigma*sigma*2.*M_PI/A;
2556 
2557  v_g = prefactor*(k5term-k2term);
2558  v += (cos(gdotb1)+cos(gdotb2))*v_g;
2559  }
2560  }
2561  }
2562 
2563  return v;
2564 
2565 }
2566 /**************************************************************************/
2571 Array<dVec,1> GraphenePotential::initialConfig(const Container *boxPtr, MTRand &random,
2572  const int numParticles) {
2573 
2574  /* The particle configuration */
2575  Array<dVec,1> initialPos(numParticles);
2576  initialPos = 0.0;
2577  double initSideCube = 1.0*numParticles;
2578 
2579  for (int i = 0; i < NDIM - 1; i++) {
2580  initSideCube /= boxPtr->side[i];
2581  }
2582 
2583  initSideCube /= ((boxPtr->side[NDIM-1] / 2.0) - 6.0);
2584 
2585  /* Get the linear size per particle, and the number of particles */
2586  double initSide = pow((initSideCube),-1.0/(1.0*NDIM));
2587 
2588  /* We determine the number of initial grid boxes there are in
2589  * in each dimension and compute their size */
2590  int totNumGridBoxes = 1;
2591  iVec numNNGrid;
2592  dVec sizeNNGrid;
2593 
2594  for (int i = 0; i < NDIM - 1; i++) {
2595  numNNGrid[i] = static_cast<int>(ceil((boxPtr->side[i] / initSide) - EPS));
2596 
2597  /* Make sure we have at least one grid box */
2598  if (numNNGrid[i] < 1)
2599  numNNGrid[i] = 1;
2600 
2601  /* Compute the actual size of the grid */
2602  sizeNNGrid[i] = boxPtr->side[i] / (1.0 * numNNGrid[i]);
2603 
2604  /* Determine the total number of grid boxes */
2605  totNumGridBoxes *= numNNGrid[i];
2606  }
2607 
2608  numNNGrid[NDIM-1] = static_cast<int>(ceil(((boxPtr->side[NDIM-1] - 12.0) / (2.0 * initSide)) - EPS));
2609 
2610  /* Make sure we have at least one grid box */
2611  if (numNNGrid[NDIM-1] < 1)
2612  numNNGrid[NDIM-1] = 1;
2613 
2614  /* Compute the actual size of the grid */
2615  sizeNNGrid[NDIM-1] = (boxPtr->side[NDIM-1] - 12.0) / (2.0 * numNNGrid[NDIM-1]);
2616 
2617  /* Determine the total number of grid boxes */
2618  totNumGridBoxes *= numNNGrid[NDIM-1];
2619 
2620  /* Now, we place the particles at the middle of each box */
2621  PIMC_ASSERT(totNumGridBoxes>=numParticles);
2622  dVec pos;
2623  for (int n = 0; n < totNumGridBoxes; n++) {
2624 
2625  iVec gridIndex;
2626  for (int i = 0; i < NDIM; i++) {
2627  int scale = 1;
2628  for (int j = i+1; j < NDIM; j++)
2629  scale *= numNNGrid[j];
2630  gridIndex[i] = (n/scale) % numNNGrid[i];
2631  }
2632 
2633  for (int i = 0; i < NDIM - 1; i++)
2634  pos[i] = (gridIndex[i]+0.5)*sizeNNGrid[i] - 0.5*boxPtr->side[i];
2635 
2636  pos[NDIM - 1] = (gridIndex[NDIM - 1]+0.5)*sizeNNGrid[NDIM - 1] + 3.0;
2637 
2638  boxPtr->putInside(pos);
2639 
2640  if (n < numParticles)
2641  initialPos(n) = pos;
2642  else
2643  break;
2644  }
2645  return initialPos;
2646 }
2647 
2648 // ---------------------------------------------------------------------------
2649 // ---------------------------------------------------------------------------
2650 // GrapheneLUTPotential Class----------------------------------------------------
2651 // ---------------------------------------------------------------------------
2652 // ---------------------------------------------------------------------------
2653 
2654 /**************************************************************************/
2657 #include <boost/math/special_functions/bessel.hpp>
2658 GrapheneLUTPotential::GrapheneLUTPotential (double _strain, double _poisson, double _a0,
2659  double _sigma, double _epsilon, const Container *_boxPtr) : PotentialBase() {
2660 
2661  double strain = _strain;
2662  double poisson = _poisson;
2663  double a0 = _a0;
2664  sigma = _sigma;
2665  epsilon = _epsilon;
2666  V_zmin = 152153.0;
2667 
2668  /* get a local copy of the system size */
2669  Lzo2 = 0.5*_boxPtr->side[NDIM-1];
2670  Lz = _boxPtr->side[NDIM-1];
2671 
2672  /* Arbitrary hard-wall cutoff 1 van der Waals radius (1.4 A) from Lz */
2673  zWall = _boxPtr->side[NDIM-1]-1.4;
2674 
2675  /* Inverse width of the wall onset, corresponding to 1/10 A here. */
2676  invWallWidth = 20.0;
2677 
2678  /* Lookup Tables */
2679  tableLength = int((zmax - zmin)/dr);
2680 
2681  /* gradvg.resize(gtot, tableLength); */
2682 
2683  /* Lattice vectors */
2684  /* @see: https://wiki.cmt.w3.uvm.edu/index.php?title=Bose_Hubbard_model_treatment_for_Helium_absorbed_on_graphene#strain */
2685  a1x = sqrt(3.0)*a0*0.5*(1-strain*poisson);
2686  a1y = 3.0*a0*0.5*(1+strain);
2687  a2x = -a1x;
2688  a2y = a1y;
2689 
2690  /* reciprocal lattice vectors */
2691  g1x = 2*M_PI/(sqrt(3.0)*a0*(1-strain*poisson));
2692  g1y = 2*M_PI/(3.0*a0*(1+strain));
2693  g2x = -g1x;
2694  g2y = g1y;
2695 
2696  /* basis vectors */
2697  b1x = sqrt(3.0)*a0*0.5*(1.0-strain*poisson);
2698  b1y = 0.5*a0*(1.0 + strain);
2699  b2x = 0.0;
2700  b2y = a0*(1.0 + strain);
2701 
2702  /* area of unit cell */
2703  A = fabs((a1x*a2y) - (a1y*a2x));
2704 
2705  double g = 0.0;
2706  double prefactor = epsilon*sigma*sigma*2.*M_PI/A;
2707  double k5term = 0.0;
2708  double k2term = 0.0;
2709  /* double dk1term = 0.0; */
2710  /* double dk2term = 0.0; */
2711  /* double dk3term = 0.0; */
2712  /* double dk4term = 0.0; */
2713  /* double dk5term = 0.0; */
2714  /* double dk6term = 0.0; */
2715  double z = 0.0;
2716 
2717  /* FULL LOOKUP TABLE OFFERS NO SPEED BENEFITS */
2718  /* vg.resize(2*(gnum+1)*gnum+1,tableLength); */
2719 
2720  /* for (int iz = 0; iz < tableLength; iz++) { */
2721  /* z = zmin + iz*dr; */
2722  /* vg(0,iz) = q*prefactor*( ((2./5.)*pow((sigma/z),10)) - pow((sigma/z),4) ); */
2723  /* } */
2724 
2725  /* int ig = 1; */
2726  /* for (int m = -gnum; m <= gnum; m++) { */
2727  /* for (int n = 1; n <= gnum; n++) { */
2728  /* g = sqrt(pow((m*g1x + n*g2x),2) + pow((m*g1y + n*g2y),2)); */
2729  /* for (int iz = 0; iz < tableLength; iz++) { */
2730  /* z = zmin + iz*dr; */
2731  /* k5term = pow((g*sigma*sigma/2./z),5)*boost::math::cyl_bessel_k(5, g*z)/30.; */
2732  /* k2term = 2.*pow((g*sigma*sigma/2./z),2)*boost::math::cyl_bessel_k(2, g*z); */
2733  /* vg(ig,iz) = prefactor*(k5term-k2term); */
2734  /* } */
2735  /* ig++; */
2736  /* } */
2737  /* } */
2738  /* for (int m=1; m <= gnum; m++) { */
2739  /* g = sqrt(pow(m*g1x,2) + pow(m*g1y,2)); */
2740  /* for (int iz = 0; iz < tableLength; iz++) { */
2741  /* z = zmin + iz*dr; */
2742  /* k5term = pow((g*sigma*sigma/2./z),5)*boost::math::cyl_bessel_k(5, g*z)/30.; */
2743  /* k2term = 2.*pow((g*sigma*sigma/2./z),2)*boost::math::cyl_bessel_k(2, g*z); */
2744  /* vg(ig,iz) = prefactor*(k5term-k2term); */
2745  /* } */
2746  /* ig++; */
2747  /* } */
2748 
2749  /* Create a unique list of all g-vector magnitudes */
2750  set<double> uniquegMag;
2751 
2752  /* This gives us the upper half-plane */
2753  uniquegMag.insert(0.0);
2754  for (int m = -gnum; m <= gnum; m++) {
2755  for (int n = 1; n <= gnum; n++) {
2756  g = sqrt(pow((m*g1x + n*g2x),2) + pow((m*g1y + n*g2y),2));
2757  uniquegMag.insert(g);
2758  }
2759  }
2760  for (int m=1; m <= gnum; m++) {
2761  g = sqrt(pow(m*g1x,2) + pow(m*g1y,2));
2762  uniquegMag.insert(g);
2763  }
2764 
2765  /* Convert the set to a vector and sort */
2766  std::vector<double> gMag(uniquegMag.begin(), uniquegMag.end());
2767  sort(gMag.begin(),gMag.end());
2768 
2769  /* Create the mapping from g-vectors to g-magnitudes */
2770  gMagID.resize(2*(gnum+1)*gnum+1);
2771  gMagID(0) = 0;
2772 
2773  int ig = 1;
2774  for (int m = -gnum; m <= gnum; m++) {
2775  for (int n = 1; n <= gnum; n++) {
2776  g = sqrt(pow((m*g1x + n*g2x),2) + pow((m*g1y + n*g2y),2));
2777  gMagID(ig) = distance(gMag.begin(), find(gMag.begin(), gMag.end(),g));
2778  ig++;
2779  }
2780  }
2781  for (int m=1; m <= gnum; m++){
2782  g = sqrt(pow(m*g1x,2) + pow(m*g1y,2));
2783  gMagID(ig) = distance(gMag.begin(), find(gMag.begin(), gMag.end(),g));
2784  ig++;
2785  }
2786 
2787  /* Fill up the lookup table based on g-magnitudes only */
2788  vg.resize(gMag.size(), tableLength);
2789 
2790  ig = 0;
2791  for (auto cg : gMag) {
2792  for (int iz = 0; iz < tableLength; iz++) {
2793  z = zmin + iz*dr;
2794  if (ig == 0)
2795  vg(ig,iz) = q*prefactor*( ((2./5.)*pow((sigma/z),10)) - pow((sigma/z),4) );
2796  else {
2797  k5term = pow((cg*sigma*sigma/2./z),5)*boost::math::cyl_bessel_k(5, cg*z)/30.;
2798  k2term = 2.*pow((cg*sigma*sigma/2./z),2)*boost::math::cyl_bessel_k(2, cg*z);
2799  vg(ig,iz) = prefactor*(k5term-k2term);
2800  }
2801  }
2802  ig++;
2803  }
2804 
2805  /* print out the potential to disk */
2806  /* int numPoints = 500; */
2807  /* double dx = _boxPtr->side[0]/numPoints; */
2808  /* double dy = _boxPtr->side[1]/numPoints; */
2809  /* dVec pos; */
2810  /* pos[2] = -0.5*_boxPtr->side[2] + 2.635; */
2811  /* for (int i = 0; i < numPoints; i++) { */
2812  /* pos[0] = -0.5*_boxPtr->side[0] + i*dx; */
2813  /* for (int j = 0; j < numPoints; j++) { */
2814  /* pos[1] = -0.5*_boxPtr->side[1] + j*dy; */
2815  /* communicate()->file("debug")->stream() << format("%24.16e %24.16e %24.16e\n") % pos[0] % pos[1] % V(pos); */
2816  /* } */
2817  /* } */
2818 
2819  /* exit(-1); */
2820 
2821  /* print out the potential to disk */
2822  /* int numPoints = 500000; */
2823  /* double dz = 2*_boxPtr->side[2]/numPoints; */
2824  /* dVec pos; */
2825  /* pos[0] = 0.0; */
2826  /* pos[1] = 0.0; */
2827  /* for (int i = 0; i < numPoints; i++) { */
2828  /* pos[2] = -0.5*_boxPtr->side[2] + 2.0 + i*dz; */
2829  /* double cz = pos[2] + 0.5*_boxPtr->side[2]; */
2830  /* communicate()->file("debug")->stream() << format("%24.16e %24.16e\n") % cz % V(pos); */
2831  /* } */
2832 
2833  /* exit(-1); */
2834 }
2835 
2836 
2837 /**************************************************************************/
2841  gMagID.free();
2842  vg.free();
2843  gradvg.free();
2844 }
2845 
2846 /**************************************************************************/
2853 double GrapheneLUTPotential::V(const dVec &r) {
2854 
2855  double z = r[2]+(Lzo2);
2856  if (z < zmin)
2857  return V_zmin;
2858 
2859  /* A sigmoid to represent the hard-wall */
2860  double VWall = V_zmin/(1.0+exp(-invWallWidth*(z-zWall)));
2861 
2862  int zindex = int((z-zmin)/dr);
2863  if (zindex >= tableLength)
2864  return q*(epsilon*sigma*sigma*2.*M_PI/A)*( ((2./5.)*pow((sigma/z),10)) - pow((sigma/z),4) );
2865 
2866  double x1 = r[0]+b1x;
2867  double x2 = r[0]+b2x;
2868  double y1 = r[1]+b1y;
2869  double y2 = r[1]+b2y;
2870 
2871  double mx,my;
2872  double v = vg(0,zindex);
2873 
2874  int ig = 1;
2875  for (int m = -gnum; m <= gnum; m++) {
2876  mx = m*g1x;
2877  my = m*g1y;
2878  for (int n = 1; n <= gnum; n++) {
2879  v += 2.0*(cos((mx + n*g2x)*x1 + (my + n*g2y)*y1) +
2880  cos((mx + n*g2x)*x2 + (my + n*g2y)*y2)) * vg(gMagID(ig),zindex);
2881  ig++;
2882  }
2883  }
2884 
2885  for (int m=1; m <= gnum; m++){
2886  mx = m*g1x;
2887  my = m*g1y;
2888  v += 2.0*(cos(mx*x1 + my*y1) + cos(mx*x2 + my*y2)) * vg(gMagID(ig),zindex);
2889  ig++;
2890  }
2891 
2892  return v + VWall;
2893 }
2894 
2895 /**************************************************************************/
2903 
2904  return 0.0;
2905 
2906  /* double x = r[0]; */
2907  /* double y = r[1]; */
2908  /* double z = r[NDIM-1]+(Lzo2); */
2909  /* int zindex = int((z-zmin)/dr); */
2910  /* dVec dv = 0.0; */
2911 
2912  /* if (z < zmin) */
2913  /* return dv; */
2914 
2915  /* if (zindex >= tableLength) { */
2916  /* dv[NDIM-1] += q*(epsilon*sigma*sigma*2.*M_PI/A)* 4. * (pow(sigma/z,4) - pow(sigma/z,10)) / z; */
2917  /* return dv; */
2918  /* } */
2919 
2920  /* dv[NDIM-1] = gradvg(0,zindex); */
2921  /* double gdotb1 = 0.; */
2922  /* double gdotb2 = 0.; */
2923 
2924  /* int k = 0; */
2925 
2926  /* /1* NB: this has not been optimized!!! *1/ */
2927  /* for (int m = -gnum; m < gnum+1; m++) { */
2928  /* for (int n = -gnum; n < gnum+1; n++) { */
2929  /* k = karr(m+gnum,n+gnum); */
2930  /* if((m != 0) or (n != 0)) { */
2931  /* gdotb1 = ((m*g1x + n*g2x)*(b1x+x)) + ((m*g1y + n*g2y)*(b1y+y)); */
2932  /* gdotb2 = ((m*g1x + n*g2x)*(b2x+x)) + ((m*g1y + n*g2y)*(b2y+y)); */
2933 
2934  /* dv[0] += -(g1x*m+g2x*n) * (sin(gdotb1) + sin(gdotb2))*vg(k,zindex); */
2935  /* dv[1] += -(g1y*m+g2y*n) * (sin(gdotb1) + sin(gdotb2))*vg(k,zindex); */
2936  /* dv[NDIM-1] += (cos(gdotb1)+cos(gdotb2))*gradvg(k,zindex); */
2937  /* } */
2938  /* } */
2939  /* } */
2940  /* return dv; */
2941 }
2942 
2943 /**************************************************************************/
2948 Array<dVec,1> GrapheneLUTPotential::initialConfig(const Container *boxPtr, MTRand &random,
2949  const int numParticles) {
2950  /* The particle configuration */
2951  Array<dVec,1> initialPos(numParticles);
2952  initialPos = 0.0;
2953 
2954  int Nlayer = round(boxPtr->side[0]*boxPtr->side[1]/a1x/a1y/4);
2955  int numX = round(boxPtr->side[0]/a1x/2);
2956  int numY = round(boxPtr->side[1]/a1y/2);
2957  double gridX = boxPtr->side[0]/numX;
2958  double gridY = boxPtr->side[1]/numY;
2959  int gridSize = 0;
2960  if (3*Nlayer < numParticles){
2961  gridSize = numParticles - 3*Nlayer;
2962  }
2963  int numInLayers = numParticles - gridSize;
2964  int layerNum = 0;
2965  int iShifted = 0;
2966  dVec pos;
2967  for (int i = 0; i < numInLayers; i++){
2968  layerNum = i/Nlayer;
2969  iShifted = i - (layerNum*Nlayer);
2970  pos[0] = ((iShifted % numX) + 0.5)*gridX - (boxPtr->side[0]/2);
2971  pos[1] = ((iShifted / numX) + 0.5)*gridY - (boxPtr->side[1]/2);
2972  pos[NDIM-1] = ((layerNum+1)*3.0)-(boxPtr->side[NDIM-1]/2);
2973  boxPtr->putInside (pos);
2974  initialPos(i) = pos;
2975  }
2976  if (gridSize) {
2977  int initCubePart = ceil(pow(1.0*gridSize,1.0/(1.0*NDIM)));
2978  dVec initSide;
2979  for (int i = 0; i < NDIM - 1; i++) {
2980  initSide[i] = boxPtr->side[i]/initCubePart;
2981  }
2982 
2983  initSide[NDIM-1] = (boxPtr->side[NDIM-1] - 12.0)/initCubePart;
2984 
2985  for (int n = 0; n < gridSize; n++) {
2986  pos[0] = (((n % initCubePart) + 0.5) * initSide[0]) - (boxPtr->side[0]/2);
2987  pos[1] = (((n / initCubePart) + 0.5) * initSide[1]) - (boxPtr->side[1]/2);
2988  pos[NDIM-1] = (((n / initCubePart / initCubePart) + 0.5) * initSide[NDIM-1]) - (boxPtr->side[NDIM-1]/2) + 12.0;
2989  boxPtr->putInside (pos);
2990  initialPos(n + numInLayers) = pos;
2991  }
2992  }
2993  return initialPos;
2994 }
2995 
2996 
2997 // ---------------------------------------------------------------------------
2998 // ---------------------------------------------------------------------------
2999 // GrapheneLUT3DPotential Class---------------------------------------------
3000 // ---------------------------------------------------------------------------
3001 // ---------------------------------------------------------------------------
3002 
3003 /**************************************************************************/
3006 GrapheneLUT3DPotential::GrapheneLUT3DPotential (string graphenelut3d_file_prefix, const Container *_boxPtr) : PotentialBase() {
3007 
3008  static auto const aflags = boost::archive::no_header | boost::archive::no_tracking;
3009  /* get a local copy of the system size */
3010  Lzo2 = _boxPtr->side[NDIM-1]/2;
3011 
3012  /* Arbitrary hard-wall cutoff 1 van der Waals radius (1.4 A) from Lz */
3013  zWall = _boxPtr->side[NDIM-1]-1.4;
3014 
3015  /* Inverse width of the wall onset, corresponding to 1/10 A here. */
3016  invWallWidth = 20.0;
3017 
3018  /* load lookup tables */
3019  {
3020  // create and open a character archive for input
3021  std::ifstream ifs(graphenelut3d_file_prefix + std::string("serialized.dat"));
3022  boost::archive::binary_iarchive ia(ifs,aflags);
3023  // write class instance to archive
3024  ia >> V3d >> gradV3d_x >> gradV3d_y >> gradV3d_z >> grad2V3d >> LUTinfo;
3025  // archive and stream closed when destructors are called
3026  }
3027 
3028  A11 = LUTinfo( 0 );
3029  A12 = LUTinfo( 1 );
3030  A21 = LUTinfo( 2 );
3031  A22 = LUTinfo( 3 );
3032  dx = LUTinfo( 4 );
3033  dy = LUTinfo( 5 );
3034  dz = LUTinfo( 6 );
3035  cell_length_a = LUTinfo( 7 );
3036  cell_length_b = LUTinfo( 8 );
3037  zmin = LUTinfo( 9 );
3038  zmax = LUTinfo( 10 );
3039  V_zmin = LUTinfo( 11 );
3040 }
3041 
3042 
3043 /**************************************************************************/
3047  V3d.free(); // Potential lookup table
3048  gradV3d_x.free(); // gradient of potential x direction lookup table
3049  gradV3d_y.free(); // gradient of potential y direction lookup table
3050  gradV3d_z.free(); // gradient of potential z direction lookup table
3051  grad2V3d.free(); // Laplacian of potential
3052  LUTinfo.free(); // Information about the 3D lookup table
3053 }
3054 
3055 /**************************************************************************/
3063  double z = r[2]+Lzo2;
3064 
3065  /* Check for extremal values of z */
3066  if (z < zmin)
3067  return V_zmin;
3068 
3069  /* A sigmoid to represent the hard-wall */
3070  double VWall = V_zmin/(1.0+exp(-invWallWidth*(z-zWall)));
3071 
3072  dVec _r = r;
3073  _r[2] += Lzo2 - zmin;
3074 
3075  cartesian_to_uc(_r, A11, A12, A21, A22);
3076  put_in_uc(_r, cell_length_a, cell_length_b);
3077  double _V = trilinear_interpolation(V3d, _r, dx, dy, dz);
3078  //double _V = direct_lookup(V3d, _r, dx, dy, dz);
3079  //std::cout << _r << " " << _V << std::endl;
3080 
3081  return _V + VWall;
3082 }
3083 
3084 /**************************************************************************/
3092  dVec _gradV = 0.0;
3093  if (r[2] + Lzo2 < zmin){
3094  return _gradV;
3095  }
3096  dVec _r = r;
3097  _r[2] += Lzo2 - zmin;
3098 
3099  cartesian_to_uc(_r, A11, A12, A21, A22);
3100  put_in_uc(_r, cell_length_a, cell_length_b);
3101  double _gradV_x = trilinear_interpolation(gradV3d_x, _r, dx, dy, dz);
3102  double _gradV_y = trilinear_interpolation(gradV3d_y, _r, dx, dy, dz);
3103  double _gradV_z = trilinear_interpolation(gradV3d_z, _r, dx, dy, dz);
3104  //double _gradV_x = direct_lookup(gradV3d_x, _r, dx, dy, dz);
3105  //double _gradV_y = direct_lookup(gradV3d_y, _r, dx, dy, dz);
3106  //double _gradV_z = direct_lookup(gradV3d_z, _r, dx, dy, dz);
3107  _gradV(0) = _gradV_x;
3108  _gradV(1) = _gradV_y;
3109  _gradV(2) = _gradV_z;
3110  return _gradV;
3111 }
3112 
3113 /**************************************************************************/
3121  double _grad2V = 0.0;
3122  if (r[2] + Lzo2 < zmin){
3123  return 0.0;
3124  }
3125  dVec _r = r;
3126  _r[2] += Lzo2 - zmin;
3127 
3128  cartesian_to_uc(_r, A11, A12, A21, A22);
3129  put_in_uc(_r, cell_length_a, cell_length_b);
3130  _grad2V = trilinear_interpolation(grad2V3d, _r, dx, dy, dz);
3131  //_grad2V = direct_lookup(grad2V3d, _r, dx, dy, dz);
3132  return _grad2V;
3133 }
3134 
3135 /**************************************************************************/
3140 Array<dVec,1> GrapheneLUT3DPotential::initialConfig1(const Container *boxPtr, MTRand &random,
3141  const int numParticles) {
3142  //FIXME We initialize on grid here above zmin, should I do C_1/3 here. I think not because might
3143  // be hard to get out of if other configurations are more favorable
3144 
3145  /* The particle configuration */
3146  Array<dVec,1> initialPos(numParticles);
3147  initialPos = 0.0;
3148 
3149  /* Get the average inter-particle separation in the accessible volume. */
3150  double initSide = pow((1.0*numParticles/boxPtr->volume),-1.0/(1.0*NDIM));
3151 
3152  /* We determine the number of initial grid boxes there are in
3153  * in each dimension and compute their size */
3154  int totNumGridBoxes = 1;
3155  iVec numNNGrid;
3156  dVec sizeNNGrid;
3157 
3158  for (int i = 0; i < NDIM; i++) {
3159  numNNGrid[i] = static_cast<int>(ceil((boxPtr->side[i] / initSide) - EPS));
3160 
3161  /* Make sure we have at least one grid box */
3162  if (numNNGrid[i] < 1)
3163  numNNGrid[i] = 1;
3164 
3165  /* Compute the actual size of the grid */
3166  sizeNNGrid[i] = boxPtr->side[i] / (1.0 * numNNGrid[i]);
3167 
3168  /* Modify due to the graphene and hard wall */
3169  if (i == NDIM - 1) {
3170  sizeNNGrid[i] = (zWall - zmin) / (1.0 * numNNGrid[i]);
3171  }
3172 
3173  /* Determine the total number of grid boxes */
3174  totNumGridBoxes *= numNNGrid[i];
3175  }
3176 
3177  /* Now, we place the particles at the middle of each box */
3178  PIMC_ASSERT(totNumGridBoxes>=numParticles);
3179  dVec pos;
3180  for (int n = 0; n < totNumGridBoxes; n++) {
3181 
3182  iVec gridIndex;
3183  for (int i = 0; i < NDIM; i++) {
3184  int scale = 1;
3185  for (int j = i+1; j < NDIM; j++)
3186  scale *= numNNGrid[j];
3187  gridIndex[i] = (n/scale) % numNNGrid[i];
3188  }
3189 
3190  /* We treat the z-direction differently due to the graphene and wall */
3191  for (int i = 0; i < NDIM; i++)
3192  pos[i] = (gridIndex[i]+0.5)*sizeNNGrid[i] - 0.5*boxPtr->side[i];
3193  pos[NDIM-1] += zmin;
3194 
3195  boxPtr->putInside(pos);
3196 
3197  if (n < numParticles)
3198  initialPos(n) = pos;
3199  else
3200  break;
3201  }
3202 
3203  /* Output for plotting in python */
3204  /* cout << "np.array(["; */
3205  /* for (int n = 0; n < numParticles; n++) { */
3206  /* cout << "["; */
3207  /* for (int i = 0; i < NDIM-1; i++) */
3208  /* cout << initialPos(n)[i] << ", "; */
3209  /* cout << initialPos(n)[NDIM-1] << "]," << endl; */
3210  /* } */
3211  /* cout << "])" << endl; */
3212  /* exit(-1); */
3213 
3214  return initialPos;
3215 }
3216 
3217 /**************************************************************************/
3222 Array<dVec,1> GrapheneLUT3DPotential::initialConfig(const Container *boxPtr, MTRand &random,
3223  const int numParticles) {
3224 
3225  /* The particle configuration */
3226  Array<dVec,1> initialPos(numParticles);
3227  initialPos = 0.0;
3228 
3229  dVec side_;
3230  side_ = boxPtr->side;
3231 
3232  /* The first particle is randomly placed inside the box */
3233  for (int i = 0; i < NDIM-1; i++)
3234  initialPos(0)[i] = side_[i]*(-0.5 + random.rand());
3235  initialPos(0)[NDIM-1] = -0.5*side_[NDIM-1]+zmin + (zWall-zmin)*random.rand();
3236 
3237  double hardCoreRadius = 1.00;
3238 
3239  /* We keep trying to insert particles with no overlap */
3240  int numInserted = 1;
3241  int numAttempts = 0;
3242  int maxNumAttempts = ipow(10,6);
3243  dVec trialPos,sep;
3244  double rsep;
3245  while ((numInserted < numParticles) && (numAttempts < maxNumAttempts)) {
3246  numAttempts += 1;
3247 
3248  /* Generate a random position inside the box */
3249  for (int i = 0; i < NDIM-1; i++)
3250  trialPos[i] = side_[i]*(-0.5 + random.rand());
3251  trialPos[NDIM-1] = -0.5*side_[NDIM-1]+zmin + (zWall-zmin)*random.rand();
3252 
3253  /* Make sure it doesn't overlap with any existing particles. If not,
3254  * add it, otherwise, break */
3255  bool acceptParticle = true;
3256  for (int n = 0; n < numInserted; n++) {
3257  sep = trialPos - initialPos(n);
3258  boxPtr->putInBC(sep);
3259  rsep = sqrt(dot(sep,sep));
3260  if (rsep < 2*hardCoreRadius) {
3261  acceptParticle = false;
3262  break;
3263  }
3264  }
3265 
3266  if (acceptParticle) {
3267  initialPos(numInserted) = trialPos;
3268  numInserted += 1;
3269  }
3270 
3271  }
3272 
3273  if (numAttempts == maxNumAttempts) {
3274  cerr << "Could not construct a valid initial configuration." << endl;
3275  exit(EXIT_FAILURE);
3276  }
3277 
3278  /* Output configuration to terminal suitable for plotting with python */
3279  /* cout << "np.array(["; */
3280  /* for (int n = 0; n < numParticles; n++) { */
3281  /* cout << "["; */
3282  /* for (int i = 0; i < NDIM-1; i++) */
3283  /* cout << initialPos(n)[i] << ", "; */
3284  /* cout << initialPos(n)[NDIM-1] << "]," << endl; */
3285  /* } */
3286  /* cout << "])" << endl; */
3287  /* exit(EXIT_FAILURE); */
3288 
3289  return initialPos;
3290 }
3291 
3292 
3293 // ---------------------------------------------------------------------------
3294 // ---------------------------------------------------------------------------
3295 // GrapheneLUT3DPotentialGenerate Class---------------------------------------
3296 // ---------------------------------------------------------------------------
3297 // ---------------------------------------------------------------------------
3298 
3299 /**************************************************************************/
3303  const double _strain, const double _poisson_ratio,
3304  const double _carbon_carbon_distance, const double _sigma,
3305  const double _epsilon, const Container *_boxPtr
3306  ) : PotentialBase() {
3307 
3308  static auto const aflags = boost::archive::no_header | boost::archive::no_tracking;
3309  // ADD FUNCTIONS TO CONVERT FROM BINARY TO TEXT AND VICE VERSA FOR SERIALIZED FILES
3310  /* get a local copy of the system size */
3311  Lzo2 = _boxPtr->side[NDIM-1]/2;
3312 
3313  //double sigma_0 = 2.74;
3314  //double epsilon_0 = 16.2463;
3315 
3316  double strain = _strain;
3317  /* double poisson_ratio = _poisson_ratio; */
3318  /* double carbon_carbon_distance = _carbon_carbon_distance; */
3319  double sigma = _sigma;
3320  double epsilon = _epsilon;
3321 
3322  /* The hard-coded resolution of the lookup table and maximum distance
3323  * above the sheet */
3324  xres = 101;
3325  yres = 101;
3326  zmax = 10.0;
3327  zres = 1001;
3328 
3329  /* xres = 5; */
3330  /* yres = 5; */
3331  /* zmax = 10.0; */
3332  /* zres = 51; */
3333 
3334  auto [ V3d, gradV3d_x, gradV3d_y, gradV3d_z, grad2V3d, xy_x, xy_y, LUTinfo ] =
3335  get_V3D_all(strain,sigma,epsilon,xres,yres,zres,zmax);
3336 
3337  string graphenelut3d_file_prefix = str( format( "graphene_%.2f_%.2f_%d_%d_%d_") %
3338  strain % zmax % xres % yres % zres );
3339 
3340  // create and open a character archive for output
3341  std::ofstream ofs(graphenelut3d_file_prefix + "serialized.dat");
3342 
3343  // save data to archive
3344  {
3345  boost::archive::binary_oarchive oa(ofs,aflags);
3346  // write class instance to archive
3347  oa << V3d << gradV3d_x << gradV3d_y << gradV3d_z << grad2V3d << LUTinfo;
3348  // archive and stream closed when destructors are called
3349  }
3350 }
3351 
3352 double GrapheneLUT3DPotentialGenerate::Vz_64(
3353  double z, double sigma, double prefactor, int atoms_in_basis ) {
3354  return atoms_in_basis*Vz_64(z,sigma,prefactor);
3355 }
3356 
3357 double GrapheneLUT3DPotentialGenerate::Vz_64(
3358  double z, double sigma, double prefactor ) {
3359  return prefactor*Vz_64(z,sigma);
3360 }
3361 
3362 double GrapheneLUT3DPotentialGenerate::Vz_64( double z, double sigma ) {
3363  return (pow(sigma,4)*(2*pow(sigma,6) - 5*pow(z,6)))/(5*pow(z,10));
3364 }
3365 
3366 double GrapheneLUT3DPotentialGenerate::gradVz_x_64(
3367  double z, double sigma, double prefactor, int atoms_in_basis ) {
3368  return 0.0;
3369 }
3370 
3371 double GrapheneLUT3DPotentialGenerate::gradVz_x_64(
3372  double z, double sigma, double prefactor ) {
3373  return 0.0;
3374 }
3375 
3376 double GrapheneLUT3DPotentialGenerate::gradVz_x_64( double z, double sigma ) {
3377  return 0.0;
3378 }
3379 
3380 double GrapheneLUT3DPotentialGenerate::gradVz_y_64(
3381  double z, double sigma, double prefactor, int atoms_in_basis ) {
3382  return 0.0;
3383 }
3384 
3385 double GrapheneLUT3DPotentialGenerate::gradVz_y_64(
3386  double z, double sigma, double prefactor ) {
3387  return 0.0;
3388 }
3389 
3390 double GrapheneLUT3DPotentialGenerate::gradVz_y_64( double z, double sigma ) {
3391  return 0.0;
3392 }
3393 
3394 double GrapheneLUT3DPotentialGenerate::gradVz_z_64(
3395  double z, double sigma, double prefactor, int atoms_in_basis ) {
3396  return atoms_in_basis*gradVz_z_64(z,sigma,prefactor);
3397 }
3398 
3399 double GrapheneLUT3DPotentialGenerate::gradVz_z_64(
3400  double z, double sigma, double prefactor ) {
3401  return prefactor*gradVz_z_64(z,sigma);
3402 }
3403 
3404 double GrapheneLUT3DPotentialGenerate::gradVz_z_64( double z, double sigma ) {
3405  return (4*pow(sigma,4)*(-pow(sigma,6) + pow(z,6)))/pow(z,11);
3406 }
3407 
3408 double GrapheneLUT3DPotentialGenerate::grad2Vz_64(
3409  double z, double sigma, double prefactor, int atoms_in_basis ) {
3410  return atoms_in_basis*grad2Vz_64(z,sigma,prefactor);
3411 }
3412 
3413 double GrapheneLUT3DPotentialGenerate::grad2Vz_64(
3414  double z, double sigma, double prefactor ) {
3415  return prefactor*grad2Vz_64(z,sigma);
3416 }
3417 
3418 double GrapheneLUT3DPotentialGenerate::grad2Vz_64( double z, double sigma ) {
3419  return (4*pow(sigma,4)*(11*pow(sigma,6) - 5*pow(z,6)))/pow(z,12);
3420 }
3421 
3422 double GrapheneLUT3DPotentialGenerate::Vg_64(
3423  double x, double y, double z, double sigma, double g, double g_i,
3424  double g_j, double b_i, double b_j, double prefactor ) {
3425  return prefactor*Vg_64(x, y, z, sigma, g, g_i, g_j, b_i, b_j);
3426 }
3427 
3428 double GrapheneLUT3DPotentialGenerate::Vg_64(
3429  double x, double y, double z, double sigma, double g, double g_i,
3430  double g_j, double b_i, double b_j ) {
3431  return (pow(g,2)*pow(sigma,4)*(-480*pow(z,3)*boost::math::cyl_bessel_k(2, g*z) +
3432  pow(g,3)*pow(sigma,6)*boost::math::cyl_bessel_k(5, g*z))*cos(g_i*(b_i + x) +
3433  g_j*(b_j + y)))/(960*pow(z,5));
3434 }
3435 
3436 double GrapheneLUT3DPotentialGenerate::gradVg_x_64(
3437  double x, double y, double z, double sigma, double g, double g_i,
3438  double g_j, double b_i, double b_j, double prefactor ) {
3439  return prefactor*gradVg_x_64(x, y, z, sigma, g, g_i, g_j, b_i, b_j);
3440 }
3441 
3442 double GrapheneLUT3DPotentialGenerate::gradVg_x_64(
3443  double x, double y, double z, double sigma, double g, double g_i,
3444  double g_j, double b_i, double b_j ) {
3445  return -(pow(g,2)*g_i*pow(sigma,4)*(-480*pow(z,3)*boost::math::cyl_bessel_k(2, g*z) +
3446  pow(g,3)*pow(sigma,6)*boost::math::cyl_bessel_k(5, g*z))*sin(g_i*(b_i + x) +
3447  g_j*(b_j + y)))/(960*pow(z,5));
3448 }
3449 
3450 double GrapheneLUT3DPotentialGenerate::gradVg_y_64(
3451  double x, double y, double z, double sigma, double g, double g_i,
3452  double g_j, double b_i, double b_j, double prefactor ) {
3453  return prefactor*gradVg_y_64(x, y, z, sigma, g, g_i, g_j, b_i, b_j);
3454 }
3455 
3456 double GrapheneLUT3DPotentialGenerate::gradVg_y_64(
3457  double x, double y, double z, double sigma, double g, double g_i,
3458  double g_j, double b_i, double b_j ) {
3459  return -(pow(g,2)*g_j*pow(sigma,4)*(-480*pow(z,3)*boost::math::cyl_bessel_k(2, g*z) +
3460  pow(g,3)*pow(sigma,6)*boost::math::cyl_bessel_k(5, g*z))*sin(g_i*(b_i + x) +
3461  g_j*(b_j + y)))/(960*pow(z,5));
3462 }
3463 
3464 double GrapheneLUT3DPotentialGenerate::gradVg_z_64(
3465  double x, double y, double z, double sigma, double g, double g_i,
3466  double g_j, double b_i, double b_j, double prefactor ) {
3467  return prefactor*gradVg_z_64(x, y, z, sigma, g, g_i, g_j, b_i, b_j);
3468 }
3469 
3470 double GrapheneLUT3DPotentialGenerate::gradVg_z_64(
3471  double x, double y, double z, double sigma, double g, double g_i,
3472  double g_j, double b_i, double b_j ) {
3473  return -(pow(g,3)*pow(sigma,4)*(10*(pow(g,2)*pow(sigma,6)*z - 48*pow(z,5))*boost::math::cyl_bessel_k(3, g*z) +
3474  g*pow(sigma,6)*(80 + pow(g*z,2))*boost::math::cyl_bessel_k(4, g*z))*cos(g_i*(b_i + x) +
3475  g_j*(b_j + y)))/(960*pow(z,7));
3476 }
3477 
3478 double GrapheneLUT3DPotentialGenerate::grad2Vg_64(
3479  double x, double y, double z, double sigma, double g, double g_i,
3480  double g_j, double b_i, double b_j, double prefactor ) {
3481  return prefactor*grad2Vg_64(x, y, z, sigma, g, g_i, g_j, b_i, b_j);
3482 }
3483 
3484 double GrapheneLUT3DPotentialGenerate::grad2Vg_64(
3485  double x, double y, double z, double sigma, double g, double g_i,
3486  double g_j, double b_i, double b_j) {
3487  return (pow(g,2)*pow(sigma,4)*(-120*g*pow(z,6)*(g*z*boost::math::cyl_bessel_k(0, g*z) +
3488  8*boost::math::cyl_bessel_k(1, g*z)) + z*(19*pow(g,4)*pow(sigma,6)*pow(z,2) +
3489  480*pow(z,4)*(-6 + (pow(g_i,2) + pow(g_j,2))*pow(z,2)) -
3490  8*pow(g,2)*(45*pow(z,6) + pow(sigma,6)*(-110 + (pow(g_i,2) +
3491  pow(g_j,2))*pow(z,2))))*boost::math::cyl_bessel_k(2, g*z) +
3492  g*(-1680*pow(z,6) + pow(sigma,6)*(5280 + pow(z,2)*(-48*(pow(g_i,2) +
3493  pow(g_j,2)) + pow(g,4)*pow(z,2) + pow(g,2)*(224 -
3494  (pow(g_i,2) + pow(g_j,2))*pow(z,2)))))*boost::math::cyl_bessel_k(3, g*z))*cos(g_i*(b_i + x) +
3495  g_j*(b_j + y)))/(960*pow(z,9));
3496 }
3497 
3498 double GrapheneLUT3DPotentialGenerate::V_64(
3499  double x, double y, double z, double sigma, double epsilon,
3500  double area_lattice, TinyVector<double,2> b_1, TinyVector<double,2> b_2,
3501  TinyVector<double,2> g_m, TinyVector<double,2> g_n, Array<int,1> g_i_array,
3502  Array<int,1> g_j_array, Array<double,1> g_magnitude_array ) {
3503  int k_max=10000;
3504  bool flag_1 = false;
3505  bool flag_2 = false;
3506  bool flag_3 = false;
3507  TinyVector <double,2> g;
3508  double pf = 2*M_PI*epsilon*pow(sigma,2)/area_lattice;
3509  double _V = 0.0;
3510  _V += 2*Vz_64(z, sigma);
3511  double _V_old = 0.0;
3512  int g_i;
3513  int g_j;
3514  double g_magnitude;
3515  for (int k = 1; k < 2*k_max; k++) {
3516  g_i = g_i_array(k);
3517  g_j = g_j_array(k);
3518  g = (g_i * g_m) + (g_j * g_n);
3519  g_magnitude = g_magnitude_array(k);
3520  _V += Vg_64(x, y, z, sigma, g_magnitude, g[0], g[1], b_1[0], b_1[1]);
3521  _V += Vg_64(x, y, z, sigma, g_magnitude, g[0], g[1], b_2[0], b_2[1]);
3522 
3523  if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_1)) {
3524  _V_old = _V;
3525  flag_1 = true;
3526  }
3527  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_2)) {
3528  _V_old = _V;
3529  flag_2 = true;
3530  }
3531  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_3)) {
3532  _V_old = _V;
3533  flag_3 = true;
3534  }
3535  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && flag_1 && flag_2 && flag_3) {
3536  _V_old = _V;
3537  break;
3538  }
3539  else if ((k >= k_max) && (g_magnitude != g_magnitude_array(k+1))) {
3540  //println("WARNING: k_max reached")
3541  break;
3542  }
3543  else {
3544  _V_old = _V;
3545  }
3546  }
3547  return pf * _V;
3548 }
3549 
3550 double GrapheneLUT3DPotentialGenerate::gradV_x_64(
3551  double x, double y, double z, double sigma, double epsilon,
3552  double area_lattice, TinyVector<double,2> b_1, TinyVector<double,2> b_2,
3553  TinyVector<double,2> g_m, TinyVector<double,2> g_n, Array<int,1> g_i_array,
3554  Array<int,1> g_j_array, Array<double,1> g_magnitude_array ) {
3555  int k_max=10000;
3556  bool flag_1 = false;
3557  bool flag_2 = false;
3558  bool flag_3 = false;
3559  TinyVector <double,2> g;
3560  double pf = 2*M_PI*epsilon*pow(sigma,2)/area_lattice;
3561  double _V = 0.0;
3562  _V += 2*gradVz_x_64(z, sigma);
3563  double _V_old = 0.0;
3564  int g_i;
3565  int g_j;
3566  double g_magnitude;
3567  for (int k = 1; k < 2*k_max; k++) {
3568  g_i = g_i_array(k);
3569  g_j = g_j_array(k);
3570  g = (g_i * g_m) + (g_j * g_n);
3571  g_magnitude = g_magnitude_array(k);
3572  _V += gradVg_x_64(x, y, z, sigma, g_magnitude, g[0], g[1], b_1[0], b_1[1]);
3573  _V += gradVg_x_64(x, y, z, sigma, g_magnitude, g[0], g[1], b_2[0], b_2[1]);
3574 
3575  if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_1)) {
3576  _V_old = _V;
3577  flag_1 = true;
3578  }
3579  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_2)) {
3580  _V_old = _V;
3581  flag_2 = true;
3582  }
3583  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_3)) {
3584  _V_old = _V;
3585  flag_3 = true;
3586  }
3587  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && flag_1 && flag_2 && flag_3) {
3588  _V_old = _V;
3589  break;
3590  }
3591  else if ((k >= k_max) && (g_magnitude != g_magnitude_array(k+1))) {
3592  //println("WARNING: k_max reached")
3593  break;
3594  }
3595  else {
3596  _V_old = _V;
3597  }
3598  }
3599  return pf * _V;
3600 }
3601 
3602 double GrapheneLUT3DPotentialGenerate::gradV_y_64(
3603  double x, double y, double z, double sigma, double epsilon,
3604  double area_lattice, TinyVector<double,2> b_1, TinyVector<double,2> b_2,
3605  TinyVector<double,2> g_m, TinyVector<double,2> g_n, Array<int,1> g_i_array,
3606  Array<int,1> g_j_array, Array<double,1> g_magnitude_array ) {
3607  int k_max=10000;
3608  bool flag_1 = false;
3609  bool flag_2 = false;
3610  bool flag_3 = false;
3611  TinyVector <double,2> g;
3612  double pf = 2*M_PI*epsilon*pow(sigma,2)/area_lattice;
3613  double _V = 0.0;
3614  _V += 2*gradVz_y_64(z, sigma);
3615  double _V_old = 0.0;
3616  int g_i;
3617  int g_j;
3618  double g_magnitude;
3619  for (int k = 1; k < 2*k_max; k++) {
3620  g_i = g_i_array(k);
3621  g_j = g_j_array(k);
3622  g = (g_i * g_m) + (g_j * g_n);
3623  g_magnitude = g_magnitude_array(k);
3624  _V += gradVg_y_64(x, y, z, sigma, g_magnitude, g[0], g[1], b_1[0], b_1[1]);
3625  _V += gradVg_y_64(x, y, z, sigma, g_magnitude, g[0], g[1], b_2[0], b_2[1]);
3626 
3627  if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_1)) {
3628  _V_old = _V;
3629  flag_1 = true;
3630  }
3631  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_2)) {
3632  _V_old = _V;
3633  flag_2 = true;
3634  }
3635  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_3)) {
3636  _V_old = _V;
3637  flag_3 = true;
3638  }
3639  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && flag_1 && flag_2 && flag_3) {
3640  _V_old = _V;
3641  break;
3642  }
3643  else if ((k >= k_max) && (g_magnitude != g_magnitude_array(k+1))) {
3644  //println("WARNING: k_max reached")
3645  break;
3646  }
3647  else {
3648  _V_old = _V;
3649  }
3650  }
3651  return pf * _V;
3652 }
3653 
3654 double GrapheneLUT3DPotentialGenerate::gradV_z_64(
3655  double x, double y, double z, double sigma, double epsilon,
3656  double area_lattice, TinyVector<double,2> b_1, TinyVector<double,2> b_2,
3657  TinyVector<double,2> g_m, TinyVector<double,2> g_n, Array<int,1> g_i_array,
3658  Array<int,1> g_j_array, Array<double,1> g_magnitude_array ) {
3659  int k_max=10000;
3660  bool flag_1 = false;
3661  bool flag_2 = false;
3662  bool flag_3 = false;
3663  TinyVector <double,2> g;
3664  double pf = 2*M_PI*epsilon*pow(sigma,2)/area_lattice;
3665  double _V = 0.0;
3666  _V += 2*gradVz_z_64(z, sigma);
3667  double _V_old = 0.0;
3668  int g_i;
3669  int g_j;
3670  double g_magnitude;
3671  for (int k = 1; k < 2*k_max; k++) {
3672  g_i = g_i_array(k);
3673  g_j = g_j_array(k);
3674  g = (g_i * g_m) + (g_j * g_n);
3675  g_magnitude = g_magnitude_array(k);
3676  _V += gradVg_z_64(x, y, z, sigma, g_magnitude, g[0], g[1], b_1[0], b_1[1]);
3677  _V += gradVg_z_64(x, y, z, sigma, g_magnitude, g[0], g[1], b_2[0], b_2[1]);
3678 
3679  if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_1)) {
3680  _V_old = _V;
3681  flag_1 = true;
3682  }
3683  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_2)) {
3684  _V_old = _V;
3685  flag_2 = true;
3686  }
3687  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_3)) {
3688  _V_old = _V;
3689  flag_3 = true;
3690  }
3691  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && flag_1 && flag_2 && flag_3) {
3692  _V_old = _V;
3693  break;
3694  }
3695  else if ((k >= k_max) && (g_magnitude != g_magnitude_array(k+1))) {
3696  //println("WARNING: k_max reached")
3697  break;
3698  }
3699  else {
3700  _V_old = _V;
3701  }
3702  }
3703  return pf * _V;
3704 }
3705 
3706 double GrapheneLUT3DPotentialGenerate::grad2V_64(
3707  double x, double y, double z, double sigma, double epsilon,
3708  double area_lattice, TinyVector<double,2> b_1, TinyVector<double,2> b_2,
3709  TinyVector<double,2> g_m, TinyVector<double,2> g_n, Array<int,1> g_i_array,
3710  Array<int,1> g_j_array, Array<double,1> g_magnitude_array ) {
3711  int k_max=10000;
3712  bool flag_1 = false;
3713  bool flag_2 = false;
3714  bool flag_3 = false;
3715  TinyVector <double,2> g;
3716  double pf = 2*M_PI*epsilon*pow(sigma,2)/area_lattice;
3717  double _V = 0.0;
3718  _V += 2*grad2Vz_64(z, sigma);
3719  double _V_old = 0.0;
3720  int g_i;
3721  int g_j;
3722  double g_magnitude;
3723  for (int k = 1; k < 2*k_max; k++) {
3724  g_i = g_i_array(k);
3725  g_j = g_j_array(k);
3726  g = (g_i * g_m) + (g_j * g_n);
3727  g_magnitude = g_magnitude_array(k);
3728  _V += grad2Vg_64(x, y, z, sigma, g_magnitude, g[0], g[1], b_1[0], b_1[1]);
3729  _V += grad2Vg_64(x, y, z, sigma, g_magnitude, g[0], g[1], b_2[0], b_2[1]);
3730 
3731  if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_1)) {
3732  _V_old = _V;
3733  flag_1 = true;
3734  }
3735  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_2)) {
3736  _V_old = _V;
3737  flag_2 = true;
3738  }
3739  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && (!flag_3)) {
3740  _V_old = _V;
3741  flag_3 = true;
3742  }
3743  else if ((_V == _V_old) && (g_magnitude != g_magnitude_array(k+1)) && flag_1 && flag_2 && flag_3) {
3744  _V_old = _V;
3745  break;
3746  }
3747  else if ((k >= k_max) && (g_magnitude != g_magnitude_array(k+1))) {
3748  //println("WARNING: k_max reached")
3749  break;
3750  }
3751  else {
3752  _V_old = _V;
3753  }
3754  }
3755  return pf * _V;
3756 }
3757 
3758 std::tuple< TinyVector<double,2>, TinyVector<double,2>, TinyVector<double,2>,
3759  TinyVector<double,2>, TinyVector<double,2>, TinyVector<double,2>
3760  > GrapheneLUT3DPotentialGenerate::get_graphene_vectors() {
3761  return get_graphene_vectors(0.00);
3762 }
3763 
3764 std::tuple< TinyVector<double,2>, TinyVector<double,2>, TinyVector<double,2>,
3765  TinyVector<double,2>, TinyVector<double,2>, TinyVector<double,2>
3766  > GrapheneLUT3DPotentialGenerate::get_graphene_vectors( double strain ) {
3767  return get_graphene_vectors(strain, 1.42, 0.165);
3768 }
3769 
3770 std::tuple< TinyVector<double,2>, TinyVector<double,2>, TinyVector<double,2>,
3771  TinyVector<double,2>, TinyVector<double,2>, TinyVector<double,2>
3772  > GrapheneLUT3DPotentialGenerate::get_graphene_vectors(
3773  double strain, double carbon_carbon_distance, double poisson_ratio) {
3774  Array<double,2> R_strain(2,2);
3775  R_strain = -strain*poisson_ratio, 0,
3776  0, strain;
3777 
3778  TinyVector <double,2> A_m_strain0( (carbon_carbon_distance/2)*sqrt(3) , (carbon_carbon_distance/2)*3 ); //isotropic
3779  TinyVector <double,2> A_m( (carbon_carbon_distance/2)*sqrt(3)*(1 - strain*poisson_ratio) ,
3780  (carbon_carbon_distance/2)*3*(1 + strain) ); //with strain
3781 
3782  TinyVector <double,2> A_n_strain0( -(carbon_carbon_distance/2)*sqrt(3) , (carbon_carbon_distance/2)*3 ); //isotropic
3783  TinyVector <double,2> A_n( -(carbon_carbon_distance/2)*sqrt(3)*(1 - strain*poisson_ratio) ,
3784  (carbon_carbon_distance/2)*3*(1 + strain) ); //with strain
3785 
3786  //FIXME might have a problem here if A_m_60_strain0 is passed by reference
3787  TinyVector <double,2> A_m_60_strain0( carbon_carbon_distance*sqrt(3) , 0 ); // A_m_strain0 rotated 60 degrees to sit on x-axis
3788  TinyVector <double,2> A_m_60 = A_m_60_strain0; // rotated with strain
3789  A_m_60(0) += R_strain(0,0)*A_m_60_strain0(0) + R_strain(0,1)*A_m_60_strain0(1);
3790  A_m_60(1) += R_strain(1,0)*A_m_60_strain0(0) + R_strain(1,1)*A_m_60_strain0(1);
3791 
3792  TinyVector <double,2> A_n_60_strain0( (carbon_carbon_distance/2)*sqrt(3) , (carbon_carbon_distance/2)*3 ); // A_n_strain0 rotated 60
3793  TinyVector <double,2> A_n_60 = A_n_60_strain0; // rotated with strain
3794  A_n_60(0) += R_strain(0,0)*A_n_60_strain0(0) + R_strain(0,1)*A_n_60_strain0(1);
3795  A_n_60(1) += R_strain(1,0)*A_n_60_strain0(0) + R_strain(1,1)*A_n_60_strain0(1);
3796 
3797 
3798  // basis vectors
3799  TinyVector <double,2> b_1( 0, carbon_carbon_distance*(1 + strain) );
3800  TinyVector <double,2> b_2( 0, carbon_carbon_distance*2*(1 + strain) );
3801 
3802  // reciprocal lattice vectors
3803  TinyVector <double,2> g_m( 3/(1 - strain*poisson_ratio) , sqrt(3)/(1 + strain) );
3804  g_m *= (2*M_PI/3/carbon_carbon_distance);
3805  TinyVector <double,2> g_n( -g_m[0] , g_m[1] );
3806 
3807  TinyVector <double,2> g_m_60( sqrt(3)/(1 - strain*poisson_ratio), -1/(1 + strain) );
3808  g_m_60 *= (2*M_PI/3/carbon_carbon_distance);
3809  TinyVector <double,2> g_n_60( 0 , 1/(1 + strain) );
3810  g_n_60 *= (4*M_PI/3/carbon_carbon_distance);
3811 
3812  return { A_m_60, A_n_60, b_1, b_2, g_m_60, g_n_60 };
3813 }
3814 
3815 std::tuple< TinyVector<double,2>, TinyVector<double,2>, TinyVector<double,2>,
3816  TinyVector<double,2>, TinyVector<double,2>, TinyVector<double,2>
3817  > GrapheneLUT3DPotentialGenerate::get_graphene_vectors_old(
3818  double strain, double carbon_carbon_distance, double poisson_ratio) {
3819  TinyVector <double,2> A_m_old( sqrt(3)*(4 + strain - 3*strain*poisson_ratio) ,
3820  3*(4 + 3*strain - strain*poisson_ratio) ); // wrong value previously used
3821  A_m_old *= (carbon_carbon_distance/8);
3822  TinyVector <double,2> A_n_old( sqrt(3)*(-4 - strain + 3*strain*poisson_ratio),
3823  3*(4 + 3*strain - strain*poisson_ratio) ); // wrong value previously used
3824  A_n_old *= (carbon_carbon_distance/8);
3825 
3826  // basis vectors
3827  TinyVector <double,2> b_1( 0 , carbon_carbon_distance*(1 + strain) );
3828  TinyVector <double,2> b_2( 0 , carbon_carbon_distance*2*(1 + strain) );
3829 
3830  // reciprocal lattice vectors
3831  TinyVector <double,2> g_m_old( sqrt(3)*(4 + 3*strain - strain*poisson_ratio)/((4 + strain - 3*strain*poisson_ratio)*(4 + 3*strain - strain*poisson_ratio)),
3832  (4 + strain - 3*strain*poisson_ratio)/((4 + strain - 3*strain*poisson_ratio)*(4 + 3*strain - strain*poisson_ratio)) );
3833  g_m_old *= (8*M_PI/3/carbon_carbon_distance);
3834  TinyVector <double,2> g_n_old( -g_m_old[0], g_m_old[1] );
3835 
3836  return { A_m_old, A_n_old, b_1, b_2, g_m_old, g_n_old };
3837 }
3838 
3839 std::tuple< Array<int,1>, Array<int,1>, Array<double,1> >
3840 GrapheneLUT3DPotentialGenerate::get_g_magnitudes(
3841  TinyVector<double,2> g_m, TinyVector<double,2> g_n ) {
3842  int number_of_g_i = 200;
3843  int number_of_g_j = 200;
3844 
3845  int size_of_arrays = pow(number_of_g_i + number_of_g_j + 1,2);
3846  Array<double,1> g_magnitude_array(size_of_arrays);
3847  Array<int,1> g_i_array(size_of_arrays);
3848  Array<int,1> g_j_array(size_of_arrays);
3849 
3850  TinyVector<double,2> g;
3851  int k = 0;
3852  for (int g_i = -number_of_g_i; g_i < number_of_g_i + 1; g_i++) {
3853  for (int g_j = -number_of_g_j; g_j < number_of_g_j + 1; g_j++){
3854  g = (g_i * g_m) + (g_j * g_n);
3855  g_magnitude_array(k) = sqrt(dot(g,g));
3856  g_i_array(k) = g_i;
3857  g_j_array(k) = g_j;
3858  k += 1;
3859  }
3860  }
3861  int sort_indeces[size_of_arrays];
3862  for (int i = 0; i < size_of_arrays; i++) {
3863  sort_indeces[i] = i;
3864  }
3865 
3866  //std::sort( sort_indeces, sort_indeces + size_of_arrays,
3867  // [] (int i, int j) {return g_magnitude_array(i) < g_magnitude_array(j);});
3868  std::stable_sort( sort_indeces, sort_indeces + size_of_arrays,
3869  [&g_magnitude_array] (int i, int j) {return g_magnitude_array(i) < g_magnitude_array(j);});
3870  // using std::stable_sort instead of std::sort
3871  // to avoid unnecessary index re-orderings
3872  // when g_magnitude_array contains elements of equal values
3873  // FIXME may need to pass reference to g_magnitude_array in labmda function
3874  // i.e. [&g_magnitude_array](...)
3875 
3876  Array<double,1> g_magnitude_array2(size_of_arrays);
3877  Array<int,1> g_i_array2(size_of_arrays);
3878  Array<int,1> g_j_array2(size_of_arrays);
3879 
3880  for (int i = 0; i < size_of_arrays; i++) {
3881  int p = sort_indeces[i];
3882 
3883  g_magnitude_array2(i) = g_magnitude_array(p);
3884  g_i_array2(i) = g_i_array(p);
3885  g_j_array2(i) = g_j_array(p);
3886  }
3887  return std::make_tuple(g_i_array2, g_j_array2, g_magnitude_array2);
3888 }
3889 
3890 
3891 void GrapheneLUT3DPotentialGenerate::calculate_V3D_64(
3892  Array<double,3> V3D, Array<double,2> xy_x, Array<double,2> xy_y,
3893  Array<double,1> z_range, double sigma, double epsilon,
3894  double area_lattice, TinyVector<double,2> b_1, TinyVector<double,2> b_2,
3895  TinyVector<double,2> g_m, TinyVector<double,2> g_n,
3896  Array<int,1> g_i_array, Array<int,1> g_j_array,
3897  Array<double,1> g_magnitude_array ) {
3898  double x;
3899  double y;
3900  double z;
3901  for (int k = 0; k < V3D.shape()[2]; k++) {
3902  z = z_range(k);
3903  for (int j = 0; j < V3D.shape()[1]; j++) {
3904  for (int i = 0; i < V3D.shape()[0]; i++) {
3905  x = xy_x(i,j);
3906  y = xy_y(i,j);
3907  V3D(i,j,k) = V_64(x,y,z,sigma,epsilon,area_lattice,b_1,b_2,g_m,g_n,g_i_array,g_j_array,g_magnitude_array);
3908  }
3909  }
3910  }
3911 }
3912 
3913 void GrapheneLUT3DPotentialGenerate::calculate_gradV3D_x_64(
3914  Array<double,3> V3D, Array<double,2> xy_x, Array<double,2> xy_y,
3915  Array<double,1> z_range, double sigma, double epsilon,
3916  double area_lattice, TinyVector<double,2> b_1, TinyVector<double,2> b_2,
3917  TinyVector<double,2> g_m, TinyVector<double,2> g_n,
3918  Array<int,1> g_i_array, Array<int,1> g_j_array,
3919  Array<double,1> g_magnitude_array ) {
3920  double x;
3921  double y;
3922  double z;
3923  for (int k = 0; k < V3D.shape()[2]; k++) {
3924  z = z_range(k);
3925  for (int j = 0; j < V3D.shape()[1]; j++) {
3926  for (int i = 0; i < V3D.shape()[0]; i++) {
3927  x = xy_x(i,j);
3928  y = xy_y(i,j);
3929  V3D(i,j,k) = gradV_x_64(x,y,z,sigma,epsilon,area_lattice,b_1,b_2,g_m,g_n,g_i_array,g_j_array,g_magnitude_array);
3930  }
3931  }
3932  }
3933 }
3934 
3935 void GrapheneLUT3DPotentialGenerate::calculate_gradV3D_y_64(
3936  Array<double,3> V3D, Array<double,2> xy_x, Array<double,2> xy_y,
3937  Array<double,1> z_range, double sigma, double epsilon,
3938  double area_lattice, TinyVector<double,2> b_1, TinyVector<double,2> b_2,
3939  TinyVector<double,2> g_m, TinyVector<double,2> g_n,
3940  Array<int,1> g_i_array, Array<int,1> g_j_array,
3941  Array<double,1> g_magnitude_array ) {
3942  double x;
3943  double y;
3944  double z;
3945  for (int k = 0; k < V3D.shape()[2]; k++) {
3946  z = z_range(k);
3947  for (int j = 0; j < V3D.shape()[1]; j++) {
3948  for (int i = 0; i < V3D.shape()[0]; i++) {
3949  x = xy_x(i,j);
3950  y = xy_y(i,j);
3951  V3D(i,j,k) = gradV_y_64(x,y,z,sigma,epsilon,area_lattice,b_1,b_2,g_m,g_n,g_i_array,g_j_array,g_magnitude_array);
3952  }
3953  }
3954  }
3955 }
3956 
3957 void GrapheneLUT3DPotentialGenerate::calculate_gradV3D_z_64(
3958  Array<double,3> V3D, Array<double,2> xy_x, Array<double,2> xy_y,
3959  Array<double,1> z_range, double sigma, double epsilon,
3960  double area_lattice, TinyVector<double,2> b_1, TinyVector<double,2> b_2,
3961  TinyVector<double,2> g_m, TinyVector<double,2> g_n,
3962  Array<int,1> g_i_array, Array<int,1> g_j_array,
3963  Array<double,1> g_magnitude_array ) {
3964  double x;
3965  double y;
3966  double z;
3967  for (int k = 0; k < V3D.shape()[2]; k++) {
3968  z = z_range(k);
3969  for (int j = 0; j < V3D.shape()[1]; j++) {
3970  for (int i = 0; i < V3D.shape()[0]; i++) {
3971  x = xy_x(i,j);
3972  y = xy_y(i,j);
3973  V3D(i,j,k) = gradV_z_64(x,y,z,sigma,epsilon,area_lattice,b_1,b_2,g_m,g_n,g_i_array,g_j_array,g_magnitude_array);
3974  }
3975  }
3976  }
3977 }
3978 
3979 void GrapheneLUT3DPotentialGenerate::calculate_grad2V3D_64(
3980  Array<double,3> V3D, Array<double,2> xy_x, Array<double,2> xy_y,
3981  Array<double,1> z_range, double sigma, double epsilon,
3982  double area_lattice, TinyVector<double,2> b_1, TinyVector<double,2> b_2,
3983  TinyVector<double,2> g_m, TinyVector<double,2> g_n,
3984  Array<int,1> g_i_array, Array<int,1> g_j_array,
3985  Array<double,1> g_magnitude_array ) {
3986  double x;
3987  double y;
3988  double z;
3989  for (int k = 0; k < V3D.shape()[2]; k++) {
3990  z = z_range(k);
3991  for (int j = 0; j < V3D.shape()[1]; j++) {
3992  for (int i = 0; i < V3D.shape()[0]; i++) {
3993  x = xy_x(i,j);
3994  y = xy_y(i,j);
3995  V3D(i,j,k) = grad2V_64(x,y,z,sigma,epsilon,area_lattice,b_1,b_2,g_m,g_n,g_i_array,g_j_array,g_magnitude_array);
3996  }
3997  }
3998  }
3999 }
4000 
4001 std::pair<double, double> GrapheneLUT3DPotentialGenerate::get_z_min_V_min(double x, double y,
4002  double sigma, double epsilon, double area_lattice,
4003  TinyVector<double,2> b_1, TinyVector<double,2> b_2,
4004  TinyVector<double,2> g_m, TinyVector<double,2> g_n,
4005  Array<int,1> g_i_array, Array<int,1> g_j_array,
4006  Array<double,1> g_magnitude_array ) {
4007 
4008  double z_min_limit = 1.0;
4009  double z_max_limit = 5.0;
4010  const int double_bits = std::numeric_limits<double>::digits/2;
4011 
4012  std::pair<double, double> r = boost::math::tools::brent_find_minima(
4013  std::bind( &GrapheneLUT3DPotentialGenerate::V_64, this, x, y, std::placeholders::_1, sigma, epsilon,
4014  area_lattice, b_1, b_2, g_m, g_n, g_i_array, g_j_array,
4015  g_magnitude_array ),
4016  z_min_limit, z_max_limit, double_bits);
4017 
4018  return r;
4019 }
4020 
4021 std::pair<double, double> GrapheneLUT3DPotentialGenerate::get_z_V_to_find(
4022  double sigma, double epsilon, double area_lattice,
4023  TinyVector<double,2> b_1, TinyVector<double,2> b_2,
4024  TinyVector<double,2> g_m, TinyVector<double,2> g_n,
4025  Array<int,1> g_i_array, Array<int,1> g_j_array,
4026  Array<double,1> g_magnitude_array ) {
4027  double V_to_find = 10000.0;
4028  TinyVector<double,2> center_of_hexagon(0.0, 0.0);
4029  TinyVector<double,2> above_A_site(b_1[0],b_1[1]);
4030  double x = center_of_hexagon(0);
4031  double y = center_of_hexagon(1);
4032 
4033  auto [_z_min, _V_min] = get_z_min_V_min(x,y,sigma,epsilon,area_lattice,b_1,b_2,g_m,g_n,g_i_array,g_j_array,g_magnitude_array);
4034 
4035  double z_min_limit = 0.0;
4036  double z_max_limit = _z_min;
4037  const int double_bits = std::numeric_limits<double>::digits/2;
4038 
4039 
4040  auto _f_V_64 = std::bind( &GrapheneLUT3DPotentialGenerate::V_64, this, x, y, std::placeholders::_1, sigma, epsilon,
4041  area_lattice, b_1, b_2, g_m, g_n, g_i_array, g_j_array,
4042  g_magnitude_array );
4043  auto _g_lambda = [&_f_V_64](double _z, double _V) { return fabs(_f_V_64(_z) - _V); };
4044  auto _g_bind = std::bind(_g_lambda, std::placeholders::_1, V_to_find);
4045  std::pair<double, double> r = boost::math::tools::brent_find_minima(
4046  _g_bind, z_min_limit, z_max_limit, double_bits);
4047 
4048  double _z_V_to_find = r.first;
4049  double _found_V = V_64(above_A_site(0), above_A_site(1), _z_V_to_find,
4050  sigma,epsilon,area_lattice,b_1,b_2,g_m,g_n,g_i_array,g_j_array,
4051  g_magnitude_array);
4052  // return z value above hex and V value above atom
4053  return { _z_V_to_find, _found_V };
4054 }
4055 
4056 Array<double,3> GrapheneLUT3DPotentialGenerate::get_V3D(
4057  double strain, double sigma, double epsilon, int x_res, int y_res,
4058  int z_res, double z_min, double z_max ) {
4059  auto [A_m, A_n, b_1, b_2, g_m, g_n] = get_graphene_vectors(strain);
4060  auto [g_i_array, g_j_array, g_magnitude_array] = get_g_magnitudes(g_m,g_n);
4061 
4062  double cell_length_a = calculate_magnitude(A_m);
4063  double cell_length_b = calculate_magnitude(A_n);
4064  //double cell_length_c = 40.0;
4065  double cell_angle_gamma= calculate_angle(A_m,A_n);
4066  //double cell_angle_gamma_degrees = cell_angle_gamma*180/M_PI;
4067  double area_lattice = cell_length_a * cell_length_b * sin(cell_angle_gamma);
4068 
4069  Array<double,1> uc_x_range(x_res);
4070  double uc_x_min = 0.0;
4071  double uc_x_max = cell_length_a;
4072  double delta_uc_x = (uc_x_max - uc_x_min) / (x_res - 1);
4073 
4074  for (int i=0; i < x_res - 1; ++i) {
4075  uc_x_range(i) = uc_x_min + (delta_uc_x * i);
4076  }
4077  uc_x_range(x_res - 1) = uc_x_max;
4078 
4079  //double uc_dx = uc_x_range(1) - uc_x_range(0);
4080 
4081  Array<double,1> uc_y_range(y_res);
4082  double uc_y_min = 0.0;
4083  double uc_y_max = cell_length_b;
4084  double delta_uc_y = (uc_y_max - uc_y_min) / (y_res - 1);
4085 
4086  for (int i=0; i < y_res - 1; ++i) {
4087  uc_y_range(i) = uc_y_min + (delta_uc_y * i);
4088  }
4089  uc_y_range(y_res - 1) = uc_y_max;
4090 
4091  //double uc_dy = uc_y_range(1) - uc_y_range(0);
4092 
4093  Array<double,2> uc_xy_x(x_res,y_res);
4094  uc_xy_x = 0;
4095  Array<double,2> uc_xy_y(x_res,y_res);
4096  uc_xy_y = 0;
4097 
4098  for (int i=0; i < x_res; ++i) {
4099  for(int j=0; j < y_res; ++j) {
4100  uc_xy_x(i,j) = uc_x_range(i);
4101  uc_xy_y(i,j) = uc_y_range(j);
4102  }
4103  }
4104 
4105 
4106  Array<double,2> B(2,2);
4107  B = 1, cos(cell_angle_gamma),
4108  0, sin(cell_angle_gamma); //transfer to Cartesian
4109 
4110  // Set up transfer matrices to transfer from unit cell coordinates to Cartesian coordinates
4111  //A = inv(B);
4112  Array<double,2> A(2,2);
4113  //A = sin(cell_angle_gamma), -cos(cell_angle_gamma),
4114  // 0, 1;
4115  //A /= sin(cell_angle_gamma);
4116  A = 1, -1/tan(cell_angle_gamma),
4117  0, 1/sin(cell_angle_gamma); //transfer to unit cell coords
4118 
4119 
4120  Array<double,2> xy_x(x_res,y_res);
4121  Array<double,2> xy_y(x_res,y_res);
4122 
4123  for (int i=0; i < x_res; ++i) {
4124  for(int j=0; j < y_res; ++j) {
4125  //xy_x(i,j) = B(0,0)*uc_xy_x(i,j) + B(0,1)*uc_xy_y(i,j);
4126  //xy_y(i,j) = B(1,0)*uc_xy_x(i,j) + B(1,1)*uc_xy_y(i,j);
4127 
4128  xy_x(i,j) = B(0,0)*uc_x_range(i) + B(0,1)*uc_y_range(j);
4129  xy_y(i,j) = B(1,0)*uc_x_range(i) + B(1,1)*uc_y_range(j);
4130  }
4131  }
4132 
4133 
4134  Array<double,1> z_range(z_res);
4135  double delta_z = (z_max - z_min) / (z_res - 1);
4136 
4137  for (int i=0; i < z_res - 1; ++i) {
4138  z_range(i) = z_min + (delta_z * i);
4139  }
4140  z_range(z_res - 1) = z_max;
4141 
4142  //double dz = z_range(1) - z_range(0);
4143 
4144  Array<double,3> _V3D(x_res,y_res,z_res);
4145  _V3D = 0;
4146  calculate_V3D_64(_V3D,xy_x,xy_y,z_range,sigma,epsilon,area_lattice,b_1,b_2,
4147  g_m,g_n,g_i_array,g_j_array,g_magnitude_array);
4148  return _V3D;
4149 }
4150 
4151 std::pair<Array<double,3> , Array<double,1>> GrapheneLUT3DPotentialGenerate::get_V3D(
4152  double strain, double sigma, double epsilon, int x_res, int y_res,
4153  int z_res, double z_max ) {
4154  auto [A_m, A_n, b_1, b_2, g_m, g_n] = get_graphene_vectors(strain);
4155  auto [g_i_array, g_j_array, g_magnitude_array] = get_g_magnitudes(g_m,g_n);
4156 
4157  double cell_length_a = calculate_magnitude(A_m);
4158  double cell_length_b = calculate_magnitude(A_n);
4159  //double cell_length_c = 40.0;
4160  double cell_angle_gamma= calculate_angle(A_m,A_n);
4161  //double cell_angle_gamma_degrees = cell_angle_gamma*180/M_PI;
4162  double area_lattice = cell_length_a * cell_length_b * sin(cell_angle_gamma);
4163 
4164  Array<double,1> uc_x_range(x_res);
4165  double uc_x_min = 0.0;
4166  double uc_x_max = cell_length_a;
4167  double delta_uc_x = (uc_x_max - uc_x_min) / (x_res - 1);
4168 
4169  for (int i=0; i < x_res - 1; ++i) {
4170  uc_x_range(i) = uc_x_min + (delta_uc_x * i);
4171  }
4172  uc_x_range(x_res - 1) = uc_x_max;
4173 
4174  double uc_dx = uc_x_range(1) - uc_x_range(0);
4175 
4176  Array<double,1> uc_y_range(y_res);
4177  double uc_y_min = 0.0;
4178  double uc_y_max = cell_length_b;
4179  double delta_uc_y = (uc_y_max - uc_y_min) / (y_res - 1);
4180 
4181  for (int i=0; i < y_res - 1; ++i) {
4182  uc_y_range(i) = uc_y_min + (delta_uc_y * i);
4183  }
4184  uc_y_range(y_res - 1) = uc_y_max;
4185 
4186  double uc_dy = uc_y_range(1) - uc_y_range(0);
4187 
4188  Array<double,2> uc_xy_x(x_res,y_res);
4189  uc_xy_x = 0;
4190  Array<double,2> uc_xy_y(x_res,y_res);
4191  uc_xy_y = 0;
4192 
4193  for (int i=0; i < x_res; ++i) {
4194  for(int j=0; j < y_res; ++j) {
4195  uc_xy_x(i,j) = uc_x_range(i);
4196  uc_xy_y(i,j) = uc_y_range(j);
4197  }
4198  }
4199 
4200 
4201  Array<double,2> B(2,2);
4202  B = 1, cos(cell_angle_gamma),
4203  0, sin(cell_angle_gamma); //transfer to Cartesian
4204 
4205  // Set up transfer matrices to transfer from unit cell coordinates to Cartesian coordinates
4206  //A = inv(B);
4207  Array<double,2> A(2,2);
4208  //A = sin(cell_angle_gamma), -cos(cell_angle_gamma),
4209  // 0, 1;
4210  //A /= sin(cell_angle_gamma);
4211  A = 1, -1/tan(cell_angle_gamma),
4212  0, 1/sin(cell_angle_gamma); //transfer to unit cell coords
4213 
4214 
4215  Array<double,2> xy_x(x_res,y_res);
4216  Array<double,2> xy_y(x_res,y_res);
4217 
4218  for (int i=0; i < x_res; ++i) {
4219  for(int j=0; j < y_res; ++j) {
4220  //xy_x(i,j) = B(0,0)*uc_xy_x(i,j) + B(0,1)*uc_xy_y(i,j);
4221  //xy_y(i,j) = B(1,0)*uc_xy_x(i,j) + B(1,1)*uc_xy_y(i,j);
4222 
4223  xy_x(i,j) = B(0,0)*uc_x_range(i) + B(0,1)*uc_y_range(j);
4224  xy_y(i,j) = B(1,0)*uc_x_range(i) + B(1,1)*uc_y_range(j);
4225  }
4226  }
4227 
4228  auto [z_min, V_z_min] = get_z_V_to_find(sigma,epsilon,area_lattice,b_1,b_2,g_m,g_n,g_i_array,g_j_array,g_magnitude_array);
4229 
4230  Array<double,1> z_range(z_res);
4231  double delta_z = (z_max - z_min) / (z_res - 1);
4232 
4233  for (int i=0; i < z_res - 1; ++i) {
4234  z_range(i) = z_min + (delta_z * i);
4235  }
4236  z_range(z_res - 1) = z_max;
4237 
4238  double dz = z_range(1) - z_range(0);
4239 
4240  Array<double,3> _V3D(x_res,y_res,z_res);
4241  _V3D = 0;
4242  calculate_V3D_64(_V3D,xy_x,xy_y,z_range,sigma,epsilon,area_lattice,b_1,b_2,
4243  g_m,g_n,g_i_array,g_j_array,g_magnitude_array);
4244 
4245  Array<double,1> _LUTinfo(12);
4246  _LUTinfo = A(0,0), A(0,1), A(1,0), A(1,1), uc_dx, uc_dy, dz, cell_length_a,
4247  cell_length_b, z_min, z_max, V_z_min;
4248 
4249  return { _V3D, _LUTinfo };
4250 }
4251 
4252 std::tuple< Array<double,3>, Array<double,3>, Array<double,3>, Array<double,3>,
4253  Array<double,3>, Array<double,2>, Array<double,2>, Array<double,1>
4254  > GrapheneLUT3DPotentialGenerate::get_V3D_all(
4255  double strain, double sigma, double epsilon, int x_res, int y_res,
4256  int z_res, double z_max ) {
4257 
4258  auto [A_m, A_n, b_1, b_2, g_m, g_n] = get_graphene_vectors(strain);
4259  auto [g_i_array, g_j_array, g_magnitude_array] = get_g_magnitudes(g_m,g_n);
4260 
4261  double cell_length_a = calculate_magnitude(A_m);
4262  double cell_length_b = calculate_magnitude(A_n);
4263  //double cell_length_c = 40.0;
4264  double cell_angle_gamma= calculate_angle(A_m,A_n);
4265  //double cell_angle_gamma_degrees = cell_angle_gamma*180/M_PI;
4266  double area_lattice = cell_length_a * cell_length_b * sin(cell_angle_gamma);
4267 
4268  Array<double,1> uc_x_range(x_res);
4269  double uc_x_min = 0.0;
4270  double uc_x_max = cell_length_a;
4271  double delta_uc_x = (uc_x_max - uc_x_min) / (x_res - 1);
4272 
4273  for (int i=0; i < x_res - 1; ++i) {
4274  uc_x_range(i) = uc_x_min + (delta_uc_x * i);
4275  }
4276  uc_x_range(x_res - 1) = uc_x_max;
4277 
4278  double uc_dx = uc_x_range(1) - uc_x_range(0);
4279 
4280  Array<double,1> uc_y_range(y_res);
4281  double uc_y_min = 0.0;
4282  double uc_y_max = cell_length_b;
4283  double delta_uc_y = (uc_y_max - uc_y_min) / (y_res - 1);
4284 
4285  for (int i=0; i < y_res - 1; ++i) {
4286  uc_y_range(i) = uc_y_min + (delta_uc_y * i);
4287  }
4288  uc_y_range(y_res - 1) = uc_y_max;
4289 
4290  double uc_dy = uc_y_range(1) - uc_y_range(0);
4291 
4292  Array<double,2> uc_xy_x(x_res,y_res);
4293  uc_xy_x = 0;
4294  Array<double,2> uc_xy_y(x_res,y_res);
4295  uc_xy_y = 0;
4296 
4297  for (int i=0; i < x_res; ++i) {
4298  for(int j=0; j < y_res; ++j) {
4299  uc_xy_x(i,j) = uc_x_range(i);
4300  uc_xy_y(i,j) = uc_y_range(j);
4301  }
4302  }
4303 
4304 
4305  Array<double,2> B(2,2);
4306  B = 1, cos(cell_angle_gamma),
4307  0, sin(cell_angle_gamma); //transfer to Cartesian
4308 
4309  // Set up transfer matrices to transfer from unit cell coordinates to Cartesian coordinates
4310  //A = inv(B);
4311  Array<double,2> A(2,2);
4312  //A = sin(cell_angle_gamma), -cos(cell_angle_gamma),
4313  // 0, 1;
4314  //A /= sin(cell_angle_gamma);
4315  A = 1, -1/tan(cell_angle_gamma),
4316  0, 1/sin(cell_angle_gamma); //transfer to unit cell coords
4317 
4318 
4319  Array<double,2> xy_x(x_res,y_res);
4320  Array<double,2> xy_y(x_res,y_res);
4321 
4322  for (int i=0; i < x_res; ++i) {
4323  for(int j=0; j < y_res; ++j) {
4324  //xy_x(i,j) = B(0,0)*uc_xy_x(i,j) + B(0,1)*uc_xy_y(i,j);
4325  //xy_y(i,j) = B(1,0)*uc_xy_x(i,j) + B(1,1)*uc_xy_y(i,j);
4326 
4327  xy_x(i,j) = B(0,0)*uc_x_range(i) + B(0,1)*uc_y_range(j);
4328  xy_y(i,j) = B(1,0)*uc_x_range(i) + B(1,1)*uc_y_range(j);
4329  }
4330  }
4331 
4332  auto [z_min, V_z_min] = get_z_V_to_find(sigma,epsilon,area_lattice,b_1,b_2,g_m,g_n,g_i_array,g_j_array,g_magnitude_array);
4333 
4334  Array<double,1> z_range(z_res);
4335  double delta_z = (z_max - z_min) / (z_res - 1);
4336 
4337  for (int i=0; i < z_res - 1; ++i) {
4338  z_range(i) = z_min + (delta_z * i);
4339  }
4340  z_range(z_res - 1) = z_max;
4341 
4342  double dz = z_range(1) - z_range(0);
4343 
4344  Array<double,3> _V3D(x_res,y_res,z_res);
4345  Array<double,3> _gradV3D_x(x_res,y_res,z_res);
4346  Array<double,3> _gradV3D_y(x_res,y_res,z_res);
4347  Array<double,3> _gradV3D_z(x_res,y_res,z_res);
4348  Array<double,3> _grad2V3D(x_res,y_res,z_res);
4349  _V3D = 0;
4350  _gradV3D_x = 0;
4351  _gradV3D_y = 0;
4352  _gradV3D_z = 0;
4353  _grad2V3D = 0;
4354 
4355  calculate_V3D_64( _V3D, xy_x, xy_y, z_range, sigma, epsilon, area_lattice,
4356  b_1, b_2, g_m, g_n, g_i_array, g_j_array, g_magnitude_array );
4357  calculate_gradV3D_x_64( _gradV3D_x, xy_x, xy_y, z_range, sigma, epsilon,
4358  area_lattice, b_1, b_2, g_m, g_n, g_i_array, g_j_array,
4359  g_magnitude_array );
4360  calculate_gradV3D_y_64( _gradV3D_y, xy_x, xy_y, z_range, sigma, epsilon,
4361  area_lattice, b_1, b_2, g_m, g_n, g_i_array, g_j_array,
4362  g_magnitude_array );
4363  calculate_gradV3D_z_64( _gradV3D_z, xy_x, xy_y, z_range, sigma, epsilon,
4364  area_lattice, b_1, b_2, g_m, g_n, g_i_array, g_j_array,
4365  g_magnitude_array );
4366  calculate_grad2V3D_64( _grad2V3D, xy_x, xy_y, z_range, sigma, epsilon,
4367  area_lattice, b_1, b_2, g_m, g_n, g_i_array, g_j_array,
4368  g_magnitude_array );
4369  Array<double,1> _LUTinfo(12);
4370  _LUTinfo = A(0,0), A(0,1), A(1,0), A(1,1), uc_dx, uc_dy, dz, cell_length_a,
4371  cell_length_b, z_min, z_max, V_z_min;
4372  return { _V3D, _gradV3D_x, _gradV3D_y, _gradV3D_z, _grad2V3D, xy_x, xy_y,
4373  _LUTinfo };
4374 }
4375 
4376 
4377 /**************************************************************************/
4381 }
4382 
4383 // ---------------------------------------------------------------------------
4384 // ---------------------------------------------------------------------------
4385 // GrapheneLUT3DPotential Class---------------------------------------------
4386 // ---------------------------------------------------------------------------
4387 // ---------------------------------------------------------------------------
4388 
4389 /**************************************************************************/
4392 GrapheneLUT3DPotentialToBinary::GrapheneLUT3DPotentialToBinary (string graphenelut3d_file_prefix, const Container *_boxPtr) : PotentialBase() {
4393 
4394  static auto const aflags = boost::archive::no_header | boost::archive::no_tracking;
4395  /* load lookup tables */
4396  {
4397  // create and open a character archive for input
4398  std::ifstream ifs_V3d(graphenelut3d_file_prefix + "V3d.txt");
4399  // save data to archive
4400  boost::archive::text_iarchive ia_V3d(ifs_V3d,aflags);
4401  // write class instance to archive
4402  ia_V3d >> V3d;
4403  // archive and stream closed when destructors are called
4404  }
4405 
4406  {
4407  // create and open a character archive for input
4408  std::ifstream ifs_gradV3d_x(graphenelut3d_file_prefix + "gradV3d_x.txt");
4409  // save data to archive
4410  boost::archive::text_iarchive ia_gradV3d_x(ifs_gradV3d_x,aflags);
4411  // write class instance to archive
4412  ia_gradV3d_x >> gradV3d_x;
4413  // archive and stream closed when destructors are called
4414  }
4415 
4416  {
4417  // create and open a character archive for input
4418  std::ifstream ifs_gradV3d_y(graphenelut3d_file_prefix + "gradV3d_y.txt");
4419  // save data to archive
4420  boost::archive::text_iarchive ia_gradV3d_y(ifs_gradV3d_y,aflags);
4421  // write class instance to archive
4422  ia_gradV3d_y >> gradV3d_y;
4423  // archive and stream closed when destructors are called
4424  }
4425 
4426  {
4427  // create and open a character archive for input
4428  std::ifstream ifs_gradV3d_z(graphenelut3d_file_prefix + "gradV3d_z.txt");
4429  // save data to archive
4430  boost::archive::text_iarchive ia_gradV3d_z(ifs_gradV3d_z,aflags);
4431  // write class instance to archive
4432  ia_gradV3d_z >> gradV3d_z;
4433  // archive and stream closed when destructors are called
4434  }
4435 
4436  {
4437  // create and open a character archive for input
4438  std::ifstream ifs_grad2V3d(graphenelut3d_file_prefix + "grad2V3d.txt");
4439  // save data to archive
4440  boost::archive::text_iarchive ia_grad2V3d(ifs_grad2V3d,aflags);
4441  // write class instance to archive
4442  ia_grad2V3d >> grad2V3d;
4443  // archive and stream closed when destructors are called
4444  }
4445 
4446  {
4447  // create and open a character archive for input
4448  std::ifstream ifs_LUTinfo(graphenelut3d_file_prefix + "LUTinfo.txt");
4449  // save data to archive
4450  boost::archive::text_iarchive ia_LUTinfo(ifs_LUTinfo,aflags);
4451  // write class instance to archive
4452  ia_LUTinfo >> LUTinfo;
4453  // archive and stream closed when destructors are called
4454  }
4455 
4456  // create and open a character archive for output
4457  std::ofstream ofs(graphenelut3d_file_prefix + "serialized.dat");
4458  {
4459  // save data to archive
4460  boost::archive::binary_oarchive oa(ofs,aflags);
4461  // write class instance to archive
4462  oa << V3d << gradV3d_x << gradV3d_y << gradV3d_z << grad2V3d << LUTinfo;
4463  // archive and stream closed when destructors are called
4464  }
4465 
4466  std::cout << "Finished converting binary file " <<
4467  graphenelut3d_file_prefix + "serialized.txt" << " to binary file " <<
4468  graphenelut3d_file_prefix + "serialized.dat" << ", exiting." <<
4469  std::endl;
4470 }
4471 
4472 /**************************************************************************/
4476  V3d.free();
4477  gradV3d_x.free();
4478  gradV3d_y.free();
4479  gradV3d_z.free();
4480  grad2V3d.free();
4481  LUTinfo.free();
4482 }
4483 
4484 // ---------------------------------------------------------------------------
4485 // ---------------------------------------------------------------------------
4486 // GrapheneLUT3DPotential Class---------------------------------------------
4487 // ---------------------------------------------------------------------------
4488 // ---------------------------------------------------------------------------
4489 
4490 /**************************************************************************/
4493 GrapheneLUT3DPotentialToText::GrapheneLUT3DPotentialToText (string graphenelut3d_file_prefix, const Container *_boxPtr) : PotentialBase() {
4494 
4495  static auto const aflags = boost::archive::no_header | boost::archive::no_tracking;
4496  /* load lookup tables */
4497  {
4498  // create and open an archive for input
4499  std::ifstream ifs(graphenelut3d_file_prefix + "serialized.dat");
4500  boost::archive::binary_iarchive ia(ifs,aflags);
4501  // read class state from archive
4502  ia >> V3d >> gradV3d_x >> gradV3d_y >> gradV3d_z >> grad2V3d >> LUTinfo;
4503  // archive and stream closed when destructors are called
4504  }
4505 
4506  // create and open a character archive for output
4507  std::ofstream ofs_V3d(graphenelut3d_file_prefix + "V3d.txt");
4508  {
4509  // save data to archive
4510  boost::archive::text_oarchive oa_V3d(ofs_V3d,aflags);
4511  // write class instance to archive
4512  oa_V3d << V3d;
4513  // archive and stream closed when destructors are called
4514  }
4515 
4516  // create and open a character archive for output
4517  std::ofstream ofs_gradV3d_x(graphenelut3d_file_prefix + "gradV3d_x.txt");
4518  {
4519  // save data to archive
4520  boost::archive::text_oarchive oa_gradV3d_x(ofs_gradV3d_x,aflags);
4521  // write class instance to archive
4522  oa_gradV3d_x << gradV3d_x;
4523  // archive and stream closed when destructors are called
4524  }
4525 
4526  // create and open a character archive for output
4527  std::ofstream ofs_gradV3d_y(graphenelut3d_file_prefix + "gradV3d_y.txt");
4528  {
4529  // save data to archive
4530  boost::archive::text_oarchive oa_gradV3d_y(ofs_gradV3d_y,aflags);
4531  // write class instance to archive
4532  oa_gradV3d_y << gradV3d_y;
4533  // archive and stream closed when destructors are called
4534  }
4535 
4536  // create and open a character archive for output
4537  std::ofstream ofs_gradV3d_z(graphenelut3d_file_prefix + "gradV3d_z.txt");
4538  {
4539  // save data to archive
4540  boost::archive::text_oarchive oa_gradV3d_z(ofs_gradV3d_z,aflags);
4541  // write class instance to archive
4542  oa_gradV3d_z << gradV3d_z;
4543  // archive and stream closed when destructors are called
4544  }
4545 
4546  // create and open a character archive for output
4547  std::ofstream ofs_grad2V3d(graphenelut3d_file_prefix + "grad2V3d.txt");
4548  {
4549  // save data to archive
4550  boost::archive::text_oarchive oa_grad2V3d(ofs_grad2V3d,aflags);
4551  // write class instance to archive
4552  oa_grad2V3d << grad2V3d;
4553  // archive and stream closed when destructors are called
4554  }
4555 
4556  // create and open a character archive for output
4557  std::ofstream ofs_LUTinfo(graphenelut3d_file_prefix + "LUTinfo.txt");
4558  {
4559  // save data to archive
4560  boost::archive::text_oarchive oa_LUTinfo(ofs_LUTinfo,aflags);
4561  // write class instance to archive
4562  oa_LUTinfo << LUTinfo;
4563  // archive and stream closed when destructors are called
4564  }
4565 
4566  std::cout << "Finished converting binary file " <<
4567  graphenelut3d_file_prefix + "serialized.dat" << " to text files " <<
4568  graphenelut3d_file_prefix + "<V3d|gradV3d_x|gradV3d_y|gradV3d_z|grad2V3d|LUTinfo>.txt" << ", exiting." <<
4569  std::endl;
4570 }
4571 
4572 /**************************************************************************/
4576  V3d.free();
4577  gradV3d_x.free();
4578  gradV3d_y.free();
4579  gradV3d_z.free();
4580  grad2V3d.free();
4581  LUTinfo.free();
4582 }
~AzizPotential()
Destructor.
Definition: potential.cpp:1515
dVec gradV(const dVec &)
Return the gradient of aziz potential for separation r using a lookup table.
Definition: potential.h:743
AzizPotential(const Container *)
Constructor.
Definition: potential.cpp:1475
double V(const dVec &)
Return the aziz potential for separation r using a lookup table.
Definition: potential.h:731
File * file(string type)
Get method returning file object.
Definition: communicator.h:85
double m() const
Get mass.
Definition: constants.h:45
double lambda() const
Get lambda = hbar^2/(2mk_B)
Definition: constants.h:46
double rc2() const
Get potential cutoff squared.
Definition: constants.h:48
double L() const
Get maximum side length.
Definition: constants.h:52
double tau() const
Get imaginary time step.
Definition: constants.h:44
The base class which holds details on the generalized box that our system will be simulated inside of...
Definition: container.h:24
double volume
The volume of the container in A^3.
Definition: container.h:35
void putInBC(dVec &r) const
Place a vector in boundary conditions.
Definition: container.h:50
virtual void putInside(dVec &) const =0
Place a vector inside the simulation cell.
dVec side
The linear dimensions of the box.
Definition: container.h:31
double maxSep
The maximum possible separation for 2 beads on the same timeslice.
Definition: container.h:37
~Delta1DPotential()
Destructor.
Definition: potential.cpp:2187
double dVdlambda(const dVec &, const dVec &)
The derivative of the effective potential with respect to lambda.
Definition: potential.cpp:2319
double dVdtau(const dVec &, const dVec &)
The derivative of the effective potential with respect to tau.
Definition: potential.cpp:2344
virtual double V(const dVec &r)
The classical potential.
Definition: potential.h:1066
Delta1DPotential(double)
Constructor.
Definition: potential.cpp:2170
DeltaPotential(double, double)
Constructor.
Definition: potential.cpp:387
~DeltaPotential()
Destructor.
Definition: potential.cpp:397
DipolePotential()
Constructor.
Definition: potential.cpp:462
~DipolePotential()
Destructor.
Definition: potential.cpp:473
dVec gradV(const dVec &r)
The gradient of the total potential coming from the interaction of a particle with all fixed particle...
Definition: potential.cpp:592
double V(const dVec &r)
The total potential coming from the interaction of a particle with all fixed particles.
Definition: potential.cpp:568
~FixedAzizPotential()
Destructor.
Definition: potential.cpp:557
FixedAzizPotential(const Container *)
Constructor.
Definition: potential.cpp:490
Array< dVec, 1 > initialConfig(const Container *, MTRand &, const int)
Initial configuration corresponding to FixedAziz potential.
Definition: potential.cpp:620
FixedPositionLJPotential(const double, const double, const Container *)
Constructor.
Definition: potential.cpp:676
~FixedPositionLJPotential()
Destructor.
Definition: potential.cpp:731
double V(const dVec &r)
Return the value of the van der Waals' interaction between a graphene sheet and a helium adatom at a ...
Definition: potential.cpp:742
FreePotential()
Constructor.
Definition: potential.cpp:277
~FreePotential()
Destructor.
Definition: potential.cpp:283
Array< double, 1 > getExcLen()
get the exclusion lengths ay and az
Definition: potential.cpp:2033
Gasparini_1_Potential(double, double, const Container *)
Constructor.
Definition: potential.cpp:1853
Array< dVec, 1 > initialConfig(const Container *, MTRand &, const int)
Initial configuration corresponding to FixedAziz potential.
Definition: potential.cpp:1874
~Gasparini_1_Potential()
Destructor.
Definition: potential.cpp:1865
~GrapheneLUT3DPotentialGenerate()
Destructor.
Definition: potential.cpp:4380
GrapheneLUT3DPotentialGenerate(const double, const double, const double, const double, const double, const Container *)
Constructor.
Definition: potential.cpp:3302
~GrapheneLUT3DPotentialToBinary()
Destructor.
Definition: potential.cpp:4475
GrapheneLUT3DPotentialToBinary(const string, const Container *)
Constructor.
Definition: potential.cpp:4392
GrapheneLUT3DPotentialToText(const string, const Container *)
Constructor.
Definition: potential.cpp:4493
~GrapheneLUT3DPotentialToText()
Destructor.
Definition: potential.cpp:4575
Array< dVec, 1 > initialConfig(const Container *, MTRand &, const int)
Initial configuration corresponding to graphene-helium vdW potential.
Definition: potential.cpp:3222
dVec gradV(const dVec &)
Return the gradient of the van der Waals' interaction between a graphene sheet and a helium adatom at...
Definition: potential.cpp:3091
double V(const dVec &)
Return the value of the van der Waals' interaction between a graphene sheet and a helium adatom at a ...
Definition: potential.cpp:3062
double grad2V(const dVec &)
Return the gradient of the van der Waals' interaction between a graphene sheet and a helium adatom at...
Definition: potential.cpp:3120
~GrapheneLUT3DPotential()
Destructor.
Definition: potential.cpp:3046
GrapheneLUT3DPotential(const string, const Container *)
Constructor.
Definition: potential.cpp:3006
Array< dVec, 1 > initialConfig1(const Container *, MTRand &, const int)
Return an initial particle configuration.
Definition: potential.cpp:3140
dVec gradV(const dVec &r)
Return the gradient of the van der Waals' interaction between a graphene sheet and a helium adatom at...
Definition: potential.cpp:2902
Array< dVec, 1 > initialConfig(const Container *, MTRand &, const int)
Initial configuration corresponding to graphene-helium vdW potential.
Definition: potential.cpp:2948
GrapheneLUTPotential(const double, const double, const double, const double, const double, const Container *)
Constructor.
Definition: potential.cpp:2658
double V(const dVec &r)
Return the value of the van der Waals' interaction between a graphene sheet and a helium adatom at a ...
Definition: potential.cpp:2853
~GrapheneLUTPotential()
Destructor.
Definition: potential.cpp:2840
GraphenePotential(const double, const double, const double, const double, const double)
Constructor.
Definition: potential.cpp:2475
double V(const dVec &r)
Return the value of the van der Waals' interaction between a graphene sheet and a helium adatom at a ...
Definition: potential.cpp:2526
~GraphenePotential()
Destructor.
Definition: potential.cpp:2515
Array< dVec, 1 > initialConfig(const Container *, MTRand &, const int)
Initial configuration corresponding to graphene-helium vdW potential.
Definition: potential.cpp:2571
~HardCylinderPotential()
Destructor.
Definition: potential.cpp:789
HardCylinderPotential(const double)
Constructor.
Definition: potential.cpp:781
HardRodPotential(double)
Constructor.
Definition: potential.cpp:2367
double dVdlambda(const dVec &, const dVec &)
The derivative of the effective potential with respect to lambda.
Definition: potential.cpp:2426
~HardRodPotential()
Destructor.
Definition: potential.cpp:2376
virtual double V(const dVec &r)
The classical potential.
Definition: potential.h:1108
double dVdtau(const dVec &, const dVec &)
The derivative of the effective potential with respect to tau.
Definition: potential.cpp:2451
virtual double V(const dVec &r)
The classical potential.
Definition: potential.h:1038
double dVdlambda(const dVec &, const dVec &)
The derivative of the effective potential with respect to lambda.
Definition: potential.cpp:2108
~HardSpherePotential()
Destructor.
Definition: potential.cpp:2062
HardSpherePotential(double)
Constructor.
Definition: potential.cpp:2053
double dVdtau(const dVec &, const dVec &)
The derivative of the effective potential with respect to tau.
Definition: potential.cpp:2140
HarmonicCylinderPotential(const double)
Constructor.
Definition: potential.cpp:358
~HarmonicCylinderPotential()
Destructor.
Definition: potential.cpp:371
Array< dVec, 1 > initialConfig(const Container *, MTRand &, const int)
Initial configuration corresponding to Harmonic potential.
Definition: potential.cpp:332
~HarmonicPotential()
Destructor.
Definition: potential.cpp:324
LJCylinderPotential(const double)
Constructor.
Definition: potential.cpp:1027
~LJCylinderPotential()
Destructor.
Definition: potential.cpp:1077
Array< dVec, 1 > initialConfig(const Container *, MTRand &, const int)
Initial configuration corresponding to the LJ cylinder potential.
Definition: potential.cpp:1191
LJHourGlassPotential(const double, const double, const double)
Constructor.
Definition: potential.cpp:1273
Array< dVec, 1 > initialConfig(const Container *, MTRand &, const int)
Initial configuration corresponding to the LJ cylinder potential.
Definition: potential.cpp:1429
double V(const dVec &)
The integrated LJ hour glass potential.
Definition: potential.cpp:1314
Array< dVec, 1 > initialConfig1(const Container *, MTRand &, const int)
Return a set of initial positions inside the hourglass.
Definition: potential.cpp:1354
~LJHourGlassPotential()
Destructor.
Definition: potential.cpp:1301
The particle (bead) lookup table.
Definition: lookuptable.h:29
int fullNumBeads
The full number of active beads in beadList;.
Definition: lookuptable.h:44
int getTotNumGridBoxes()
Return the total number of grid boxes.
Definition: lookuptable.h:49
void updateGrid(const Path &)
We update the full nearest neighbor grid by filling up all data arrays at all time slices.
void updateFullInteractionList(const beadLocator &, const int)
Fill up the fullBeadList array with a list of beads in the same grid box as the supplied beadIndex an...
Array< beadLocator, 1 > fullBeadList
The full dynamic list of interacting beads.
Definition: lookuptable.h:41
const Container * boxPtr
The simulation cell.
Definition: lookuptable.h:35
int gridNumber(const dVec &)
Given a particle position, return the grid number for the nearest neighbor lookup table.
Definition: lookuptable.h:154
~LorentzianPotential()
Destructor.
Definition: potential.cpp:424
LorentzianPotential(double, double)
Constructor.
Definition: potential.cpp:413
PlatedLJCylinderPotential(const double, const double, const double, const double, const double)
Constructor.
Definition: potential.cpp:800
~PlatedLJCylinderPotential()
Destructor.
Definition: potential.cpp:846
Array< dVec, 1 > initialConfig(const Container *, MTRand &, const int)
Initial configuration corresponding to the LJ cylinder potential.
Definition: potential.cpp:948
The base class from which all specific potentials are derived from.
Definition: potential.h:32
double deltaSeparation(double sep1, double sep2) const
Return the minimum image difference for 1D separations.
Definition: potential.cpp:119
PotentialBase()
Constructor.
Definition: potential.cpp:25
virtual Array< dVec, 1 > initialConfig(const Container *, MTRand &, const int)
Default Initial configuration of particles.
Definition: potential.cpp:41
virtual ~PotentialBase()
Destructor.
Definition: potential.cpp:32
double tailV
Tail correction factor.
Definition: potential.h:61
void output(const double)
A debug method that output's the potential to a supplied separation.
Definition: potential.cpp:106
virtual Array< double, 1 > getExcLen()
Array to hold data elements.
Definition: potential.cpp:135
virtual double V(const dVec &)
The potential.
Definition: potential.h:39
~SingleWellPotential()
Destructor.
Definition: potential.cpp:301
SingleWellPotential()
Constructor.
Definition: potential.cpp:295
~SutherlandPotential()
Destructor.
Definition: potential.cpp:448
SutherlandPotential(double)
Constructor.
Definition: potential.cpp:439
~SzalewiczPotential()
Destructor.
Definition: potential.cpp:1694
SzalewiczPotential(const Container *)
Constructor.
Definition: potential.cpp:1625
Pre-tabulated potential for complicated functions.
Definition: potential.h:80
Array< double, 1 > lookupd2Vdr2
A lookup table for d2Vint/dr2.
Definition: potential.h:88
TabulatedPotential()
Constructor.
Definition: potential.cpp:149
int tableLength
The number of elements in the lookup table.
Definition: potential.h:91
TinyVector< double, 2 > extdVdr
Extremal value of dV/dr.
Definition: potential.h:94
virtual ~TabulatedPotential()
Destructor.
Definition: potential.cpp:158
virtual double valuedVdr(const double)=0
The functional value of dV/dr.
TinyVector< double, 2 > extV
Extremal value of V.
Definition: potential.h:93
void initLookupTable(const double, const double)
Given a discretization factor and the system size, create and fill the lookup tables for the potentia...
Definition: potential.cpp:168
Array< double, 1 > lookupdVdr
A lookup table for dVint/dr.
Definition: potential.h:87
virtual double valued2Vdr2(const double)=0
The functional value of d2V/dr2.
virtual double valueV(const double)=0
The functional value of V.
virtual double direct(const Array< double, 1 > &, const TinyVector< double, 2 > &, const double)
Use a direct lookup for the potential table.
Definition: potential.cpp:255
virtual double newtonGregory(const Array< double, 1 > &, const TinyVector< double, 2 > &, const double)
Use the Newton-Gregory forward difference method to do a 2-point lookup on the potential table.
Definition: potential.cpp:226
TinyVector< double, 2 > extd2Vdr2
Extremal value of d2V/dr2.
Definition: potential.h:95
Array< double, 1 > lookupV
A potential lookup table.
Definition: potential.h:86
double dr
The discretization for the lookup table.
Definition: potential.h:90
#define NDIM
Number of spatial dimnsions.
Definition: common.h:71
int ipow(int base, int power)
Return the integer value of a number raised to a power.
Definition: common.h:136
#define EPS
A small number.
Definition: common.h:94
#define LBIG
The log of a big number.
Definition: common.h:97
TinyVector< double, NDIM > dVec
A NDIM-vector of type double.
Definition: common.h:111
#define PIMC_ASSERT(X)
Rename assert method.
Definition: common.h:64
TinyVector< int, NDIM > iVec
A NDIM-vector of type integer.
Definition: common.h:114
#define XXX
Used to refer to a nonsense beadIndex.
Definition: common.h:98
Class definitions for all file input/output.
Communicator * communicate()
Global public access to the communcator singleton.
Definition: communicator.h:121
ConstantParameters * constants()
Global public access to the constants.
Definition: constants.h:201
LookupTable class definition.
Path class definition.
All possible potential classes.