QMCPACK
FiniteDifference Class Reference
+ Inheritance diagram for FiniteDifference:
+ Collaboration diagram for FiniteDifference:

Classes

struct  PositionChange
 

Public Types

enum  FiniteDiffType { FiniteDiff_LowOrder, FiniteDiff_Richardson }
 
using PosChangeVector = std::vector< PositionChange >
 
using ValueVector = std::vector< ValueType >
 
- Public Types inherited from QMCTraits
enum  { DIM = OHMMS_DIM, DIM_VGL = OHMMS_DIM + 2 }
 
using QTBase = QMCTypes< OHMMS_PRECISION, DIM >
 
using QTFull = QMCTypes< OHMMS_PRECISION_FULL, DIM >
 
using RealType = QTBase::RealType
 
using ComplexType = QTBase::ComplexType
 
using ValueType = QTBase::ValueType
 
using PosType = QTBase::PosType
 
using GradType = QTBase::GradType
 
using TensorType = QTBase::TensorType
 
using IndexType = OHMMS_INDEXTYPE
 define other types More...
 
using FullPrecRealType = QTFull::RealType
 
using FullPrecValueType = QTFull::ValueType
 
using PropertySetType = RecordNamedProperty< FullPrecRealType >
 define PropertyList_t More...
 
using PtclGrpIndexes = std::vector< std::pair< int, int > >
 

Public Member Functions

 FiniteDifference (size_t ndim_in, FiniteDiffType fd_type=FiniteDiff_Richardson)
 
void finiteDifferencePoints (RealType delta, MCWalkerConfiguration &W, PosChangeVector &positions)
 Generate points to evaluate. More...
 
void computeFiniteDiff (RealType delta, PosChangeVector &positions, ValueVector &values, ParticleSet::ParticleGradient &G_fd, ParticleSet::ParticleLaplacian &L_fd)
 Compute finite difference after log psi is computed for each point. More...
 
void computeFiniteDiffLowOrder (RealType delta, PosChangeVector &positions, ValueVector &values, ParticleSet::ParticleGradient &G_fd, ParticleSet::ParticleLaplacian &L_fd)
 
void computeFiniteDiffRichardson (RealType delta, PosChangeVector &positions, ValueVector &values, ParticleSet::ParticleGradient &G_fd, ParticleSet::ParticleLaplacian &L_fd)
 

Public Attributes

const size_t ndim
 
int m_RichardsonSize
 
FiniteDiffType m_fd_type
 

Detailed Description

Definition at line 308 of file WaveFunctionTester.cpp.


Class Documentation

◆ qmcplusplus::FiniteDifference::PositionChange

struct qmcplusplus::FiniteDifference::PositionChange

Definition at line 325 of file WaveFunctionTester.cpp.

+ Collaboration diagram for FiniteDifference::PositionChange:
Class Members
int index
PosType r

Member Typedef Documentation

◆ PosChangeVector

using PosChangeVector = std::vector<PositionChange>

Definition at line 330 of file WaveFunctionTester.cpp.

◆ ValueVector

using ValueVector = std::vector<ValueType>

Definition at line 331 of file WaveFunctionTester.cpp.

Member Enumeration Documentation

◆ FiniteDiffType

Enumerator
FiniteDiff_LowOrder 
FiniteDiff_Richardson 

Definition at line 311 of file WaveFunctionTester.cpp.

Constructor & Destructor Documentation

◆ FiniteDifference()

FiniteDifference ( size_t  ndim_in,
FiniteDiffType  fd_type = FiniteDiff_Richardson 
)
inline

Definition at line 316 of file WaveFunctionTester.cpp.

Member Function Documentation

◆ computeFiniteDiff()

void computeFiniteDiff ( RealType  delta,
PosChangeVector positions,
ValueVector values,
ParticleSet::ParticleGradient G_fd,
ParticleSet::ParticleLaplacian L_fd 
)

Compute finite difference after log psi is computed for each point.

Definition at line 404 of file WaveFunctionTester.cpp.

References FiniteDifference::computeFiniteDiffLowOrder(), FiniteDifference::computeFiniteDiffRichardson(), FiniteDifference::FiniteDiff_LowOrder, FiniteDifference::FiniteDiff_Richardson, and FiniteDifference::m_fd_type.

Referenced by WaveFunctionTester::checkGradientAtConfiguration(), and WaveFunctionTester::computeNumericalGrad().

409 {
410  assert(positions.size() == values.size());
411  if (positions.size() == 0)
412  return;
413 
414  ValueType logpsi = values[0];
415 
417  {
418  computeFiniteDiffLowOrder(delta, positions, values, G_fd, L_fd);
419  }
420  else if (m_fd_type == FiniteDiff_Richardson)
421  {
422  computeFiniteDiffRichardson(delta, positions, values, G_fd, L_fd);
423  }
424 }
QTBase::ValueType ValueType
Definition: Configuration.h:60
void computeFiniteDiffLowOrder(RealType delta, PosChangeVector &positions, ValueVector &values, ParticleSet::ParticleGradient &G_fd, ParticleSet::ParticleLaplacian &L_fd)
void computeFiniteDiffRichardson(RealType delta, PosChangeVector &positions, ValueVector &values, ParticleSet::ParticleGradient &G_fd, ParticleSet::ParticleLaplacian &L_fd)

◆ computeFiniteDiffLowOrder()

void computeFiniteDiffLowOrder ( RealType  delta,
PosChangeVector positions,
ValueVector values,
ParticleSet::ParticleGradient G_fd,
ParticleSet::ParticleLaplacian L_fd 
)

Definition at line 426 of file WaveFunctionTester.cpp.

References FiniteDifference::ndim.

Referenced by FiniteDifference::computeFiniteDiff().

431 {
432  ValueType logpsi = values[0];
433 
434  // lowest order derivative formula
435  ValueType c1 = 1.0 / delta / 2.0;
436  ValueType c2 = 1.0 / delta / delta;
437 
438  const RealType twoD(2 * ndim);
439  const int pt_per_deriv = 2; // number of points per derivative
440  for (int pt_i = 1; pt_i < values.size(); pt_i += pt_per_deriv * ndim)
441  {
442  GradType g0;
443  ValueType lap0 = 0.0;
444  for (int idim = 0; idim < ndim; idim++)
445  {
446  int idx = pt_i + idim * pt_per_deriv;
447  ValueType logpsi_m = values[idx];
448  ValueType logpsi_p = values[idx + 1];
449 
450  g0[idim] = logpsi_p - logpsi_m;
451  lap0 += logpsi_p + logpsi_m;
452  }
453 
454  int iat = positions[pt_i].index;
455  GradType g = c1 * g0;
456  G_fd[iat] = g;
457 
458  ValueType lap = c2 * (lap0 - twoD * logpsi);
459  //ValueType lap = c2*(lap0 - 2.0*OHMMS_DIM*logpsi);
460  L_fd[iat] = lap;
461  }
462 }
QTBase::GradType GradType
Definition: Configuration.h:62
QTBase::ValueType ValueType
Definition: Configuration.h:60
QMCTraits::RealType RealType

◆ computeFiniteDiffRichardson()

void computeFiniteDiffRichardson ( RealType  delta,
PosChangeVector positions,
ValueVector values,
ParticleSet::ParticleGradient G_fd,
ParticleSet::ParticleLaplacian L_fd 
)

Definition at line 468 of file WaveFunctionTester.cpp.

References qmcplusplus::abs(), qmcplusplus::Units::charge::e, FiniteDifference::m_RichardsonSize, FiniteDifference::ndim, and norm().

Referenced by FiniteDifference::computeFiniteDiff().

