Path Integral Quantum Monte Carlo
wavefunction.cpp
Go to the documentation of this file.
1 
9 #include "wavefunction.h"
10 #include "path.h"
11 
12 // ---------------------------------------------------------------------------
13 // ---------------------------------------------------------------------------
14 // WAVEFUNCTION BASE CLASS ---------------------------------------------------
15 // ---------------------------------------------------------------------------
16 // ---------------------------------------------------------------------------
17 
18 /**************************************************************************/
22  string _name) :
23  name(_name),
24  path(_path),
25  lookup(_lookup)
26 {
27  // empty constructor
28 }
29 
30 /**************************************************************************/
34  // empty destructor
35 }
36 
37 // ---------------------------------------------------------------------------
38 // ---------------------------------------------------------------------------
39 // SECH WAVEFUNCTION CLASS ---------------------------------------------------
40 // ---------------------------------------------------------------------------
41 // ---------------------------------------------------------------------------
42 
43 /**************************************************************************/
46 SechWaveFunction::SechWaveFunction(const Path &_path,LookupTable &_lookup, string _name) :
47  WaveFunctionBase(_path,_lookup,_name)
48 {
49  /* Set the parameter to its optimized value */
50  a = sqrt(0.5*M_PI);
51 }
52 
53 /**************************************************************************/
57  // empty destructor
58 }
59 
60 /**************************************************************************/
63 double SechWaveFunction::PsiTrial(const int slice) {
64 
65  /* The cumulative value */
66  double psiT = 1.0;
67  beadLocator beadIndex;
68  dVec pos;
69  double r;
70 
71  for (int ptcl = 0; ptcl < path.numBeadsAtSlice(slice); ptcl++) {
72  beadIndex = slice,ptcl;
73  pos = path(beadIndex);
74  r = sqrt(dot(pos,pos));
75  psiT *= 1.0/cosh(a*r);
76  }
77 
78  return psiT;
79 }
80 
81 
82 // ---------------------------------------------------------------------------
83 // ---------------------------------------------------------------------------
84 // JASTROW WAVEFUNCTION CLASS ---------------------------------------------------
85 // ---------------------------------------------------------------------------
86 // ---------------------------------------------------------------------------
87 
88 /**************************************************************************/
91 JastrowWaveFunction::JastrowWaveFunction(const Path &_path,LookupTable &_lookup, string _name) :
92 WaveFunctionBase(_path,_lookup,_name)
93 {
94  /* Set the parameter to its optimized value */
95  alpha = 19.0;
96  beta = 0.12;
97  //beta = 3.07;
98 }
99 
100 /**************************************************************************/
104  // empty destructor
105 }
106 
107 /**************************************************************************/
110 double JastrowWaveFunction::PsiTrial(const double r) {
111 
112  double psiT = exp( -0.5*alpha/(1.0+beta*pow(r,5.0)));
113  return psiT;
114 }
115 
116 /**************************************************************************/
119 double JastrowWaveFunction::delPsiTrial(const double r) {
120 
121  double delpsiT = 2.5*alpha*beta*pow(r,4.0)/pow((1.0+beta*pow(r,5.0)),2.0);
122  return delpsiT;
123 }
124 
125 /**************************************************************************/
128 double JastrowWaveFunction::delSqPsiTrial(const double r) {
129 
130  double delSqPsiT = -1.25*alpha*beta*pow(r,3.0)*(-8.0+
131  beta*(4.0-5.0*alpha)*pow(r,5.0)+12.0*pow(beta,2.0)*pow(r,10.0) )/
132  ( pow((1.0+beta*pow(r,5.0)),4.0) );
133  return delSqPsiT;
134 }
135 
136 /**************************************************************************/
139 double JastrowWaveFunction::PsiTrial(const int slice) {
140 
141  /* The cumulative value */
142  double psiT = 1.0;
143  int numParticles = path.numBeadsAtSlice(slice);
144  dVec sep; // The spatial separation between beads.
145  double r; // Distance between beads
146  beadLocator bead1,bead2;
147  bead1[0] = bead2[0] = slice;
148 
149  for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
150  /* The loop over all other particles, to find the total interaction
151  * potential */
152  for (bead2[1] = bead1[1]+1; bead2[1] < numParticles; bead2[1]++) {
153  sep = path.getSeparation(bead2,bead1);
154  r = sqrt(dot(sep,sep));
155  psiT *= PsiTrial(r);
156  } // bead2
157 
158  } // bead1
159 
160  return psiT;
161 }
162 
163 /**************************************************************************/
166 double JastrowWaveFunction::gradSqPsiTrial(const int slice) {
167 
168  /* The cumulative value */
169  double gradSqPsiT = 1.0;
170  double delPsi12;
171  int numParticles = path.numBeadsAtSlice(slice);
172  dVec sep; // The spatial separation between beads.
173  double r; // Distance between beads
174  beadLocator bead1,bead2,bead3;
175  bead1[0] = bead2[0] = bead3[0] = slice;
176 
177  for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
178  /* The loop over all other particles, to find the total interaction
179  * potential */
180  for (bead2[1] = bead1[1]+1; bead2[1] < numParticles; bead2[1]++) {
181  sep = path.getSeparation(bead2,bead1);
182  r = sqrt(dot(sep,sep));
183  gradSqPsiT += 2.0*delSqPsiTrial(r);
184  delPsi12 = delPsiTrial(r);
185  for (bead3[1] = bead2[1]+1; bead3[1] < numParticles; bead3[1]++) {
186  sep = path.getSeparation(bead3,bead1);
187  r = sqrt(dot(sep,sep));
188  gradSqPsiT += 2.0*delPsi12*delPsiTrial(r);
189  }// bead 3
190  } // bead2
191 
192  } // bead1
193 
194  return gradSqPsiT;
195 }
196 
197 
198 // ---------------------------------------------------------------------------
199 // ---------------------------------------------------------------------------
200 // LiebLiniger WAVEFUNCTION CLASS --------------------------------------------
201 // ---------------------------------------------------------------------------
202 // ---------------------------------------------------------------------------
203 
204 /**************************************************************************/
208 WaveFunctionBase(_path,_lookup,_name)
209 {
210  /* Set the parameter to its optimized value */
211  R = constants()->R_LL_wfn();
212  k = constants()->k_LL_wfn();
213 }
214 
215 /**************************************************************************/
219  // empty destructor
220 }
221 
222 /**************************************************************************/
225 double LiebLinigerWaveFunction::PsiTrial(const double r) {
226 
227  double psiT = 1.0;
228  if(r<R)
229  psiT = cos(k*(abs(r)-R));
230  return psiT;
231 }
232 
233 /**************************************************************************/
236 double LiebLinigerWaveFunction::delPsiTrial(const double r) {
237 
238  return 0.0;
239 }
240 
241 /**************************************************************************/
245 
246  return 0.0;
247 }
248 
249 /**************************************************************************/
252 double LiebLinigerWaveFunction::PsiTrial(const int slice) {
253 
254  /* The cumulative value */
255  double psiT = 1.0;
256  Array <bool,1> doParticles(path.numBeadsAtSlice(slice));
257  doParticles = true;
258 
259  dVec sep; // The spatial separation between beads.
260  double r; // Distance between beads
261  beadLocator bead1,bead2;
262  bead1[0] = bead2[0] = slice;
263 
264  /* No cutoff */
265 // for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
266 // /* The loop over all other particles, to find the total interaction
267 // * potential */
268 // for (bead2[1] = bead1[1]+1; bead2[1] < numParticles; bead2[1]++) {
269 // sep = path.getSeparation(bead2,bead1);
270 // r = sqrt(dot(sep,sep));
271 // psiT *= PsiTrial(r);
272 // } // bead2
273 //
274 // } // bead1
275 
276  /* Using cutoff */
277  for (bead1[1] = 0; bead1[1] < path.numBeadsAtSlice(slice); bead1[1]++) {
278 
279  doParticles(bead1[1]) = false;
280 
281  /* Get the interaction list */
283 
284  /* Sum the interaction potential over all NN beads */
285  for (int n = 0; n < lookup.numBeads; n++) {
286  bead2 = lookup.beadList(n);
287  if (doParticles(bead2[1])) {
288  sep = path.getSeparation(bead2,bead1);
289  r = sqrt(dot(sep,sep));
290  psiT *= PsiTrial(r);
291  }
292  } // n
293 
294  } // bead1
295 
296  return psiT;
297 }
298 
299 /**************************************************************************/
303 
304  /* The cumulative value */
305  double psiT = 1.0;
306 
307  dVec sep; // The spatial separation between beads.
308  double r; // Distance between beads
309 
310  /* No cutoff */
311  // for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
312  // /* The loop over all other particles, to find the total interaction
313  // * potential */
314  // for (bead2[1] = bead1[1]+1; bead2[1] < numParticles; bead2[1]++) {
315  // sep = path.getSeparation(bead2,bead1);
316  // r = sqrt(dot(sep,sep));
317  // psiT *= PsiTrial(r);
318  // } // bead2
319  //
320  // } // bead1
321 
322  /* We only continue if bead1 is turned on */
323  if (path.worm.beadOn(bead1)) {
324 
325  /* Fill up th nearest neighbor list */
327 
328  /* Sum the interaction potential over all NN beads */
329  for (int n = 0; n < lookup.numBeads; n++) {
330 
331  sep = path.getSeparation(bead1,lookup.beadList(n));
332  r = sqrt(dot(sep,sep));
333  psiT *= PsiTrial(r);
334  }
335  }
336 
337  return psiT;
338 }
339 
340 
341 /**************************************************************************/
345 
346  return 0.0;
347 }
348 
349 
350 // ---------------------------------------------------------------------------
351 // ---------------------------------------------------------------------------
352 // SUTHERLAND WAVEFUNCTION CLASS ---------------------------------------------
353 // ---------------------------------------------------------------------------
354 // ---------------------------------------------------------------------------
355 
356 /**************************************************************************/
360  double _lambda, string _name) :
361 WaveFunctionBase(_path,_lookup,_name)
362 {
363  // The Sutherland model value of the interaction paramter \lambda
364  lambda = _lambda;
365 
366  // pi/L
367  pioL = M_PI/constants()->L();
368 }
369 
370 /**************************************************************************/
374  // empty destructor
375 }
376 
377 /**************************************************************************/
380 double SutherlandWaveFunction::PsiTrial(const int slice) {
381 
382  /* The cumulative value */
383  double psiT = 1.0;
384 
385  /* The number of particles */
386  int numParticles = path.numBeadsAtSlice(slice);
387 
388  dVec sep; // The spatial separation between beads.
389  double r; // Distance between beads
390  beadLocator bead1,bead2;
391  bead1[0] = bead2[0] = slice;
392 
393  /* No cutoff */
394  for (bead1[1] = 0; bead1[1] < numParticles; bead1[1]++) {
395  /* The loop over all other particles, to find the total wavefunction */
396  for (bead2[1] = bead1[1]+1; bead2[1] < numParticles; bead2[1]++) {
397  sep = path.getSeparation(bead2,bead1);
398  r = sqrt(dot(sep,sep));
399  psiT *= PsiTrial(r);
400  } // bead2
401  } // bead1
402 
403  return psiT;
404 }
double R_LL_wfn() const
Get Lieb-Liniger length scale.
Definition: constants.h:117
double k_LL_wfn() const
Get Lieb-Liniger wave number.
Definition: constants.h:118
double L() const
Get maximum side length.
Definition: constants.h:52
JastrowWaveFunction(const Path &, LookupTable &_lookup, string _name="Jastrow")
Constructor.
double gradSqPsiTrial(const int)
The value of the N-body trial wave function.
double delSqPsiTrial(const double r)
The derivative of psi over psi.
double PsiTrial(const double)
The value of the 2-body trial wave function.
double delPsiTrial(const double r)
The derivative of psi over psi.
~JastrowWaveFunction()
Destructor.
double gradSqPsiTrial(const int)
The value of the N-body trial wave function.
double PsiTrial(const double)
The value of the 2-body trial wave function.
double delPsiTrial(const double r)
The derivative of psi over psi.
LiebLinigerWaveFunction(const Path &, LookupTable &_lookup, string _name="LiebLiniger")
Constructor.
~LiebLinigerWaveFunction()
Destructor.
double delSqPsiTrial(const double r)
The derivative of psi over psi.
The particle (bead) lookup table.
Definition: lookuptable.h:29
Array< beadLocator, 1 > beadList
The cutoff dynamic list of interacting beads.
Definition: lookuptable.h:40
void updateInteractionList(const Path &, const beadLocator &)
Update the NN lookup table and the array of beadLocators containing all beads which 'interact' with t...
int numBeads
The cutoff number of active beads in beadList;.
Definition: lookuptable.h:43
The space-time trajectories.
Definition: path.h:29
dVec getSeparation(const beadLocator &, const beadLocator &) const
Return the separation vector between two particles in the same timeslice.
Definition: path.h:173
Array< int, 1 > numBeadsAtSlice
The number of active beads at a given time slice.
Definition: path.h:48
Worm worm
Details on the worm.
Definition: path.h:44
double PsiTrial(const int)
The value of the trial wave function.
SechWaveFunction(const Path &, LookupTable &_lookup, string _name="SHO sech")
Constructor.
~SechWaveFunction()
Destructor.
SutherlandWaveFunction(const Path &, LookupTable &_lookup, double, string _name="Sutherland")
Constructor.
~SutherlandWaveFunction()
Destructor.
double PsiTrial(const double r)
The 2-body trial wavefunction.
Definition: wavefunction.h:143
Holds a base class that all trial wave function classes will be derived from.
Definition: wavefunction.h:26
LookupTable & lookup
We need a non-constant reference for updates.
Definition: wavefunction.h:45
virtual ~WaveFunctionBase()
Empty base constructor.
WaveFunctionBase(const Path &, LookupTable &_lookup, string _name="constant")
Setup the path data members for the constant trial wavefunction.
const Path & path
A reference to the paths.
Definition: wavefunction.h:41
int beadOn(int, int) const
Safely get a bead (int indexed)
Definition: worm.h:148
TinyVector< int, 2 > beadLocator
time-slice,bead-number world line index
Definition: common.h:117
TinyVector< double, NDIM > dVec
A NDIM-vector of type double.
Definition: common.h:111
ConstantParameters * constants()
Global public access to the constants.
Definition: constants.h:201
Path class definition.
Action class definitions.