QMCPACK
Public Types | Public Member Functions | Public Attributes | List of all members
SphericalTensor< T, Point_t, Tensor_t, GGG_t > Class Template Reference

SphericalTensor that evaluates the Real Spherical Harmonics. More...

+ Collaboration diagram for SphericalTensor< T, Point_t, Tensor_t, GGG_t >:

Public Types

using value_type = T
 
using pos_type = Point_t
 
using hess_type = Tensor_t
 
using ggg_type = GGG_t
 
using This_t = SphericalTensor< T, Point_t >
 

Public Member Functions

 SphericalTensor (const int l_max, bool addsign=false)
 constructor More...
 
void evaluate (const Point_t &p)
 makes a table of $ r^l S_l^m $ and their gradients up to Lmax. More...
 
void evaluateAll (const Point_t &p)
 makes a table of $ r^l S_l^m $ and their gradients up to Lmax. More...
 
void evaluateTest (const Point_t &p)
 makes a table of $ r^l S_l^m $ and their gradients up to Lmax. More...
 
void evaluateWithHessian (const Point_t &p)
 makes a table of $ r^l S_l^m $ and their gradients and hessians up to Lmax. More...
 
void evaluateThirdDerivOnly (const Point_t &p)
 makes a table of $ r^l S_l^m $ and their gradients and hessians and third derivatives up to Lmax. More...
 
int index (int l, int m) const
 returns the index/locator for ( $l,m$) combo, $ l(l+1)+m $ More...
 
value_type getYlm (int l, int m) const
 returns the value of $ r^l S_l^m $ given $ (l,m) $ More...
 
Point_t getGradYlm (int l, int m) const
 returns the gradient of $ r^l S_l^m $ given $ (l,m) $ More...
 
Tensor_t getHessYlm (int l, int m) const
 returns the hessian of $ r^l S_l^m $ given $ (l,m) $ More...
 
GGG_t getGGGYlm (int lm) const
 returns the matrix of third derivatives of $ r^l S_l^m $ given $ (l,m) $ More...
 
value_type getYlm (int lm) const
 returns the value of $ r^l S_l^m $ given $ (l,m) $ More...
 
Point_t getGradYlm (int lm) const
 returns the gradient of $ r^l S_l^m $ given $ (l,m) $ More...
 
Tensor_t getHessYlm (int lm) const
 returns the hessian of $ r^l S_l^m $ given $ (l,m) $ More...
 
int size () const
 
int lmax () const
 

Public Attributes

int Lmax
 maximum angular momentum for the center More...
 
std::vector< value_type > Ylm
 values Ylm $=r^l S_l^m(x,y,z)$ More...
 
std::vector< value_type > NormFactor
 Normalization factors. More...
 
std::vector< value_type > FactorLM
 pre-evaluated factor $1/\sqrt{(l+m)\times(l+1-m)}$ More...
 
std::vector< value_type > FactorL
 pre-evaluated factor $\sqrt{(2l+1)/(4\pi)}$ More...
 
std::vector< value_type > Factor2L
 pre-evaluated factor $(2l+1)/(2l-1)$ More...
 
std::vector< Point_t > gradYlm
 gradients gradYlm $=\nabla r^l S_l^m(x,y,z)$ More...
 
std::vector< hess_type > hessYlm
 hessian hessYlm $=(\nabla X \nabla) r^l S_l^m(x,y,z)$ More...
 
std::vector< value_type > laplYlm
 mmorales: HACK HACK HACK, to avoid having to rewrite QMCWaveFunctions/SphericalBasisSet.h More...
 
std::vector< ggg_type > gggYlm
 

Detailed Description

template<class T, class Point_t, class Tensor_t = qmcplusplus::Tensor<T, 3>, class GGG_t = qmcplusplus::TinyVector<Tensor_t, 3>>
class SphericalTensor< T, Point_t, Tensor_t, GGG_t >

SphericalTensor that evaluates the Real Spherical Harmonics.

The template parameters

  • T, the value_type, e.g. double
  • Point_t, a vector type to provide xyz coordinate. Point_t must have the operator[] defined, e.g., TinyVector<double,3>.

Real Spherical Harmonics Ylm $=r^l S_l^m(x,y,z) $ is stored in an array ordered as [0,-1 0 1,-2 -1 0 1 2, -Lmax,-Lmax+1,..., Lmax-1,Lmax] where Lmax is the maximum angular momentum of a center. All the data members, e.g, Ylm and pre-calculated factors, can be accessed by index(l,m) which returns the locator of the combination for l and m.

Definition at line 40 of file SphericalTensor.h.

Member Typedef Documentation

◆ ggg_type

using ggg_type = GGG_t

Definition at line 46 of file SphericalTensor.h.

◆ hess_type

using hess_type = Tensor_t

Definition at line 45 of file SphericalTensor.h.

◆ pos_type

using pos_type = Point_t

Definition at line 44 of file SphericalTensor.h.

◆ This_t

using This_t = SphericalTensor<T, Point_t>

Definition at line 47 of file SphericalTensor.h.

◆ value_type

using value_type = T

Definition at line 43 of file SphericalTensor.h.

Constructor & Destructor Documentation

◆ SphericalTensor()

SphericalTensor ( const int  l_max,
bool  addsign = false 
)
explicit

constructor

Parameters
l_maxmaximum angular momentum
addsignflag to determine what convention to use

Evaluate all the constants and prefactors. The spherical harmonics is defined as

\[ Y_l^m (\theta,\phi) = \sqrt{\frac{(2l+1)(l-m)!}{4\pi(l+m)!}} P_l^m(\cos\theta)e^{im\phi}\]

Note that the data member Ylm is a misnomer and should not be confused with "spherical harmonics" $Y_l^m$.

  • When addsign == true, e.g., Gaussian packages

    \begin{eqnarray*} S_l^m &=& (-1)^m \sqrt{2}\Re(Y_l^{|m|}), \;\;\;m > 0 \\ &=& Y_l^0, \;\;\;m = 0 \\ &=& (-1)^m \sqrt{2}\Im(Y_l^{|m|}),\;\;\;m < 0 \end{eqnarray*}

  • When addsign == false, e.g., SIESTA package,

    \begin{eqnarray*} S_l^m &=& \sqrt{2}\Re(Y_l^{|m|}), \;\;\;m > 0 \\ &=& Y_l^0, \;\;\;m = 0 \\ &=&\sqrt{2}\Im(Y_l^{|m|}),\;\;\;m < 0 \end{eqnarray*}

Definition at line 141 of file SphericalTensor.h.

References qmcplusplus::atan(), SphericalTensor< T, Point_t, Tensor_t, GGG_t >::Factor2L, SphericalTensor< T, Point_t, Tensor_t, GGG_t >::FactorL, SphericalTensor< T, Point_t, Tensor_t, GGG_t >::FactorLM, SphericalTensor< T, Point_t, Tensor_t, GGG_t >::gggYlm, SphericalTensor< T, Point_t, Tensor_t, GGG_t >::gradYlm, SphericalTensor< T, Point_t, Tensor_t, GGG_t >::hessYlm, SphericalTensor< T, Point_t, Tensor_t, GGG_t >::index(), SphericalTensor< T, Point_t, Tensor_t, GGG_t >::laplYlm, SphericalTensor< T, Point_t, Tensor_t, GGG_t >::Lmax, qmcplusplus::Units::distance::m, SphericalTensor< T, Point_t, Tensor_t, GGG_t >::NormFactor, pi, qmcplusplus::pow(), qmcplusplus::sqrt(), and SphericalTensor< T, Point_t, Tensor_t, GGG_t >::Ylm.

141  : Lmax(l_max)
142 {
143  int ntot = (Lmax + 1) * (Lmax + 1);
144  const value_type pi = 4.0 * atan(1.0);
145  Ylm.resize(ntot);
146  gradYlm.resize(ntot);
147  gradYlm[0] = 0.0;
148  laplYlm.resize(ntot);
149  for (int i = 0; i < ntot; i++)
150  laplYlm[i] = 0.0;
151  hessYlm.resize(ntot);
152  gggYlm.resize(ntot);
153  hessYlm[0] = 0.0;
154  if (ntot >= 4)
155  {
156  hessYlm[1] = 0.0;
157  hessYlm[2] = 0.0;
158  hessYlm[3] = 0.0;
159  }
160  if (ntot >= 9)
161  {
162  // m=-2
163  hessYlm[4] = 0.0;
164  hessYlm[4](1, 0) = std::sqrt(15.0 / 4.0 / pi);
165  hessYlm[4](0, 1) = hessYlm[4](1, 0);
166  // m=-1
167  hessYlm[5] = 0.0;
168  hessYlm[5](1, 2) = std::sqrt(15.0 / 4.0 / pi);
169  hessYlm[5](2, 1) = hessYlm[5](1, 2);
170  // m=0
171  hessYlm[6] = 0.0;
172  hessYlm[6](0, 0) = -std::sqrt(5.0 / 4.0 / pi);
173  hessYlm[6](1, 1) = hessYlm[6](0, 0);
174  hessYlm[6](2, 2) = -2.0 * hessYlm[6](0, 0);
175  // m=1
176  hessYlm[7] = 0.0;
177  hessYlm[7](0, 2) = std::sqrt(15.0 / 4.0 / pi);
178  hessYlm[7](2, 0) = hessYlm[7](0, 2);
179  // m=2
180  hessYlm[8] = 0.0;
181  hessYlm[8](0, 0) = std::sqrt(15.0 / 4.0 / pi);
182  hessYlm[8](1, 1) = -hessYlm[8](0, 0);
183  }
184  for (int i = 0; i < ntot; i++)
185  {
186  gggYlm[i][0] = 0.0;
187  gggYlm[i][1] = 0.0;
188  gggYlm[i][2] = 0.0;
189  }
190  NormFactor.resize(ntot, 1);
191  const value_type sqrt2 = sqrt(2.0);
192  if (addsign)
193  {
194  for (int l = 0; l <= Lmax; l++)
195  {
196  NormFactor[index(l, 0)] = 1.0;
197  for (int m = 1; m <= l; m++)
198  {
199  NormFactor[index(l, m)] = std::pow(-1.0e0, m) * sqrt2;
200  NormFactor[index(l, -m)] = std::pow(-1.0e0, -m) * sqrt2;
201  }
202  }
203  }
204  else
205  {
206  for (int l = 0; l <= Lmax; l++)
207  {
208  for (int m = 1; m <= l; m++)
209  {
210  NormFactor[index(l, m)] = sqrt2;
211  NormFactor[index(l, -m)] = sqrt2;
212  }
213  }
214  }
215  FactorL.resize(Lmax + 1);
216  const value_type omega = 1.0 / sqrt(16.0 * atan(1.0));
217  for (int l = 1; l <= Lmax; l++)
218  FactorL[l] = sqrt(static_cast<T>(2 * l + 1)) * omega;
219  Factor2L.resize(Lmax + 1);
220  for (int l = 1; l <= Lmax; l++)
221  Factor2L[l] = static_cast<T>(2 * l + 1) / static_cast<T>(2 * l - 1);
222  FactorLM.resize(ntot);
223  for (int l = 1; l <= Lmax; l++)
224  for (int m = 1; m <= l; m++)
225  {
226  T fac2 = 1.0 / sqrt(static_cast<T>((l + m) * (l + 1 - m)));
227  FactorLM[index(l, m)] = fac2;
228  FactorLM[index(l, -m)] = fac2;
229  }
230 }
SphericalTensor::Ylm
std::vector< value_type > Ylm
values Ylm
Definition: SphericalTensor.h:120
SphericalTensor::gggYlm
std::vector< ggg_type > gggYlm
Definition: SphericalTensor.h:137
SphericalTensor::FactorLM
std::vector< value_type > FactorLM
pre-evaluated factor
Definition: SphericalTensor.h:124
pi
const double pi
Definition: Standard.h:56
SphericalTensor::gradYlm
std::vector< Point_t > gradYlm
gradients gradYlm
Definition: SphericalTensor.h:130
SphericalTensor::Lmax
int Lmax
maximum angular momentum for the center
Definition: SphericalTensor.h:117
SphericalTensor::index
int index(int l, int m) const
returns the index/locator for ( ) combo,
Definition: SphericalTensor.h:89
SphericalTensor::FactorL
std::vector< value_type > FactorL
pre-evaluated factor
Definition: SphericalTensor.h:126
qmcplusplus::pow
MakeReturn< BinaryNode< FnPow, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t, typename CreateLeaf< Vector< T2, C2 > >::Leaf_t > >::Expression_t pow(const Vector< T1, C1 > &l, const Vector< T2, C2 > &r)
Definition: OhmmsVectorOperators.h:316
qmcplusplus::atan
MakeReturn< UnaryNode< FnArcTan, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t atan(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:48
SphericalTensor::hessYlm
std::vector< hess_type > hessYlm
hessian hessYlm
Definition: SphericalTensor.h:132
qmcplusplus::sqrt
MakeReturn< UnaryNode< FnSqrt, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t sqrt(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:136
SphericalTensor::Factor2L
std::vector< value_type > Factor2L
pre-evaluated factor
Definition: SphericalTensor.h:128
SphericalTensor::laplYlm
std::vector< value_type > laplYlm
mmorales: HACK HACK HACK, to avoid having to rewrite QMCWaveFunctions/SphericalBasisSet.h
Definition: SphericalTensor.h:135
SphericalTensor::NormFactor
std::vector< value_type > NormFactor
Normalization factors.
Definition: SphericalTensor.h:122

Member Function Documentation

◆ evaluate()

void evaluate ( const Point_t &  p)

makes a table of $ r^l S_l^m $ and their gradients up to Lmax.

Definition at line 233 of file SphericalTensor.h.

References qmcplusplus::atan(), qmcplusplus::Units::distance::m, pi, qmcplusplus::sqrt(), and Ylm().

234 {
235  value_type x = p[0], y = p[1], z = p[2];
236  const value_type pi = 4.0 * atan(1.0);
237  const value_type pi4 = 4.0 * pi;
238  const value_type omega = 1.0 / sqrt(pi4);
239  const value_type sqrt2 = sqrt(2.0);
240  /* Calculate r, cos(theta), sin(theta), cos(phi), sin(phi) from input
241  coordinates. Check here the coordinate singularity at cos(theta) = +-1.
242  This also takes care of r=0 case. */
243  value_type cphi, sphi, ctheta;
244  value_type r2xy = x * x + y * y;
245  value_type r = sqrt(r2xy + z * z);
246  if (r2xy < std::numeric_limits<T>::epsilon())
247  {
248  cphi = 0.0;
249  sphi = 1.0;
250  ctheta = (z < 0) ? -1.0 : 1.0;
251  }
252  else
253  {
254  ctheta = z / r;
255  value_type rxyi = 1.0 / sqrt(r2xy);
256  cphi = x * rxyi;
257  sphi = y * rxyi;
258  }
259  value_type stheta = sqrt(1.0 - ctheta * ctheta);
260  /* Now to calculate the associated legendre functions P_lm from the
261  recursion relation from l=0 to Lmax. Conventions of J.D. Jackson,
262  Classical Electrodynamics are used. */
263  Ylm[0] = 1.0;
264  // calculate P_ll and P_l,l-1
265  value_type fac = 1.0;
266  int j = -1;
267  for (int l = 1; l <= Lmax; l++)
268  {
269  j += 2;
270  fac *= -j * stheta;
271  int ll = index(l, l);
272  int l1 = index(l, l - 1);
273  int l2 = index(l - 1, l - 1);
274  Ylm[ll] = fac;
275  Ylm[l1] = j * ctheta * Ylm[l2];
276  }
277  // Use recurence to get other plm's //
278  for (int m = 0; m < Lmax - 1; m++)
279  {
280  int j = 2 * m + 1;
281  for (int l = m + 2; l <= Lmax; l++)
282  {
283  j += 2;
284  int lm = index(l, m);
285  int l1 = index(l - 1, m);
286  int l2 = index(l - 2, m);
287  Ylm[lm] = (ctheta * j * Ylm[l1] - (l + m - 1) * Ylm[l2]) / (l - m);
288  }
289  }
290  // Now to calculate r^l Y_lm. //
291  value_type sphim, cphim, temp;
292  Ylm[0] = omega; //1.0/sqrt(pi4);
293  value_type rpow = 1.0;
294  for (int l = 1; l <= Lmax; l++)
295  {
296  rpow *= r;
297  //fac = rpow*sqrt(static_cast<T>(2*l+1))*omega;//rpow*sqrt((2*l+1)/pi4);
298  //FactorL[l] = sqrt(2*l+1)/sqrt(4*pi)
299  fac = rpow * FactorL[l];
300  int l0 = index(l, 0);
301  Ylm[l0] *= fac;
302  cphim = 1.0;
303  sphim = 0.0;
304  for (int m = 1; m <= l; m++)
305  {
306  temp = cphim * cphi - sphim * sphi;
307  sphim = sphim * cphi + cphim * sphi;
308  cphim = temp;
309  int lm = index(l, m);
310  fac *= FactorLM[lm];
311  temp = fac * Ylm[lm];
312  Ylm[lm] = temp * cphim;
313  lm = index(l, -m);
314  Ylm[lm] = temp * sphim;
315  }
316  }
317  for (int i = 0; i < Ylm.size(); i++)
318  Ylm[i] *= NormFactor[i];
319 }
SphericalTensor::Ylm
std::vector< value_type > Ylm
values Ylm
Definition: SphericalTensor.h:120
SphericalTensor::FactorLM
std::vector< value_type > FactorLM
pre-evaluated factor
Definition: SphericalTensor.h:124
pi
const double pi
Definition: Standard.h:56
SphericalTensor::Lmax
int Lmax
maximum angular momentum for the center
Definition: SphericalTensor.h:117
SphericalTensor::index
int index(int l, int m) const
returns the index/locator for ( ) combo,
Definition: SphericalTensor.h:89
SphericalTensor::FactorL
std::vector< value_type > FactorL
pre-evaluated factor
Definition: SphericalTensor.h:126
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::sqrt
MakeReturn< UnaryNode< FnSqrt, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t sqrt(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:136
SphericalTensor::NormFactor
std::vector< value_type > NormFactor
Normalization factors.
Definition: SphericalTensor.h:122

◆ evaluateAll()

void evaluateAll ( const Point_t &  p)

makes a table of $ r^l S_l^m $ and their gradients up to Lmax.

Definition at line 322 of file SphericalTensor.h.

References qmcplusplus::abs(), qmcplusplus::atan(), qmcplusplus::Units::distance::m, pi, qmcplusplus::sqrt(), and Ylm().

323 {
324  value_type x = p[0], y = p[1], z = p[2];
325  const value_type pi = 4.0 * atan(1.0);
326  const value_type pi4 = 4.0 * pi;
327  const value_type omega = 1.0 / sqrt(pi4);
328  const value_type sqrt2 = sqrt(2.0);
329  /* Calculate r, cos(theta), sin(theta), cos(phi), sin(phi) from input
330  coordinates. Check here the coordinate singularity at cos(theta) = +-1.
331  This also takes care of r=0 case. */
332  value_type cphi, sphi, ctheta;
333  value_type r2xy = x * x + y * y;
334  value_type r = sqrt(r2xy + z * z);
335  if (r2xy < std::numeric_limits<T>::epsilon())
336  {
337  cphi = 0.0;
338  sphi = 1.0;
339  ctheta = (z < 0) ? -1.0 : 1.0;
340  }
341  else
342  {
343  ctheta = z / r;
344  value_type rxyi = 1.0 / sqrt(r2xy);
345  cphi = x * rxyi;
346  sphi = y * rxyi;
347  }
348  value_type stheta = sqrt(1.0 - ctheta * ctheta);
349  /* Now to calculate the associated legendre functions P_lm from the
350  recursion relation from l=0 to Lmax. Conventions of J.D. Jackson,
351  Classical Electrodynamics are used. */
352  Ylm[0] = 1.0;
353  // calculate P_ll and P_l,l-1
354  value_type fac = 1.0;
355  int j = -1;
356  for (int l = 1; l <= Lmax; l++)
357  {
358  j += 2;
359  fac *= -j * stheta;
360  int ll = index(l, l);
361  int l1 = index(l, l - 1);
362  int l2 = index(l - 1, l - 1);
363  Ylm[ll] = fac;
364  Ylm[l1] = j * ctheta * Ylm[l2];
365  }
366  // Use recurence to get other plm's //
367  for (int m = 0; m < Lmax - 1; m++)
368  {
369  int j = 2 * m + 1;
370  for (int l = m + 2; l <= Lmax; l++)
371  {
372  j += 2;
373  int lm = index(l, m);
374  int l1 = index(l - 1, m);
375  int l2 = index(l - 2, m);
376  Ylm[lm] = (ctheta * j * Ylm[l1] - (l + m - 1) * Ylm[l2]) / (l - m);
377  }
378  }
379  // Now to calculate r^l Y_lm. //
380  value_type sphim, cphim, temp;
381  Ylm[0] = omega; //1.0/sqrt(pi4);
382  value_type rpow = 1.0;
383  for (int l = 1; l <= Lmax; l++)
384  {
385  rpow *= r;
386  //fac = rpow*sqrt(static_cast<T>(2*l+1))*omega;//rpow*sqrt((2*l+1)/pi4);
387  //FactorL[l] = sqrt(2*l+1)/sqrt(4*pi)
388  fac = rpow * FactorL[l];
389  int l0 = index(l, 0);
390  Ylm[l0] *= fac;
391  cphim = 1.0;
392  sphim = 0.0;
393  for (int m = 1; m <= l; m++)
394  {
395  //fac2 = (l+m)*(l+1-m);
396  //fac = fac/sqrt(fac2);
397  temp = cphim * cphi - sphim * sphi;
398  sphim = sphim * cphi + cphim * sphi;
399  cphim = temp;
400  int lm = index(l, m);
401  //fac2,fac use precalculated FactorLM
402  fac *= FactorLM[lm];
403  temp = fac * Ylm[lm];
404  Ylm[lm] = temp * cphim;
405  lm = index(l, -m);
406  Ylm[lm] = temp * sphim;
407  }
408  }
409  // Calculating Gradient now//
410  for (int l = 1; l <= Lmax; l++)
411  {
412  //value_type fac = ((value_type) (2*l+1))/(2*l-1);
413  value_type fac = Factor2L[l];
414  for (int m = -l; m <= l; m++)
415  {
416  int lm = index(l - 1, 0);
417  value_type gx, gy, gz, dpr, dpi, dmr, dmi;
418  int ma = std::abs(m);
419  value_type cp = sqrt(fac * (l - ma - 1) * (l - ma));
420  value_type cm = sqrt(fac * (l + ma - 1) * (l + ma));
421  value_type c0 = sqrt(fac * (l - ma) * (l + ma));
422  gz = (l > ma) ? c0 * Ylm[lm + m] : 0.0;
423  if (l > ma + 1)
424  {
425  dpr = cp * Ylm[lm + ma + 1];
426  dpi = cp * Ylm[lm - ma - 1];
427  }
428  else
429  {
430  dpr = 0.0;
431  dpi = 0.0;
432  }
433  if (l > 1)
434  {
435  switch (ma)
436  {
437  case 0:
438  dmr = -cm * Ylm[lm + 1];
439  dmi = cm * Ylm[lm - 1];
440  break;
441  case 1:
442  dmr = cm * Ylm[lm];
443  dmi = 0.0;
444  break;
445  default:
446  dmr = cm * Ylm[lm + ma - 1];
447  dmi = cm * Ylm[lm - ma + 1];
448  }
449  }
450  else
451  {
452  dmr = cm * Ylm[lm];
453  dmi = 0.0;
454  //dmr = (l==1) ? cm*Ylm[lm]:0.0;
455  //dmi = 0.0;
456  }
457  if (m < 0)
458  {
459  gx = 0.5 * (dpi - dmi);
460  gy = -0.5 * (dpr + dmr);
461  }
462  else
463  {
464  gx = 0.5 * (dpr - dmr);
465  gy = 0.5 * (dpi + dmi);
466  }
467  lm = index(l, m);
468  if (ma)
469  gradYlm[lm] = NormFactor[lm] * Point_t(gx, gy, gz);
470  else
471  gradYlm[lm] = Point_t(gx, gy, gz);
472  }
473  }
474  for (int i = 0; i < Ylm.size(); i++)
475  Ylm[i] *= NormFactor[i];
476  //for (int i=0; i<Ylm.size(); i++) gradYlm[i]*= NormFactor[i];
477 }
SphericalTensor::Ylm
std::vector< value_type > Ylm
values Ylm
Definition: SphericalTensor.h:120
qmcplusplus::abs
MakeReturn< UnaryNode< FnFabs, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t abs(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:88
SphericalTensor::FactorLM
std::vector< value_type > FactorLM
pre-evaluated factor
Definition: SphericalTensor.h:124
pi
const double pi
Definition: Standard.h:56
SphericalTensor::gradYlm
std::vector< Point_t > gradYlm
gradients gradYlm
Definition: SphericalTensor.h:130
SphericalTensor::Lmax
int Lmax
maximum angular momentum for the center
Definition: SphericalTensor.h:117
SphericalTensor::index
int index(int l, int m) const
returns the index/locator for ( ) combo,
Definition: SphericalTensor.h:89
SphericalTensor::FactorL
std::vector< value_type > FactorL
pre-evaluated factor
Definition: SphericalTensor.h:126
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::sqrt
MakeReturn< UnaryNode< FnSqrt, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t sqrt(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:136
SphericalTensor::Factor2L
std::vector< value_type > Factor2L
pre-evaluated factor
Definition: SphericalTensor.h:128
SphericalTensor::NormFactor
std::vector< value_type > NormFactor
Normalization factors.
Definition: SphericalTensor.h:122

◆ evaluateTest()

void evaluateTest ( const Point_t &  p)

makes a table of $ r^l S_l^m $ and their gradients up to Lmax.

Definition at line 713 of file SphericalTensor.h.

References SCTFunctor< SCT, L >::apply(), qmcplusplus::atan(), norm(), pi, qmcplusplus::sqrt(), and Ylm().

714 {
715  const value_type pi = 4.0 * atan(1.0);
716  const value_type norm = 1.0 / sqrt(4.0 * pi);
717  Ylm[0] = 1.0;
718  gradYlm[0] = 0.0;
719  switch (Lmax)
720  {
721  case (0):
722  break;
723  case (1):
724  SCTFunctor<This_t, 1>::apply(Ylm, gradYlm, p);
725  break;
726  case (2):
727  SCTFunctor<This_t, 2>::apply(Ylm, gradYlm, p);
728  break;
729  //case(3): SCTFunctor<This_t,3>::apply(Ylm,gradYlm,p); break;
730  //case(4): SCTFunctor<This_t,4>::apply(Ylm,gradYlm,p); break;
731  //case(5): SCTFunctor<This_t,5>::apply(Ylm,gradYlm,p); break;
732  defaults:
733  std::cerr << "Lmax>2 is not valid." << std::endl;
734  break;
735  }
736  for (int i = 0; i < Ylm.size(); i++)
737  Ylm[i] *= norm;
738  for (int i = 0; i < Ylm.size(); i++)
739  gradYlm[i] *= norm;
740 }
SphericalTensor::Ylm
std::vector< value_type > Ylm
values Ylm
Definition: SphericalTensor.h:120
pi
const double pi
Definition: Standard.h:56
SphericalTensor::gradYlm
std::vector< Point_t > gradYlm
gradients gradYlm
Definition: SphericalTensor.h:130
SphericalTensor::Lmax
int Lmax
maximum angular momentum for the center
Definition: SphericalTensor.h:117
norm
double norm(const zVec &c)
Definition: VectorOps.h:118
SCTFunctor::apply
static void apply(std::vector< value_type > &Ylm, std::vector< pos_type > &gYlm, const pos_type &p)
Definition: SphericalTensor.h:650
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::sqrt
MakeReturn< UnaryNode< FnSqrt, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t sqrt(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:136

◆ evaluateThirdDerivOnly()

void evaluateThirdDerivOnly ( const Point_t &  p)

makes a table of $ r^l S_l^m $ and their gradients and hessians and third derivatives up to Lmax.

Definition at line 641 of file SphericalTensor.h.

642 {}

◆ evaluateWithHessian()

void evaluateWithHessian ( const Point_t &  p)

makes a table of $ r^l S_l^m $ and their gradients and hessians up to Lmax.

Definition at line 481 of file SphericalTensor.h.

References qmcplusplus::abs(), qmcplusplus::atan(), qmcplusplus::Units::distance::m, pi, qmcplusplus::sqrt(), and Ylm().

482 {
483  value_type x = p[0], y = p[1], z = p[2];
484  const value_type pi = 4.0 * atan(1.0);
485  const value_type pi4 = 4.0 * pi;
486  const value_type omega = 1.0 / sqrt(pi4);
487  const value_type sqrt2 = sqrt(2.0);
488  /* Calculate r, cos(theta), sin(theta), cos(phi), sin(phi) from input
489  coordinates. Check here the coordinate singularity at cos(theta) = +-1.
490  This also takes care of r=0 case. */
491  value_type cphi, sphi, ctheta;
492  value_type r2xy = x * x + y * y;
493  value_type r = sqrt(r2xy + z * z);
494  if (r2xy < std::numeric_limits<T>::epsilon())
495  {
496  cphi = 0.0;
497  sphi = 1.0;
498  ctheta = (z < 0) ? -1.0 : 1.0;
499  }
500  else
501  {
502  ctheta = z / r;
503  value_type rxyi = 1.0 / sqrt(r2xy);
504  cphi = x * rxyi;
505  sphi = y * rxyi;
506  }
507  value_type stheta = sqrt(1.0 - ctheta * ctheta);
508  /* Now to calculate the associated legendre functions P_lm from the
509  recursion relation from l=0 to Lmax. Conventions of J.D. Jackson,
510  Classical Electrodynamics are used. */
511  Ylm[0] = 1.0;
512  // calculate P_ll and P_l,l-1
513  value_type fac = 1.0;
514  int j = -1;
515  for (int l = 1; l <= Lmax; l++)
516  {
517  j += 2;
518  fac *= -j * stheta;
519  int ll = index(l, l);
520  int l1 = index(l, l - 1);
521  int l2 = index(l - 1, l - 1);
522  Ylm[ll] = fac;
523  Ylm[l1] = j * ctheta * Ylm[l2];
524  }
525  // Use recurence to get other plm's //
526  for (int m = 0; m < Lmax - 1; m++)
527  {
528  int j = 2 * m + 1;
529  for (int l = m + 2; l <= Lmax; l++)
530  {
531  j += 2;
532  int lm = index(l, m);
533  int l1 = index(l - 1, m);
534  int l2 = index(l - 2, m);
535  Ylm[lm] = (ctheta * j * Ylm[l1] - (l + m - 1) * Ylm[l2]) / (l - m);
536  }
537  }
538  // Now to calculate r^l Y_lm. //
539  value_type sphim, cphim, temp;
540  Ylm[0] = omega; //1.0/sqrt(pi4);
541  value_type rpow = 1.0;
542  for (int l = 1; l <= Lmax; l++)
543  {
544  rpow *= r;
545  //fac = rpow*sqrt(static_cast<T>(2*l+1))*omega;//rpow*sqrt((2*l+1)/pi4);
546  //FactorL[l] = sqrt(2*l+1)/sqrt(4*pi)
547  fac = rpow * FactorL[l];
548  int l0 = index(l, 0);
549  Ylm[l0] *= fac;
550  cphim = 1.0;
551  sphim = 0.0;
552  for (int m = 1; m <= l; m++)
553  {
554  //fac2 = (l+m)*(l+1-m);
555  //fac = fac/sqrt(fac2);
556  temp = cphim * cphi - sphim * sphi;
557  sphim = sphim * cphi + cphim * sphi;
558  cphim = temp;
559  int lm = index(l, m);
560  //fac2,fac use precalculated FactorLM
561  fac *= FactorLM[lm];
562  temp = fac * Ylm[lm];
563  Ylm[lm] = temp * cphim;
564  lm = index(l, -m);
565  Ylm[lm] = temp * sphim;
566  }
567  }
568  // Calculating Gradient now//
569  for (int l = 1; l <= Lmax; l++)
570  {
571  //value_type fac = ((value_type) (2*l+1))/(2*l-1);
572  value_type fac = Factor2L[l];
573  for (int m = -l; m <= l; m++)
574  {
575  int lm = index(l - 1, 0);
576  value_type gx, gy, gz, dpr, dpi, dmr, dmi;
577  int ma = std::abs(m);
578  value_type cp = sqrt(fac * (l - ma - 1) * (l - ma));
579  value_type cm = sqrt(fac * (l + ma - 1) * (l + ma));
580  value_type c0 = sqrt(fac * (l - ma) * (l + ma));
581  gz = (l > ma) ? c0 * Ylm[lm + m] : 0.0;
582  if (l > ma + 1)
583  {
584  dpr = cp * Ylm[lm + ma + 1];
585  dpi = cp * Ylm[lm - ma - 1];
586  }
587  else
588  {
589  dpr = 0.0;
590  dpi = 0.0;
591  }
592  if (l > 1)
593  {
594  switch (ma)
595  {
596  case 0:
597  dmr = -cm * Ylm[lm + 1];
598  dmi = cm * Ylm[lm - 1];
599  break;
600  case 1:
601  dmr = cm * Ylm[lm];
602  dmi = 0.0;
603  break;
604  default:
605  dmr = cm * Ylm[lm + ma - 1];
606  dmi = cm * Ylm[lm - ma + 1];
607  }
608  }
609  else
610  {
611  dmr = cm * Ylm[lm];
612  dmi = 0.0;
613  //dmr = (l==1) ? cm*Ylm[lm]:0.0;
614  //dmi = 0.0;
615  }
616  if (m < 0)
617  {
618  gx = 0.5 * (dpi - dmi);
619  gy = -0.5 * (dpr + dmr);
620  }
621  else
622  {
623  gx = 0.5 * (dpr - dmr);
624  gy = 0.5 * (dpi + dmi);
625  }
626  lm = index(l, m);
627  if (ma)
628  gradYlm[lm] = NormFactor[lm] * Point_t(gx, gy, gz);
629  else
630  gradYlm[lm] = Point_t(gx, gy, gz);
631  }
632  }
633  for (int i = 0; i < Ylm.size(); i++)
634  Ylm[i] *= NormFactor[i];
635  //for (int i=0; i<Ylm.size(); i++) gradYlm[i]*= NormFactor[i];
636 }
SphericalTensor::Ylm
std::vector< value_type > Ylm
values Ylm
Definition: SphericalTensor.h:120
qmcplusplus::abs
MakeReturn< UnaryNode< FnFabs, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t abs(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:88
SphericalTensor::FactorLM
std::vector< value_type > FactorLM
pre-evaluated factor
Definition: SphericalTensor.h:124
pi
const double pi
Definition: Standard.h:56
SphericalTensor::gradYlm
std::vector< Point_t > gradYlm
gradients gradYlm
Definition: SphericalTensor.h:130
SphericalTensor::Lmax
int Lmax
maximum angular momentum for the center
Definition: SphericalTensor.h:117
SphericalTensor::index
int index(int l, int m) const
returns the index/locator for ( ) combo,
Definition: SphericalTensor.h:89
SphericalTensor::FactorL
std::vector< value_type > FactorL
pre-evaluated factor
Definition: SphericalTensor.h:126
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::sqrt
MakeReturn< UnaryNode< FnSqrt, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t sqrt(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:136
SphericalTensor::Factor2L
std::vector< value_type > Factor2L
pre-evaluated factor
Definition: SphericalTensor.h:128
SphericalTensor::NormFactor
std::vector< value_type > NormFactor
Normalization factors.
Definition: SphericalTensor.h:122

◆ getGGGYlm()

GGG_t getGGGYlm ( int  lm) const
inline

returns the matrix of third derivatives of $ r^l S_l^m $ given $ (l,m) $

Definition at line 101 of file SphericalTensor.h.

References SphericalTensor< T, Point_t, Tensor_t, GGG_t >::gggYlm.

101 { return gggYlm[lm]; }
SphericalTensor::gggYlm
std::vector< ggg_type > gggYlm
Definition: SphericalTensor.h:137

◆ getGradYlm() [1/2]

Point_t getGradYlm ( int  l,
int  m 
) const
inline

returns the gradient of $ r^l S_l^m $ given $ (l,m) $

Definition at line 95 of file SphericalTensor.h.

References SphericalTensor< T, Point_t, Tensor_t, GGG_t >::gradYlm, SphericalTensor< T, Point_t, Tensor_t, GGG_t >::index(), and qmcplusplus::Units::distance::m.

95 { return gradYlm[index(l, m)]; }
SphericalTensor::gradYlm
std::vector< Point_t > gradYlm
gradients gradYlm
Definition: SphericalTensor.h:130
SphericalTensor::index
int index(int l, int m) const
returns the index/locator for ( ) combo,
Definition: SphericalTensor.h:89

◆ getGradYlm() [2/2]

Point_t getGradYlm ( int  lm) const
inline

returns the gradient of $ r^l S_l^m $ given $ (l,m) $

Definition at line 107 of file SphericalTensor.h.

References SphericalTensor< T, Point_t, Tensor_t, GGG_t >::gradYlm.

107 { return gradYlm[lm]; }
SphericalTensor::gradYlm
std::vector< Point_t > gradYlm
gradients gradYlm
Definition: SphericalTensor.h:130

◆ getHessYlm() [1/2]

Tensor_t getHessYlm ( int  l,
int  m 
) const
inline

returns the hessian of $ r^l S_l^m $ given $ (l,m) $

Definition at line 98 of file SphericalTensor.h.

References SphericalTensor< T, Point_t, Tensor_t, GGG_t >::hessYlm, SphericalTensor< T, Point_t, Tensor_t, GGG_t >::index(), and qmcplusplus::Units::distance::m.

98 { return hessYlm[index(l, m)]; }
SphericalTensor::index
int index(int l, int m) const
returns the index/locator for ( ) combo,
Definition: SphericalTensor.h:89
SphericalTensor::hessYlm
std::vector< hess_type > hessYlm
hessian hessYlm
Definition: SphericalTensor.h:132

◆ getHessYlm() [2/2]

Tensor_t getHessYlm ( int  lm) const
inline

returns the hessian of $ r^l S_l^m $ given $ (l,m) $

Definition at line 110 of file SphericalTensor.h.

References SphericalTensor< T, Point_t, Tensor_t, GGG_t >::hessYlm.

110 { return hessYlm[lm]; }
SphericalTensor::hessYlm
std::vector< hess_type > hessYlm
hessian hessYlm
Definition: SphericalTensor.h:132

◆ getYlm() [1/2]

value_type getYlm ( int  l,
int  m 
) const
inline

returns the value of $ r^l S_l^m $ given $ (l,m) $

Definition at line 92 of file SphericalTensor.h.

References SphericalTensor< T, Point_t, Tensor_t, GGG_t >::index(), qmcplusplus::Units::distance::m, and SphericalTensor< T, Point_t, Tensor_t, GGG_t >::Ylm.

92 { return Ylm[index(l, m)]; }
SphericalTensor::Ylm
std::vector< value_type > Ylm
values Ylm
Definition: SphericalTensor.h:120
SphericalTensor::index
int index(int l, int m) const
returns the index/locator for ( ) combo,
Definition: SphericalTensor.h:89

◆ getYlm() [2/2]

value_type getYlm ( int  lm) const
inline

returns the value of $ r^l S_l^m $ given $ (l,m) $

Definition at line 104 of file SphericalTensor.h.

References SphericalTensor< T, Point_t, Tensor_t, GGG_t >::Ylm.

104 { return Ylm[lm]; }
SphericalTensor::Ylm
std::vector< value_type > Ylm
values Ylm
Definition: SphericalTensor.h:120

◆ index()

int index ( int  l,
int  m 
) const
inline

returns the index/locator for ( $l,m$) combo, $ l(l+1)+m $

Definition at line 89 of file SphericalTensor.h.

References qmcplusplus::Units::distance::m.

Referenced by SphericalTensor< T, Point_t, Tensor_t, GGG_t >::getGradYlm(), SphericalTensor< T, Point_t, Tensor_t, GGG_t >::getHessYlm(), SphericalTensor< T, Point_t, Tensor_t, GGG_t >::getYlm(), and SphericalTensor< T, Point_t, Tensor_t, GGG_t >::SphericalTensor().

89 { return (l * (l + 1)) + m; }

◆ lmax()

int lmax ( ) const
inline

Definition at line 114 of file SphericalTensor.h.

References SphericalTensor< T, Point_t, Tensor_t, GGG_t >::Lmax.

114 { return Lmax; }
SphericalTensor::Lmax
int Lmax
maximum angular momentum for the center
Definition: SphericalTensor.h:117

◆ size()

int size ( void  ) const
inline

Definition at line 112 of file SphericalTensor.h.

References SphericalTensor< T, Point_t, Tensor_t, GGG_t >::Ylm.

112 { return Ylm.size(); }
SphericalTensor::Ylm
std::vector< value_type > Ylm
values Ylm
Definition: SphericalTensor.h:120

Member Data Documentation

◆ Factor2L

std::vector<value_type> Factor2L

pre-evaluated factor $(2l+1)/(2l-1)$

Definition at line 128 of file SphericalTensor.h.

Referenced by SphericalTensor< T, Point_t, Tensor_t, GGG_t >::SphericalTensor().

◆ FactorL

std::vector<value_type> FactorL

pre-evaluated factor $\sqrt{(2l+1)/(4\pi)}$

Definition at line 126 of file SphericalTensor.h.

Referenced by SphericalTensor< T, Point_t, Tensor_t, GGG_t >::SphericalTensor().

◆ FactorLM

std::vector<value_type> FactorLM

pre-evaluated factor $1/\sqrt{(l+m)\times(l+1-m)}$

Definition at line 124 of file SphericalTensor.h.

Referenced by SphericalTensor< T, Point_t, Tensor_t, GGG_t >::SphericalTensor().

◆ gggYlm

std::vector<ggg_type> gggYlm

Definition at line 137 of file SphericalTensor.h.

Referenced by SphericalTensor< T, Point_t, Tensor_t, GGG_t >::getGGGYlm(), and SphericalTensor< T, Point_t, Tensor_t, GGG_t >::SphericalTensor().

◆ gradYlm

std::vector<Point_t> gradYlm

gradients gradYlm $=\nabla r^l S_l^m(x,y,z)$

Definition at line 130 of file SphericalTensor.h.

Referenced by SphericalTensor< T, Point_t, Tensor_t, GGG_t >::getGradYlm(), and SphericalTensor< T, Point_t, Tensor_t, GGG_t >::SphericalTensor().

◆ hessYlm

std::vector<hess_type> hessYlm

hessian hessYlm $=(\nabla X \nabla) r^l S_l^m(x,y,z)$

Definition at line 132 of file SphericalTensor.h.

Referenced by SphericalTensor< T, Point_t, Tensor_t, GGG_t >::getHessYlm(), and SphericalTensor< T, Point_t, Tensor_t, GGG_t >::SphericalTensor().

◆ laplYlm

std::vector<value_type> laplYlm

mmorales: HACK HACK HACK, to avoid having to rewrite QMCWaveFunctions/SphericalBasisSet.h

Definition at line 135 of file SphericalTensor.h.

Referenced by SphericalTensor< T, Point_t, Tensor_t, GGG_t >::SphericalTensor().

◆ Lmax

int Lmax

maximum angular momentum for the center

Definition at line 117 of file SphericalTensor.h.

Referenced by SphericalTensor< T, Point_t, Tensor_t, GGG_t >::lmax(), and SphericalTensor< T, Point_t, Tensor_t, GGG_t >::SphericalTensor().

◆ NormFactor

std::vector<value_type> NormFactor

Normalization factors.

Definition at line 122 of file SphericalTensor.h.

Referenced by SphericalTensor< T, Point_t, Tensor_t, GGG_t >::SphericalTensor().

◆ Ylm

std::vector<value_type> Ylm

values Ylm $=r^l S_l^m(x,y,z)$

Definition at line 120 of file SphericalTensor.h.

Referenced by SphericalTensor< T, Point_t, Tensor_t, GGG_t >::getYlm(), SphericalTensor< T, Point_t, Tensor_t, GGG_t >::size(), and SphericalTensor< T, Point_t, Tensor_t, GGG_t >::SphericalTensor().

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