QMCPACK
test_ylm.cpp
Go to the documentation of this file.
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: Mark Dewing, markdewing@gmail.com, University of Illinois at Urbana-Champaign
8 //
9 // File created by: Mark Dewing, markdewing@gmail.com, University of Illinois at Urbana-Champaign
10 //////////////////////////////////////////////////////////////////////////////////////
11 
12 
13 #include "catch.hpp"
14 
15 #include "Numerics/Ylm.h"
16 
17 #include <stdio.h>
18 #include <string>
19 
20 using std::string;
21 
22 namespace qmcplusplus
23 {
24 TEST_CASE("Legendre", "[numerics]")
25 {
26  // l=0 should be 1.0 for all values
27  double v = LegendrePll(0, 0.5);
28  CHECK(v == Approx(1.0));
29 
30  // check endpoints
31  for (int i = 0; i < 10; i++)
32  {
33  double vp1 = LegendrePll(0, 2.0);
34  CHECK(vp1 == Approx(1.0));
35  double vm1 = LegendrePll(0, 2.0);
36  CHECK(std::abs(vm1) == Approx(1.0));
37  }
38 }
39 
40 TEST_CASE("Spherical Harmonics", "[numerics]")
41 {
42  // l=0 should be 1.0 for all values
43  using vec_t = TinyVector<double, 3>;
44  vec_t v;
45  v[0] = 1.0;
46  v[1] = 0.0;
47  v[2] = 0.0;
48  std::complex<double> out = Ylm(0, 0, v);
49  CHECK(out.real() == Approx(0.28209479177387814)); // Y00 = 1/2 1/sqrt(pi)
50 }
51 
52 TEST_CASE("Spherical Harmonics Many", "[numerics]")
53 {
54  struct Point
55  {
56  double x;
57  double y;
58  double z;
59  };
60 
61  struct YlmValue
62  {
63  Point p;
64  int l;
65  int m;
66  double y_re;
67  double y_im;
68  };
69 
70 #if 0
71  // Use gen_ylm.py to create this file to test more values
72 #include "ylm.inc"
73 #else
74  // Test a small subset of values
75  const int N = 11;
76  YlmValue Vals[N] = {
77  {{0, 0, 1}, 1, -1, 0, 0},
78  {{0.587785, 0, 0.809017}, 1, -1, 0.203076, 0},
79  {{0.587785, 0, 0.809017}, 1, 0, 0.395288, 0},
80  {{0.951057, 0, 0.309017}, 1, 1, -0.328584, 0},
81  {{0.587785, 0, 0.809017}, 2, -2, 0.133454, 0},
82  {{0.293893, 0.904508, 0.309017}, 2, -1, 0.0701612, -0.215934},
83  {{0.587785, 0, -0.809017}, 2, 0, 0.303888, 0},
84  {{0.293893, 0.904508, -0.309017}, 2, 1, 0.0701612, 0.215934},
85  {{0.293893, 0.904508, 0.309017}, 2, 2, -0.282661, 0.205365},
86  {{-0.475528, -0.345492, -0.809017}, 3, -3, 0.0261823, 0.0805808},
87  {{-0.475528, -0.345492, 0.809017}, 4, 4, -0.0427344, 0.0310484},
88  };
89 #endif
90 
91  using vec_t = TinyVector<double, 3>;
92  for (int i = 0; i < N; i++)
93  {
94  YlmValue& v = Vals[i];
95  vec_t w;
96  w[0] = v.p.z; // first component appears to be aligned along the z-axis
97  w[1] = v.p.x;
98  w[2] = v.p.y;
99 
100  std::complex<double> out = Ylm(v.l, v.m, w);
101  //printf("%d %d expected %g %g actual %g %g\n",v.l,v.m,v.y_re,v.y_im, out.real(), out.imag());
102  CHECK(v.y_re == Approx(out.real()));
103  CHECK(v.y_im == Approx(out.imag()));
104  }
105 }
106 
107 TEST_CASE("Derivatives of Spherical Harmonics", "[numerics]")
108 {
109  struct Point
110  {
111  double x;
112  double y;
113  double z;
114  };
115 
116  struct YlmDerivValue
117  {
118  Point p;
119  int l;
120  int m;
121  double th_re;
122  double th_im;
123  double ph_re;
124  double ph_im;
125  };
126 
127  //reference values from sympy
128  //d/dtheta Ylm and d/dphi Ylm
129  const int N = 10;
130  YlmDerivValue Vals[N] = {
131  {{0.587785, 0, 0.809017}, 1, -1, 0.27951007, 0.0, 0.0, -0.2030763},
132  {{0.587785, 0, 0.809017}, 1, 0, -0.2871933, 0.0, 0.0, 0.0},
133  {{0.951057, 0, 0.309017}, 1, 1, -0.1067635, 0.0, 0.0, -0.3285845},
134  {{0.587785, 0, 0.809017}, 2, -2, 0.3673685, 0.0, 0.0, -0.2669088},
135  {{0.293893, 0.904508, 0.309017}, 2, -1, -0.1931373, 0.5944147, -0.2159338, -0.0701613},
136  {{0.587785, 0, -0.809017}, 2, 0, 0.8998655, 0.0, 0.0, 0.0},
137  {{0.293893, 0.904508, -0.309017}, 2, 1, 0.1931373, 0.5944146, -0.2159339, 0.0701613},
138  {{0.293893, 0.904508, 0.309017}, 2, 2, -0.1836842, 0.1334547, -0.4107311, -0.5653217},
139  {{-0.475528, -0.345492, -0.809017}, 3, -3, -0.1081114, -0.3327295, 0.2417422, -0.0785476},
140  {{-0.475528, -0.345492, 0.809017}, 4, 4, -0.2352762, 0.1709368, -0.1241928, -0.1709382},
141  };
142 
143  using vec_t = TinyVector<double, 3>;
144  for (int i = 0; i < N; i++)
145  {
146  YlmDerivValue& v = Vals[i];
147  vec_t w;
148  w[0] = v.p.z; // first component appears to be aligned along the z-axis
149  w[1] = v.p.x;
150  w[2] = v.p.y;
151 
152  std::complex<double> theta, phi;
153  derivYlmSpherical(v.l, v.m, w, theta, phi, false);
154  //printf("%d %d expected %g %g %g %g actual %g %g %g %g\n",v.l,v.m,v.th_re,v.th_im, v.ph_re, v.ph_im, theta.real(), theta.imag(), phi.real(), phi.imag());
155  CHECK(std::real(theta) == Approx(v.th_re));
156  CHECK(std::imag(theta) == Approx(v.th_im));
157  CHECK(std::real(phi) == Approx(v.ph_re));
158  CHECK(std::imag(phi) == Approx(v.ph_im));
159  }
160 }
161 
162 TEST_CASE("Spherical Harmonics Wrapper", "[numerics]")
163 {
164  struct Point
165  {
166  double x;
167  double y;
168  double z;
169  };
170 
171  struct YlmValue
172  {
173  Point p;
174  int l;
175  int m;
176  double y_re;
177  double y_im;
178  };
179 
180  // Test a small subset of values, same test as Spherical Harmonics except some are scaled to not be unit vectors
181  const int N = 11;
182  YlmValue Vals[N] = {
183  {{0, 0, 5}, 1, -1, 0, 0},
184  {{2.938925, 0, 4.045085}, 1, -1, 0.203076, 0},
185  {{2.938925, 0, 4.045085}, 1, 0, 0.395288, 0},
186  {{0.951057, 0, 0.309017}, 1, 1, -0.328584, 0},
187  {{0.587785, 0, 0.809017}, 2, -2, 0.133454, 0},
188  {{1.469465, 4.52254, 1.545085}, 2, -1, 0.0701612, -0.215934},
189  {{0.587785, 0, -0.809017}, 2, 0, 0.303888, 0},
190  {{0.293893, 0.904508, -0.309017}, 2, 1, 0.0701612, 0.215934},
191  {{1.469465, 4.52254, 1.545085}, 2, 2, -0.282661, 0.205365},
192  {{-0.475528, -0.345492, -0.809017}, 3, -3, 0.0261823, 0.0805808},
193  {{-2.37764, -1.72746, 4.045085}, 4, 4, -0.0427344, 0.0310484},
194  };
195 
196  using vec_t = TinyVector<double, 3>;
197  for (int i = 0; i < N; i++)
198  {
199  YlmValue& v = Vals[i];
200  vec_t w;
201  w[0] = v.p.x;
202  w[1] = v.p.y;
203  w[2] = v.p.z;
204 
205  std::complex<double> out = sphericalHarmonic(v.l, v.m, w);
206  //printf("%d %d expected %g %g actual %g %g\n",v.l,v.m,v.y_re,v.y_im, out.real(), out.imag());
207  CHECK(v.y_re == Approx(out.real()));
208  CHECK(v.y_im == Approx(out.imag()));
209  }
210 }
211 
212 TEST_CASE("Cartesian derivatives of Spherical Harmonics", "[numerics]")
213 {
214  struct Point
215  {
216  double x;
217  double y;
218  double z;
219  };
220 
221  struct YlmDerivValue
222  {
223  Point p;
224  int l;
225  int m;
226  double dx_re;
227  double dx_im;
228  double dy_re;
229  double dy_im;
230  double dz_re;
231  double dz_im;
232  };
233 
234  //reference values from sympy
235  // (d/dx, d/dy, d/dz)Ylm
236  const int N = 10;
237  YlmDerivValue Vals[N] =
238  {{{0.587785, 0, 0.809017}, 1, -1, 0.2261289, 0.0, 0.0, -0.3454942, -0.1642922, 0.0},
239  {{2.938925, 0, 4.045085}, 1, 0, -0.0464689, 0.0, 0.0, 0.0, 0.0337616, 0.0},
240  {{4.755285, 0, 1.545085}, 1, 1, -0.0065983, 0.0, 0.0, -0.0690987, 0.020307, 0.0},
241  {{0.587785, 0, 0.809017}, 2, -2, 0.2972075, 0.0, 0.0, -0.4540925, -0.2159337, 0.0},
242  {{1.469465, 4.52254, 1.545085}, 2, -1, 0.0394982, 0.0253845, -0.0253846, 0.0303795, 0.0367369, -0.1130644},
243  {{0.587785, 0, -0.809017}, 2, 0, -0.7280067, 0.0, 0.0, 0.0, -0.5289276, 0.0},
244  {{0.293893, 0.904508, -0.309017}, 2, 1, 0.1974909, -0.1269229, -0.1269229, -0.1518973, -0.183685, -0.5653221},
245  {{1.469465, 4.52254, 1.545085}, 2, 2, 0.0786382, 0.1156131, -0.0374877, -0.0288926, 0.0349388, -0.0253846},
246  {{-0.475528, -0.345492, -0.809017}, 3, -3, 0.1709827, -0.2963218, -0.3841394, -0.0501111, 0.0635463, 0.1955735},
247  {{-2.37764, -1.72746, 4.045085}, 4, 4, 0.0059594, -0.0565636, 0.0565635, 0.0307981, 0.0276584, -0.0200945}};
248 
249  using vec_t = TinyVector<double, 3>;
250  for (int i = 0; i < N; i++)
251  {
252  YlmDerivValue& v = Vals[i];
253  vec_t w;
254  w[0] = v.p.x; // first component appears to be aligned along the z-axis
255  w[1] = v.p.y;
256  w[2] = v.p.z;
257 
258  TinyVector<std::complex<double>, 3> grad;
259  sphericalHarmonicGrad(v.l, v.m, w, grad);
260  //printf("%d %d expected %g %g %g %g actual %g %g %g %g\n",v.l,v.m,v.th_re,v.th_im, v.ph_re, v.ph_im, theta.real(), theta.imag(), phi.real(), phi.imag());
261  CHECK(std::real(grad[0]) == Approx(v.dx_re));
262  CHECK(std::imag(grad[0]) == Approx(v.dx_im));
263  CHECK(std::real(grad[1]) == Approx(v.dy_re));
264  CHECK(std::imag(grad[1]) == Approx(v.dy_im));
265  CHECK(std::real(grad[2]) == Approx(v.dz_re));
266  CHECK(std::imag(grad[2]) == Approx(v.dz_im));
267  }
268 }
269 
270 } // namespace qmcplusplus
qmcplusplus::Units::real
QMCTraits::RealType real
Definition: unit_conversion.h:23
qmcplusplus
helper functions for EinsplineSetBuilder
Definition: Configuration.h:43
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::TEST_CASE
TEST_CASE("complex_helper", "[type_traits]")
Definition: test_complex_helper.cpp:38
qmcplusplus::sphericalHarmonicGrad
void sphericalHarmonicGrad(const int l, const int m, const TinyVector< T, 3 > &r, TinyVector< std::complex< T >, 3 > &grad)
get cartesian derivatives of spherical Harmonics.
Definition: Ylm.h:158
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::derivYlmSpherical
void derivYlmSpherical(const int l, const int m, const TinyVector< T, 3 > &r, std::complex< T > &theta_deriv, std::complex< T > &phi_deriv, const bool conj)
calculate the derivative of a Ylm with respect to theta and phi param[in] l: angular momentum param[i...
Definition: Ylm.h:115
qmcplusplus::CHECK
CHECK(log_values[0]==ComplexApprox(std::complex< double >{ 5.603777579195571, -6.1586603331188225 }))
qmcplusplus::LegendrePll
T LegendrePll(int l, T x)
Definition: Ylm.h:25
qmcplusplus::sphericalHarmonic
std::complex< T > sphericalHarmonic(const int l, const int m, const TinyVector< T, 3 > &r)
wrapper for Ylm, which can take any normal position vector as an argument param[in] l angular momentu...
Definition: Ylm.h:141