Path Integral Quantum Monte Carlo
Public Member Functions | Data Fields
Gasparini_1_Potential Class Reference

Computes potential energy for Gasparini potential. More...

#include <potential.h>

+ Inheritance diagram for Gasparini_1_Potential:
+ Collaboration diagram for Gasparini_1_Potential:

Public Member Functions

 Gasparini_1_Potential (double, double, const Container *)
 Constructor.
 
 ~Gasparini_1_Potential ()
 Destructor.
 
double V (const dVec &r)
 The potential.
 
dVec gradV (const dVec &)
 The gradient of the potential.
 
double grad2V (const dVec &r)
 Laplacian of the potential.
 
Array< dVec, 1 > initialConfig (const Container *, MTRand &, const int)
 Initial configuration corresponding to FixedAziz potential. More...
 
Array< double, 1 > getExcLen ()
 get the exclusion lengths ay and az More...
 
- Public Member Functions inherited from PotentialBase
 PotentialBase ()
 Constructor.
 
virtual ~PotentialBase ()
 Destructor.
 
virtual double V (const dVec &, const dVec &)
 The effective potential for the pair product approximation.
 
virtual double dVdlambda (const dVec &, const dVec &)
 The derivative of the effective potential with respect to lambda and tau.
 
virtual double dVdtau (const dVec &, const dVec &)
 
void output (const double)
 A debug method that output's the potential to a supplied separation. More...
 

Data Fields

const double excZ
 
const double excY
 
const double V0
 
- Data Fields inherited from PotentialBase
double tailV
 Tail correction factor.
 

Additional Inherited Members

- Protected Member Functions inherited from PotentialBase
double deltaSeparation (double sep1, double sep2) const
 Return the minimum image difference for 1D separations.
 

Detailed Description

Computes potential energy for Gasparini potential.

This is the case of gradV = 0 for all positions in cell.

Definition at line 990 of file potential.h.

Member Function Documentation

◆ getExcLen()

Array< double, 1 > Gasparini_1_Potential::getExcLen ( )
virtual

get the exclusion lengths ay and az

Returns the exclusion lengths ay and az.

Reimplemented from PotentialBase.

Definition at line 2033 of file potential.cpp.

2033  {
2034 
2035  Array<double, 1> excLens(2);
2036  excLens(0) = excY; // was ay
2037  excLens(1) = excZ; // was az
2038 
2039  return (excLens);
2040 }

◆ initialConfig()

Array< dVec, 1 > Gasparini_1_Potential::initialConfig ( const Container boxPtr,
MTRand &  random,
const int  numParticles 
)
virtual

Initial configuration corresponding to FixedAziz potential.

Return initial particle configuration.

Return initial positions to exclude volume in the simulation cell.

Reimplemented from PotentialBase.

Definition at line 1874 of file potential.cpp.

1875  {
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 }
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
#define NDIM
Number of spatial dimnsions.
Definition: common.h:71
#define EPS
A small number.
Definition: common.h:94
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
+ Here is the call graph for this function:

The documentation for this class was generated from the following files: