1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
8 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
9 //
10 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
11 //////////////////////////////////////////////////////////////////////////////////////
12
13
14 #ifndef QMCPLUSPLUS_YLM_H
15 #define QMCPLUSPLUS_YLM_H
16
17 #include <cmath>
18 #include <complex>
19 #include "OhmmsPETE/TinyVector.h"
20 #include "config/stdlib/Constants.h"
21
22 namespace qmcplusplus
23 {
24 template<typename T>
LegendrePll(int l,T x)25 inline T LegendrePll(int l, T x)
26 {
27 if (l == 0)
28 return 1.0;
29 else
30 {
31 T sqt = std::sqrt(1.0 - x) * std::sqrt(1.0 + x);
32 T val = 1.0;
33 T dblfact = 1.0;
34 for (int i = 1; i <= l; i++)
35 {
36 val *= -sqt;
37 val *= dblfact;
38 dblfact += 2.0;
39 }
40 return val;
41 }
42 }
43
44 template<typename T>
LegendrePlm(int l,int m,T x)45 inline T LegendrePlm(int l, int m, T x)
46 {
47 if (m < 0)
48 {
49 m = std::abs(m);
50 T posval = LegendrePlm(l, m, x);
51 T sign = (m % 2 == 0) ? 1.0 : -1.0;
52 T mfact = 1.0;
53 for (int i = 2; i <= (l - m); i++)
54 mfact *= static_cast<T>(i);
55 T pfact = 1.0;
56 for (int i = 2; i <= (l + m); i++)
57 pfact *= static_cast<T>(i);
58 return posval * sign * mfact / pfact;
59 }
60 // Now we can assume that m is not negative
61 T pmm = LegendrePll(m, x);
62 T pmp1m = x * (2 * m + 1) * pmm;
63 if (m == l)
64 return pmm;
65 else if (l == (m + 1))
66 return pmp1m;
67 else
68 // Use recursive formula
69 {
70 T Plm2m = pmm;
71 T Plm1m = pmp1m;
72 T Pl = 0;
73 for (int i = m + 2; i <= l; i++)
74 {
75 Pl = (1.0 / static_cast<T>(i - m)) * (x * (2 * i - 1) * Plm1m - (i + m - 1) * Plm2m);
76 Plm2m = Plm1m;
77 Plm1m = Pl;
78 }
79 return Pl;
80 }
81 }
82
83 template<typename T>
Ylm(int l,int m,const TinyVector<T,3> & r)84 inline std::complex<T> Ylm(int l, int m, const TinyVector<T, 3>& r)
85 {
86 T costheta = r[0];
87 T phi = std::atan2(r[2], r[1]);
88 int lmm = l - m;
89 int lpm = l + m;
90 T mfact = 1.0;
91 T pfact = 1.0;
92 for (int i = lmm; i > 0; i--)
93 mfact *= static_cast<T>(i);
94 for (int i = lpm; i > 0; i--)
95 pfact *= static_cast<T>(i);
96 T prefactor = std::sqrt(static_cast<T>(2 * l + 1) * mfact / (4.0 * M_PI * pfact));
97 prefactor *= LegendrePlm(l, m, costheta);
98 return std::complex<T>(prefactor * std::cos(m * phi), prefactor * std::sin(m * phi));
99 }
100 } // namespace qmcplusplus
101 #endif
102