473 {
474  RealType tol = 1e-7;
475  ValueType logpsi = values[0];
476 
477  const int pt_per_deriv = 2 * (m_RichardsonSize + 1); // number of points per derivative
478  for (int pt_i = 1; pt_i < values.size(); pt_i += pt_per_deriv * ndim)
479  {
480  GradType g0;
481  GradType gmin;
482  ValueType lmin;
483  std::vector<GradType> g_base(m_RichardsonSize + 1);
484  std::vector<GradType> g_rich(m_RichardsonSize + 1);
485  std::vector<GradType> g_prev(m_RichardsonSize + 1);
486 
487  std::vector<ValueType> l_base(m_RichardsonSize + 1);
488  std::vector<ValueType> l_rich(m_RichardsonSize + 1);
489  std::vector<ValueType> l_prev(m_RichardsonSize + 1);
490 
491  // Initial gradients and Laplacians at different deltas.
492  RealType dd = delta;
493  const ValueType ctwo(2);
494  for (int inr = 0; inr < m_RichardsonSize + 1; inr++)
495  {
496  RealType twodd = 2 * dd;
497  RealType ddsq = dd * dd;
498  l_base[inr] = 0.0;
499  for (int idim = 0; idim < ndim; idim++)
500  {
501  int idx = pt_i + idim * pt_per_deriv + 2 * inr;
502  ValueType logpsi_m = values[idx];
503  ValueType logpsi_p = values[idx + 1];
504  g_base[inr][idim] = (logpsi_p - logpsi_m) / twodd;
505  l_base[inr] += (logpsi_p + logpsi_m - ctwo * logpsi) / ddsq;
506  //g_base[inr][idim] = (logpsi_p - logpsi_m)/dd/2.0;
507  //l_base[inr] += (logpsi_p + logpsi_m - 2.0*logpsi)/(dd*dd);
508  }
509  dd = dd / 2;
510  }
511 
512  // Gradient
513 
514  g_prev[0] = g_base[0];
515  RealType fac = 1;
516  bool found = false;
517  for (int inr = 1; inr < m_RichardsonSize + 1; inr++)
518  {
519  g_rich[0] = g_base[inr];
520 
521  fac *= 4;
522  for (int j = 1; j < inr + 1; j++)
523  {
524  g_rich[j] = g_rich[j - 1] + (g_rich[j - 1] - g_prev[j - 1]) / (fac - 1);
525  }
526 
527  RealType err1 = 0.0;
528  RealType norm = 0.0;
529  for (int idim = 0; idim < ndim; idim++)
530  {
531  err1 += std::abs(g_rich[inr][idim] - g_prev[inr - 1][idim]);
532  norm += std::abs(g_prev[inr - 1][idim]);
533  }
534 
535  RealType err_rel = err1 / norm;
536 
537  // Not sure about the best stopping criteria
538  if (err_rel < tol)
539  {
540  gmin = g_rich[inr];
541  found = true;
542  break;
543  }
544  g_prev = g_rich;
545  }
546 
547  if (!found)
548  {
549  gmin = g_rich[m_RichardsonSize];
550  }
551 
552 
553  // Laplacian
554  // TODO: eliminate the copied code between the gradient and Laplacian
555  // computations.
556 
557  l_prev[0] = l_base[0];
558 
559  fac = 1;
560  found = false;
561  for (int inr = 1; inr < m_RichardsonSize + 1; inr++)
562  {
563  l_rich[0] = l_base[inr];
564 
565  fac *= 4;
566  for (int j = 1; j < inr + 1; j++)
567  {
568  l_rich[j] = l_rich[j - 1] + (l_rich[j - 1] - l_prev[j - 1]) / (fac - 1);
569  }
570 
571  RealType err1 = std::abs(l_rich[inr] - l_prev[inr - 1]);
572  RealType err_rel = std::abs(err1 / l_prev[inr - 1]);
573 
574  if (err_rel < tol)
575  {
576  lmin = l_rich[inr];
577  found = true;
578  break;
579  }
580  l_prev = l_rich;
581  }
582 
583  if (!found)
584  {
585  lmin = l_rich[m_RichardsonSize];
586  }
587 
588  int iat = positions[pt_i].index;
589  G_fd[iat] = gmin;
590  L_fd[iat] = lmin;
591  }
592 }
QTBase::GradType GradType
Definition: Configuration.h:62
MakeReturn< UnaryNode< FnFabs, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t abs(const Vector< T1, C1 > &l)
double norm(const zVec &c)
Definition: VectorOps.h:118
QTBase::ValueType ValueType
Definition: Configuration.h:60
QMCTraits::RealType RealType

◆ finiteDifferencePoints()

void finiteDifferencePoints ( RealType  delta,
MCWalkerConfiguration W,
PosChangeVector positions 
)

Generate points to evaluate.

Definition at line 358 of file WaveFunctionTester.cpp.

References FiniteDifference::FiniteDiff_Richardson, ParticleSet::getTotalNum(), FiniteDifference::PositionChange::index, FiniteDifference::m_fd_type, FiniteDifference::m_RichardsonSize, FiniteDifference::ndim, ParticleSet::R, and FiniteDifference::PositionChange::r.

Referenced by WaveFunctionTester::checkGradientAtConfiguration(), and WaveFunctionTester::computeNumericalGrad().

359 {
360  // First position is the central point
361  PositionChange p;
362  p.index = 0;
363  p.r = W.R[0];
364  positions.push_back(p);
365 
366  int nat = W.getTotalNum();
367  for (int iat = 0; iat < nat; iat++)
368  {
369  PositionChange p;
370  p.index = iat;
371  PosType r0 = W.R[iat];
372 
373  for (int idim = 0; idim < ndim; idim++)
374  {
375  p.r = r0;
376  p.r[idim] = r0[idim] - delta;
377  positions.push_back(p);
378 
379  p.r = r0;
380  p.r[idim] = r0[idim] + delta;
381  positions.push_back(p);
382 
384  {
385  RealType dd = delta / 2;
386  for (int nr = 0; nr < m_RichardsonSize; nr++)
387  {
388  p.r = r0;
389  p.r[idim] = r0[idim] - dd;
390  positions.push_back(p);
391 
392  p.r = r0;
393  p.r[idim] = r0[idim] + dd;
394  positions.push_back(p);
395 
396  dd = dd / 2;
397  }
398  }
399  }
400  }
401 }
QMCTraits::PosType PosType
QMCTraits::RealType RealType

Member Data Documentation

◆ m_fd_type

◆ m_RichardsonSize

◆ ndim


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