QMCPACK
Public Types | Public Member Functions | Public Attributes | List of all members
Quadrature3D< T > Struct Template Reference
+ Collaboration diagram for Quadrature3D< T >:

Public Types

enum  SymmType { SINGLE, TETRA, OCTA, ICOSA }
 
using RealType = T
 
using PosType = TinyVector< T, 3 >
 
using mRealType = OHMMS_PRECISION_FULL
 

Public Member Functions

 Quadrature3D (int rule, bool request_abort=true)
 
void CheckQuadratureRule (int lexact)
 
void CheckQuadratureRuleReal (int lexact)
 

Public Attributes

int nk
 
SymmType symmetry
 
int lexact
 
RealType A
 
RealType B
 
RealType C
 
RealType D
 
std::vector< PosType > xyz_m
 
std::vector< RealType > weight_m
 
bool quad_ok
 
const bool fail_abort
 

Detailed Description

template<class T>
struct qmcplusplus::Quadrature3D< T >

Definition at line 26 of file Quadrature.h.

Member Typedef Documentation

◆ mRealType

using mRealType = OHMMS_PRECISION_FULL

Definition at line 30 of file Quadrature.h.

◆ PosType

using PosType = TinyVector<T, 3>

Definition at line 29 of file Quadrature.h.

◆ RealType

using RealType = T

Definition at line 28 of file Quadrature.h.

Member Enumeration Documentation

◆ SymmType

enum SymmType
Enumerator
SINGLE 
TETRA 
OCTA 
ICOSA 

Definition at line 33 of file Quadrature.h.

34  {
35  SINGLE,
36  TETRA,
37  OCTA,
38  ICOSA
39  } SymmType;

Constructor & Destructor Documentation

◆ Quadrature3D()

Quadrature3D ( int  rule,
bool  request_abort = true 
)
inline

Definition at line 48 of file Quadrature.h.

References Quadrature3D< T >::A, qmcplusplus::abs(), qmcplusplus::acos(), qmcplusplus::atan(), Quadrature3D< T >::B, Quadrature3D< T >::C, Quadrature3D< T >::CheckQuadratureRule(), Quadrature3D< T >::CheckQuadratureRuleReal(), qmcplusplus::cos(), Quadrature3D< T >::D, qmcplusplus::dot(), qmcplusplus::Units::charge::e, ERRORMSG, Quadrature3D< T >::ICOSA, Quadrature3D< T >::lexact, Quadrature3D< T >::nk, Quadrature3D< T >::OCTA, qmcplusplus::Units::time::s, qmcplusplus::sin(), Quadrature3D< T >::SINGLE, qmcplusplus::sqrt(), Quadrature3D< T >::symmetry, Quadrature3D< T >::TETRA, Quadrature3D< T >::weight_m, and Quadrature3D< T >::xyz_m.

48  : quad_ok(true), fail_abort(request_abort)
49  {
50  A = B = C = D = 0;
51  switch (rule)
52  {
53  case 1:
54  nk = 1;
55  symmetry = SINGLE;
56  lexact = 0;
57  A = 1.0;
58  break;
59  case 2:
60  nk = 4;
61  symmetry = TETRA;
62  lexact = 2;
63  A = 0.25;
64  break;
65  case 3:
66  nk = 6;
67  symmetry = OCTA;
68  lexact = 3;
69  A = 1.0 / 6.0;
70  break;
71  case 4:
72  nk = 12;
73  symmetry = ICOSA;
74  lexact = 5;
75  A = 1.0 / 12.0;
76  B = 1.0 / 12.0;
77  break;
78  case 5:
79  nk = 18;
80  symmetry = OCTA;
81  lexact = 5;
82  A = 1.0 / 30.0;
83  B = 1.0 / 15.0;
84  break;
85  case 6:
86  nk = 26;
87  symmetry = OCTA;
88  lexact = 7;
89  A = 1.0 / 21.0;
90  B = 4.0 / 105.0;
91  C = 27.0 / 840.0;
92  break;
93  case 7:
94  nk = 32;
95  symmetry = ICOSA;
96  lexact = 9;
97  A = 5.0 / 168.0;
98  B = 5.0 / 168.0;
99  C = 27.0 / 840.0;
100  break;
101  case 8:
102  nk = 50;
103  symmetry = OCTA;
104  lexact = 11;
105  A = 4.0 / 315.0;
106  B = 64.0 / 2835.0;
107  C = 27.0 / 1280.0;
108  D = 14641.0 / 725760.0;
109  break;
110  default:
111  ERRORMSG("Unrecognized spherical quadrature rule " << rule << ".");
112  abort();
113  }
114  // First, build a_i, b_i, and c_i points
115  std::vector<PosType> a, b, c, d;
116  RealType p = 1.0 / std::sqrt(2.0);
117  RealType q = 1.0 / std::sqrt(3.0);
118  RealType r = 1.0 / std::sqrt(11.0);
119  RealType s = 3.0 / std::sqrt(11.0);
120  if (symmetry == SINGLE)
121  {
122  a.push_back(PosType(1.0, 0.0, 0.0));
123  }
124  else if (symmetry == TETRA)
125  {
126  a.push_back(PosType(q, q, q));
127  a.push_back(PosType(q, -q, -q));
128  a.push_back(PosType(-q, q, -q));
129  a.push_back(PosType(-q, -q, q));
130  }
131  else if (symmetry == OCTA)
132  {
133  a.push_back(PosType(1.0, 0.0, 0.0));
134  a.push_back(PosType(-1.0, 0.0, 0.0));
135  a.push_back(PosType(0.0, 1.0, 0.0));
136  a.push_back(PosType(0.0, -1.0, 0.0));
137  a.push_back(PosType(0.0, 0.0, 1.0));
138  a.push_back(PosType(0.0, 0.0, -1.0));
139  b.push_back(PosType(p, p, 0.0));
140  b.push_back(PosType(p, -p, 0.0));
141  b.push_back(PosType(-p, p, 0.0));
142  b.push_back(PosType(-p, -p, 0.0));
143  b.push_back(PosType(p, 0.0, p));
144  b.push_back(PosType(p, 0.0, -p));
145  b.push_back(PosType(-p, 0.0, p));
146  b.push_back(PosType(-p, 0.0, -p));
147  b.push_back(PosType(0.0, p, p));
148  b.push_back(PosType(0.0, p, -p));
149  b.push_back(PosType(0.0, -p, p));
150  b.push_back(PosType(0.0, -p, -p));
151  c.push_back(PosType(q, q, q));
152  c.push_back(PosType(q, q, -q));
153  c.push_back(PosType(q, -q, q));
154  c.push_back(PosType(q, -q, -q));
155  c.push_back(PosType(-q, q, q));
156  c.push_back(PosType(-q, q, -q));
157  c.push_back(PosType(-q, -q, q));
158  c.push_back(PosType(-q, -q, -q));
159  d.push_back(PosType(r, r, s));
160  d.push_back(PosType(r, r, -s));
161  d.push_back(PosType(r, -r, s));
162  d.push_back(PosType(r, -r, -s));
163  d.push_back(PosType(-r, r, s));
164  d.push_back(PosType(-r, r, -s));
165  d.push_back(PosType(-r, -r, s));
166  d.push_back(PosType(-r, -r, -s));
167  d.push_back(PosType(r, s, r));
168  d.push_back(PosType(r, s, -r));
169  d.push_back(PosType(r, -s, r));
170  d.push_back(PosType(r, -s, -r));
171  d.push_back(PosType(-r, s, r));
172  d.push_back(PosType(-r, s, -r));
173  d.push_back(PosType(-r, -s, r));
174  d.push_back(PosType(-r, -s, -r));
175  d.push_back(PosType(s, r, r));
176  d.push_back(PosType(s, r, -r));
177  d.push_back(PosType(s, -r, r));
178  d.push_back(PosType(s, -r, -r));
179  d.push_back(PosType(-s, r, r));
180  d.push_back(PosType(-s, r, -r));
181  d.push_back(PosType(-s, -r, r));
182  d.push_back(PosType(-s, -r, -r));
183  }
184  else if (symmetry == ICOSA)
185  {
186  mRealType t, p; // theta and phi
187  // a points
188  t = 0.0;
189  p = 0.0;
190  a.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
191  t = M_PI;
192  p = 0.0;
193  a.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
194  // b points
195  for (int k = 0; k < 5; k++)
196  {
197  t = std::atan(2.0);
198  p = (mRealType)(2 * k + 0) * M_PI / 5.0;
199  b.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
200  t = M_PI - std::atan(2.0);
201  p = (mRealType)(2 * k + 1) * M_PI / 5.0;
202  b.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
203  }
204  // c points
205  mRealType t1 = std::acos((2.0 + std::sqrt(5.0)) / std::sqrt(15.0 + 6.0 * std::sqrt(5.0)));
206  mRealType t2 = std::acos(1.0 / std::sqrt(15.0 + 6.0 * std::sqrt(5.0)));
207  for (int k = 0; k < 5; k++)
208  {
209  t = t1;
210  p = (mRealType)(2 * k + 1) * M_PI / 5.0;
211  c.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
212  t = t2;
213  p = (mRealType)(2 * k + 1) * M_PI / 5.0;
214  c.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
215  t = M_PI - t1;
216  p = (mRealType)(2 * k + 0) * M_PI / 5.0;
217  c.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
218  t = M_PI - t2;
219  p = (mRealType)(2 * k + 0) * M_PI / 5.0;
220  c.push_back(PosType(std::cos(t), std::sin(t) * std::cos(p), std::sin(t) * std::sin(p)));
221  }
222  }
223  // Now, construct rule
224  if (std::abs(A) > 1.0e-10)
225  for (int i = 0; i < a.size(); i++)
226  {
227  xyz_m.push_back(a[i]);
228  weight_m.push_back(A);
229  }
230  if (std::abs(B) > 1.0e-10)
231  for (int i = 0; i < b.size(); i++)
232  {
233  xyz_m.push_back(b[i]);
234  weight_m.push_back(B);
235  }
236  if (std::abs(C) > 1.0e-10)
237  for (int i = 0; i < c.size(); i++)
238  {
239  xyz_m.push_back(c[i]);
240  weight_m.push_back(C);
241  }
242  if (std::abs(D) > 1.0e-10)
243  for (int i = 0; i < d.size(); i++)
244  {
245  xyz_m.push_back(d[i]);
246  weight_m.push_back(D);
247  }
248  // Finally, check the rule for correctness
249  assert(xyz_m.size() == nk);
250  assert(weight_m.size() == nk);
251  double wSum = 0.0;
252 
253  for (int k = 0; k < nk; k++)
254  {
255  PosType r = xyz_m[k];
256  double nrm = dot(r, r);
257  assert(std::abs(nrm - 1.0) < 2 * std::numeric_limits<float>::epsilon());
258  wSum += weight_m[k];
259  //cout << pp_nonloc->xyz_m[k] << " " << pp_nonloc->weight_m[k] << std::endl;
260  }
261  assert(std::abs(wSum - 1.0) < 2 * std::numeric_limits<float>::epsilon());
262  // Check the quadrature rule
263  // using complex spherical harmonics
264  CheckQuadratureRule(lexact);
265  // using real spherical harmonics
266  CheckQuadratureRuleReal(lexact);
267  }
ERRORMSG
#define ERRORMSG(msg)
Definition: OutputManager.h:86
qmcplusplus::Quadrature3D::CheckQuadratureRuleReal
void CheckQuadratureRuleReal(int lexact)
Definition: Quadrature.h:303
qmcplusplus::abs
MakeReturn< UnaryNode< FnFabs, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t abs(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:88
qmcplusplus::sin
MakeReturn< UnaryNode< FnSin, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t sin(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:120
PosType
QMCTraits::PosType PosType
Definition: qmcfinitesize.cpp:48
qmcplusplus::cos
MakeReturn< UnaryNode< FnCos, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t cos(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:64
qmcplusplus::acos
MakeReturn< UnaryNode< FnArcCos, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t acos(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:32
qmcplusplus::atan
MakeReturn< UnaryNode< FnArcTan, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t atan(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:48
qmcplusplus::Quadrature3D::mRealType
OHMMS_PRECISION_FULL mRealType
Definition: Quadrature.h:30
qmcplusplus::Quadrature3D::CheckQuadratureRule
void CheckQuadratureRule(int lexact)
Definition: Quadrature.h:269
qmcplusplus::Quadrature3D::xyz_m
std::vector< PosType > xyz_m
Definition: Quadrature.h:43
qmcplusplus::Quadrature3D::PosType
TinyVector< T, 3 > PosType
Definition: Quadrature.h:29
qmcplusplus::sqrt
MakeReturn< UnaryNode< FnSqrt, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t sqrt(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:136
RealType
QMCTraits::RealType RealType
Definition: qmcfinitesize.cpp:47
qmcplusplus::Quadrature3D::weight_m
std::vector< RealType > weight_m
Definition: Quadrature.h:44
qmcplusplus::dot
Tensor< typename BinaryReturn< T1, T2, OpMultiply >::Type_t, D > dot(const AntiSymTensor< T1, D > &lhs, const AntiSymTensor< T2, D > &rhs)
Definition: AntiSymTensor.h:520

Member Function Documentation

◆ CheckQuadratureRule()

void CheckQuadratureRule ( int  lexact)
inline

Definition at line 269 of file Quadrature.h.

References qmcplusplus::abs(), APP_ABORT, qmcplusplus::app_error(), qmcplusplus::conj(), Quadrature3D< T >::fail_abort, qmcplusplus::imag(), Quadrature3D< T >::lexact, Quadrature3D< T >::quad_ok, qmcplusplus::real(), Quadrature3D< T >::weight_m, Quadrature3D< T >::xyz_m, and qmcplusplus::Ylm().

Referenced by Quadrature3D< T >::Quadrature3D().

270  {
271  std::vector<PosType>& grid = xyz_m;
272  std::vector<RealType>& w = weight_m;
273  for (int l1 = 0; l1 <= lexact; l1++)
274  for (int l2 = 0; l2 <= (lexact - l1); l2++)
275  for (int m1 = -l1; m1 <= l1; m1++)
276  for (int m2 = -l2; m2 <= l2; m2++)
277  {
278  std::complex<mRealType> sum(0.0, 0.0);
279  for (int k = 0; k < grid.size(); k++)
280  {
281  std::complex<mRealType> v1 = Ylm(l1, m1, grid[k]);
282  std::complex<mRealType> v2 = Ylm(l2, m2, grid[k]);
283  sum += 4.0 * M_PI * w[k] * qmcplusplus::conj(v1) * v2;
284  }
285  mRealType re = real(sum);
286  mRealType im = imag(sum);
287  if ((l1 == l2) && (m1 == m2))
288  re -= 1.0;
289  if ((std::abs(im) > 7 * std::numeric_limits<float>::epsilon()) ||
290  (std::abs(re) > 7 * std::numeric_limits<float>::epsilon()))
291  {
292  app_error() << "Broken spherical quadrature for " << grid.size() << "-point rule.\n" << std::endl;
293  app_error() << " Should be zero: Real part = " << re << " Imaginary part = " << im << std::endl;
294  quad_ok = false;
295  if (fail_abort)
296  APP_ABORT("Give up");
297  }
298  // fprintf (stderr, "(l1,m1,l2m,m2) = (%2d,%2d,%2d,%2d) sum = (%20.16f %20.16f)\n",
299  // l1, m1, l2, m2, real(sum), imag(sum));
300  }
301  }
qmcplusplus::Ylm
std::complex< T > Ylm(int l, int m, const TinyVector< T, 3 > &r)
calculates Ylm param[in] l angular momentum param[in] m magnetic quantum number param[in] r position ...
Definition: Ylm.h:89
qmcplusplus::abs
MakeReturn< UnaryNode< FnFabs, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t abs(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:88
qmcplusplus::app_error
std::ostream & app_error()
Definition: OutputManager.h:67
qmcplusplus::real
float real(const float &c)
real part of a scalar. Cannot be replaced by std::real due to AFQMC specific needs.
Definition: complex_help.hpp:86
qmcplusplus::imag
float imag(const float &c)
imaginary part of a scalar. Cannot be replaced by std::imag due to AFQMC specific needs...
Definition: complex_help.hpp:91
qmcplusplus::Quadrature3D::mRealType
OHMMS_PRECISION_FULL mRealType
Definition: Quadrature.h:30
APP_ABORT
#define APP_ABORT(msg)
Widely used but deprecated fatal error macros from legacy code.
Definition: AppAbort.h:27
qmcplusplus::Quadrature3D::xyz_m
std::vector< PosType > xyz_m
Definition: Quadrature.h:43
qmcplusplus::conj
float conj(const float &c)
Workaround to allow conj on scalar to return real instead of complex.
Definition: complex_help.hpp:96
qmcplusplus::Quadrature3D::weight_m
std::vector< RealType > weight_m
Definition: Quadrature.h:44

◆ CheckQuadratureRuleReal()

void CheckQuadratureRuleReal ( int  lexact)
inline

Definition at line 303 of file Quadrature.h.

References qmcplusplus::abs(), APP_ABORT, qmcplusplus::app_error(), Quadrature3D< T >::fail_abort, Quadrature3D< T >::lexact, Quadrature3D< T >::quad_ok, Quadrature3D< T >::weight_m, Quadrature3D< T >::xyz_m, and qmcplusplus::Ylm().

Referenced by Quadrature3D< T >::Quadrature3D().

304  {
305  std::vector<PosType>& grid = xyz_m;
306  std::vector<RealType>& w = weight_m;
307  SoaSphericalTensor<RealType> Ylm(lexact);
308  for (int l1 = 0; l1 <= lexact; l1++)
309  for (int l2 = 0; l2 <= (lexact - l1); l2++)
310  for (int m1 = -l1; m1 <= l1; m1++)
311  for (int m2 = -l2; m2 <= l2; m2++)
312  {
313  mRealType sum(0.0);
314  for (int k = 0; k < grid.size(); k++)
315  {
316  Ylm.evaluateV(grid[k][0], grid[k][1], grid[k][2]);
317  const RealType* Ylm_v = Ylm[0];
318  RealType v1 = Ylm_v[Ylm.index(l1, m1)];
319  RealType v2 = Ylm_v[Ylm.index(l2, m2)];
320  sum += 4.0 * M_PI * w[k] * v1 * v2;
321  }
322  if ((l1 == l2) && (m1 == m2))
323  sum -= 1.0;
324  if (std::abs(sum) > 16 * std::numeric_limits<float>::epsilon())
325  {
326  app_error() << "Broken real spherical quadrature for " << grid.size() << "-point rule.\n" << std::endl;
327  app_error() << " Should be zero: " << sum << std::endl;
328  quad_ok = false;
329  if (fail_abort)
330  APP_ABORT("Give up");
331  }
332  }
333  }
qmcplusplus::Ylm
std::complex< T > Ylm(int l, int m, const TinyVector< T, 3 > &r)
calculates Ylm param[in] l angular momentum param[in] m magnetic quantum number param[in] r position ...
Definition: Ylm.h:89
qmcplusplus::abs
MakeReturn< UnaryNode< FnFabs, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t abs(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:88
qmcplusplus::app_error
std::ostream & app_error()
Definition: OutputManager.h:67
qmcplusplus::Quadrature3D::mRealType
OHMMS_PRECISION_FULL mRealType
Definition: Quadrature.h:30
APP_ABORT
#define APP_ABORT(msg)
Widely used but deprecated fatal error macros from legacy code.
Definition: AppAbort.h:27
qmcplusplus::Quadrature3D::xyz_m
std::vector< PosType > xyz_m
Definition: Quadrature.h:43
RealType
QMCTraits::RealType RealType
Definition: qmcfinitesize.cpp:47
qmcplusplus::Quadrature3D::weight_m
std::vector< RealType > weight_m
Definition: Quadrature.h:44

Member Data Documentation

◆ A

RealType A

Definition at line 42 of file Quadrature.h.

Referenced by Quadrature3D< T >::Quadrature3D().

◆ B

RealType B

Definition at line 42 of file Quadrature.h.

Referenced by Quadrature3D< T >::Quadrature3D().

◆ C

RealType C

Definition at line 42 of file Quadrature.h.

Referenced by Quadrature3D< T >::Quadrature3D().

◆ D

RealType D

Definition at line 42 of file Quadrature.h.

Referenced by Quadrature3D< T >::Quadrature3D().

◆ fail_abort

const bool fail_abort

Definition at line 46 of file Quadrature.h.

Referenced by Quadrature3D< T >::CheckQuadratureRule(), and Quadrature3D< T >::CheckQuadratureRuleReal().

◆ lexact

int lexact

Definition at line 41 of file Quadrature.h.

Referenced by Quadrature3D< T >::CheckQuadratureRule(), Quadrature3D< T >::CheckQuadratureRuleReal(), and Quadrature3D< T >::Quadrature3D().

◆ nk

int nk

Definition at line 32 of file Quadrature.h.

Referenced by Quadrature3D< T >::Quadrature3D(), and ECPComponentBuilder::SetQuadratureRule().

◆ quad_ok

bool quad_ok

Definition at line 45 of file Quadrature.h.

Referenced by Quadrature3D< T >::CheckQuadratureRule(), Quadrature3D< T >::CheckQuadratureRuleReal(), and qmcplusplus::TEST_CASE().

◆ symmetry

SymmType symmetry

Definition at line 40 of file Quadrature.h.

Referenced by Quadrature3D< T >::Quadrature3D().

◆ weight_m

std::vector<RealType> weight_m

Definition at line 44 of file Quadrature.h.

Referenced by Quadrature3D< T >::CheckQuadratureRule(), Quadrature3D< T >::CheckQuadratureRuleReal(), Quadrature3D< T >::Quadrature3D(), and ECPComponentBuilder::SetQuadratureRule().

◆ xyz_m

std::vector<PosType> xyz_m

Definition at line 43 of file Quadrature.h.

Referenced by Quadrature3D< T >::CheckQuadratureRule(), Quadrature3D< T >::CheckQuadratureRuleReal(), Quadrature3D< T >::Quadrature3D(), and ECPComponentBuilder::SetQuadratureRule().

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