QMCPACK
test_spline_applyrotation.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: Joshua Townsend, jptowns@sandia.gov, Sandia National Laboratories
8 //
9 // File created by: Joshua Townsend, jptowns@sandia.gov, Sandia National Laboratories
10 //////////////////////////////////////////////////////////////////////////////////////
11 
12 #include "catch.hpp"
13 
14 #include "OhmmsData/Libxml2Doc.h"
15 #include "OhmmsPETE/OhmmsMatrix.h"
16 #include "Particle/ParticleSet.h"
17 #include "Particle/ParticleSetPool.h"
18 #include "QMCWaveFunctions/WaveFunctionComponent.h"
19 #include "BsplineFactory/EinsplineSetBuilder.h"
20 #include "BsplineFactory/EinsplineSpinorSetBuilder.h"
21 #include "QMCWaveFunctions/BsplineFactory/SplineC2C.h"
22 #include "Utilities/for_testing/checkMatrix.hpp"
23 
24 #include <stdio.h>
25 #include <string>
26 #include <limits>
27 
28 using std::string;
29 
30 namespace qmcplusplus
31 {
32 TEST_CASE("Spline applyRotation zero rotation", "[wavefunction]")
33 {
34  // How to get rid of all this annoying boilerplate?
35  Communicate* c = OHMMS::Controller;
36 
37  // diamondC_1x1x1
38  ParticleSet::ParticleLayout lattice;
39  lattice.R = {3.37316115, 3.37316115, 0.0, 0.0, 3.37316115, 3.37316115, 3.37316115, 0.0, 3.37316115};
40 
41  ParticleSetPool ptcl = ParticleSetPool(c);
42  ptcl.setSimulationCell(lattice);
43  auto ions_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
44  auto elec_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
45  ParticleSet& ions_(*ions_uptr);
46  ParticleSet& elec_(*elec_uptr);
47 
48  ions_.setName("ion");
49  ptcl.addParticleSet(std::move(ions_uptr));
50  ions_.create({2});
51  ions_.R[0] = {0.0, 0.0, 0.0};
52  ions_.R[1] = {1.68658058, 1.68658058, 1.68658058};
53 
54  elec_.setName("elec");
55  ptcl.addParticleSet(std::move(elec_uptr));
56  elec_.create({2});
57  elec_.R[0] = {0.0, 0.0, 0.0};
58  elec_.R[1] = {0.0, 1.0, 0.0};
59 
60  SpeciesSet& tspecies = elec_.getSpeciesSet();
61  int upIdx = tspecies.addSpecies("u");
62  int chargeIdx = tspecies.addAttribute("charge");
63  tspecies(chargeIdx, upIdx) = -1;
64 
65  // Load diamondC_1x1x1 wfn and explicitly construct a SplineC2C object with 7 orbitals
66  // This results in padding of the spline coefs table and thus is a more stringent test.
67  const char* particles = R"(<tmp>
68 <determinantset type="einspline" href="diamondC_1x1x1.pwscf.h5" tilematrix="1 0 0 0 1 0 0 0 1" twistnum="0" source="ion" gpu="no" precision="double" size="7"/>
69 </tmp>)";
70 
71  Libxml2Document doc;
72  bool okay = doc.parseFromString(particles);
73  REQUIRE(okay);
74  xmlNodePtr root = doc.getRoot();
75  xmlNodePtr ein1 = xmlFirstElementChild(root);
76  EinsplineSetBuilder einSet(elec_, ptcl.getPool(), c, ein1);
77  auto spo = einSet.createSPOSetFromXML(ein1);
78  REQUIRE(spo);
79 
80  const auto orbitalsetsize = spo->getOrbitalSetSize();
81  REQUIRE(orbitalsetsize == 7);
82  SPOSet::ValueMatrix psiM_bare(elec_.R.size(), orbitalsetsize);
83  SPOSet::GradMatrix dpsiM_bare(elec_.R.size(), orbitalsetsize);
84  SPOSet::ValueMatrix d2psiM_bare(elec_.R.size(), orbitalsetsize);
85  spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM_bare, dpsiM_bare, d2psiM_bare);
86 
87  // Check before rotation is applied. From 'test_einset.cpp'
88  CHECK(std::real(psiM_bare[1][0]) == Approx(-0.8886948824));
89  CHECK(std::real(psiM_bare[1][1]) == Approx(1.4194120169));
90  // grad
91  CHECK(std::real(dpsiM_bare[1][0][0]) == Approx(-0.0000183403));
92  CHECK(std::real(dpsiM_bare[1][0][1]) == Approx(0.1655139178));
93  CHECK(std::real(dpsiM_bare[1][0][2]) == Approx(-0.0000193077));
94  CHECK(std::real(dpsiM_bare[1][1][0]) == Approx(-1.3131694794));
95  CHECK(std::real(dpsiM_bare[1][1][1]) == Approx(-1.1174004078));
96  CHECK(std::real(dpsiM_bare[1][1][2]) == Approx(-0.8462534547));
97  // lapl
98  CHECK(std::real(d2psiM_bare[1][0]) == Approx(1.3313053846));
99  CHECK(std::real(d2psiM_bare[1][1]) == Approx(-4.712583065));
100 
101  // zero rotation = identity matrix
102  SPOSet::ValueMatrix rot_mat(orbitalsetsize, orbitalsetsize);
103  rot_mat = 0;
104  for (int i = 0; i < orbitalsetsize; i++)
105  rot_mat[i][i] = 1;
106 
107  spo->storeParamsBeforeRotation(); // avoids coefs_copy_ nullptr segfault
108  spo->applyRotation(rot_mat, false);
109 
110  // Get the data for rotated orbitals
111  SPOSet::ValueMatrix psiM_rot(elec_.R.size(), orbitalsetsize);
112  SPOSet::GradMatrix dpsiM_rot(elec_.R.size(), orbitalsetsize);
113  SPOSet::ValueMatrix d2psiM_rot(elec_.R.size(), orbitalsetsize);
114  spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM_rot, dpsiM_rot, d2psiM_rot);
115 
116  // Compute the expected value, gradient, and laplacian by hand using the transformation above
117  // Check value
118  SPOSet::ValueMatrix psiM_rot_manual(elec_.R.size(), orbitalsetsize);
119  for (int i = 0; i < elec_.R.size(); i++)
120  for (int j = 0; j < orbitalsetsize; j++)
121  {
122  psiM_rot_manual[i][j] = 0.;
123  for (int k = 0; k < orbitalsetsize; k++)
124  psiM_rot_manual[i][j] += psiM_bare[i][k] * rot_mat[k][j];
125  }
126  auto check = checkMatrix(psiM_rot_manual, psiM_rot, true);
127  CHECKED_ELSE(check.result) { FAIL(check.result_message); }
128 
129  // Check grad
130  SPOSet::GradMatrix dpsiM_rot_manual(elec_.R.size(), orbitalsetsize);
131  for (int i = 0; i < elec_.R.size(); i++)
132  for (int j = 0; j < orbitalsetsize; j++)
133  {
134  dpsiM_rot_manual[i][j][0] = 0.;
135  dpsiM_rot_manual[i][j][1] = 0.;
136  dpsiM_rot_manual[i][j][2] = 0.;
137  for (int k = 0; k < orbitalsetsize; k++)
138  for (int l = 0; l < 3; l++)
139  dpsiM_rot_manual[i][j][l] += dpsiM_bare[i][k][l] * rot_mat[k][j];
140  }
141 
142  // No checkMatrix for tensors? Gotta do it by hand...
143  double res = 0.;
144  for (int i = 0; i < elec_.R.size(); i++)
145  for (int j = 0; j < orbitalsetsize; j++)
146  for (int k = 0; k < 3; k++)
147  res += std::sqrt(std::norm(dpsiM_rot_manual[i][j][k] - dpsiM_rot[i][j][k]));
148 
149  CHECK(res == Approx(0.).epsilon(2e-4)); // increase threshold to accomodate mixed precision.
150 
151  // Check laplacian
152  SPOSet::ValueMatrix d2psiM_rot_manual(elec_.R.size(), orbitalsetsize);
153  for (int i = 0; i < elec_.R.size(); i++)
154  for (int j = 0; j < orbitalsetsize; j++)
155  {
156  d2psiM_rot_manual[i][j] = 0.;
157  for (int k = 0; k < orbitalsetsize; k++)
158  d2psiM_rot_manual[i][j] += d2psiM_bare[i][k] * rot_mat[k][j];
159  }
160 
161  check = checkMatrix(d2psiM_rot_manual, d2psiM_rot, true, 2e-4);
162  CHECKED_ELSE(check.result) { FAIL(check.result_message); }
163 
164 } // TEST_CASE
165 
166 
167 TEST_CASE("Spline applyRotation one rotation", "[wavefunction]")
168 {
169  // How to get rid of all this annoying boilerplate?
170  Communicate* c = OHMMS::Controller;
171 
172  // diamondC_1x1x1
173  ParticleSet::ParticleLayout lattice;
174  lattice.R = {3.37316115, 3.37316115, 0.0, 0.0, 3.37316115, 3.37316115, 3.37316115, 0.0, 3.37316115};
175 
176  ParticleSetPool ptcl = ParticleSetPool(c);
177  ptcl.setSimulationCell(lattice);
178  auto ions_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
179  auto elec_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
180  ParticleSet& ions_(*ions_uptr);
181  ParticleSet& elec_(*elec_uptr);
182 
183  ions_.setName("ion");
184  ptcl.addParticleSet(std::move(ions_uptr));
185  ions_.create({2});
186  ions_.R[0] = {0.0, 0.0, 0.0};
187  ions_.R[1] = {1.68658058, 1.68658058, 1.68658058};
188 
189  elec_.setName("elec");
190  ptcl.addParticleSet(std::move(elec_uptr));
191  elec_.create({2});
192  elec_.R[0] = {0.0, 0.0, 0.0};
193  elec_.R[1] = {0.0, 1.0, 0.0};
194 
195  SpeciesSet& tspecies = elec_.getSpeciesSet();
196  int upIdx = tspecies.addSpecies("u");
197  int chargeIdx = tspecies.addAttribute("charge");
198  tspecies(chargeIdx, upIdx) = -1;
199 
200  // Load diamondC_1x1x1 wfn and explicitly construct a SplineC2C object with 7 orbitals
201  // This results in padding of the spline coefs table and thus is a more stringent test.
202  const char* particles = R"(<tmp>
203 <determinantset type="einspline" href="diamondC_1x1x1.pwscf.h5" tilematrix="1 0 0 0 1 0 0 0 1" twistnum="0" source="ion" gpu="no" precision="double" size="7"/>
204 </tmp>)";
205 
206  Libxml2Document doc;
207  bool okay = doc.parseFromString(particles);
208  REQUIRE(okay);
209  xmlNodePtr root = doc.getRoot();
210  xmlNodePtr ein1 = xmlFirstElementChild(root);
211  EinsplineSetBuilder einSet(elec_, ptcl.getPool(), c, ein1);
212  auto spo = einSet.createSPOSetFromXML(ein1);
213  REQUIRE(spo);
214 
215  const auto orbitalsetsize = spo->getOrbitalSetSize();
216  REQUIRE(orbitalsetsize == 7);
217  SPOSet::ValueMatrix psiM_bare(elec_.R.size(), orbitalsetsize);
218  SPOSet::GradMatrix dpsiM_bare(elec_.R.size(), orbitalsetsize);
219  SPOSet::ValueMatrix d2psiM_bare(elec_.R.size(), orbitalsetsize);
220  spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM_bare, dpsiM_bare, d2psiM_bare);
221 
222  // Check before rotation is applied. From 'test_einset.cpp'
223  CHECK(std::real(psiM_bare[1][0]) == Approx(-0.8886948824));
224  CHECK(std::real(psiM_bare[1][1]) == Approx(1.4194120169));
225  // grad
226  CHECK(std::real(dpsiM_bare[1][0][0]) == Approx(-0.0000183403));
227  CHECK(std::real(dpsiM_bare[1][0][1]) == Approx(0.1655139178));
228  CHECK(std::real(dpsiM_bare[1][0][2]) == Approx(-0.0000193077));
229  CHECK(std::real(dpsiM_bare[1][1][0]) == Approx(-1.3131694794));
230  CHECK(std::real(dpsiM_bare[1][1][1]) == Approx(-1.1174004078));
231  CHECK(std::real(dpsiM_bare[1][1][2]) == Approx(-0.8462534547));
232  // lapl
233  CHECK(std::real(d2psiM_bare[1][0]) == Approx(1.3313053846));
234  CHECK(std::real(d2psiM_bare[1][1]) == Approx(-4.712583065));
235 
236  /* Apply a rotation, U = exp(-K) as shown below, to the splines.
237  Because K = -K.T, U is unitary.
238  K=
239  0.00000000e+00 -2.33449314e-01 5.21735139e-01 1.67744276e-01 -1.68565994e-01 1.26632312e-02 2.29272040e-01
240  2.33449314e-01 0.00000000e+00 -3.81982621e-01 -6.76995896e-01 7.15727948e-02 9.10926961e-01 -6.94864205e-01
241  -5.21735139e-01 3.81982621e-01 0.00000000e+00 2.37433139e-01 -7.09878105e-01 -6.33841289e-01 -7.49009582e-01
242  -1.67744276e-01 6.76995896e-01 -2.37433139e-01 0.00000000e+00 -6.72002172e-01 4.27152921e-01 -4.82983743e-03
243  1.68565994e-01 -7.15727948e-02 7.09878105e-01 6.72002172e-01 0.00000000e+00 -7.93860012e-01 -7.76624731e-01
244  -1.26632312e-02 -9.10926961e-01 6.33841289e-01 -4.27152921e-01 7.93860012e-01 0.00000000e+00 -2.74931052e-01
245  -2.29272040e-01 6.94864205e-01 7.49009582e-01 4.82983743e-03 7.76624731e-01 2.74931052e-01 0.00000000e+00
246  U=
247  8.06061880e-01 3.54921598e-01 -3.40706426e-01 -5.90163619e-02 1.11650454e-01 -1.99768450e-01 -2.28818375e-01
248  1.58069821e-02 2.10363421e-01 2.74922448e-01 2.12581764e-01 2.64602356e-01 -6.34971914e-01 6.01265560e-01
249  4.43696646e-01 -6.06912539e-02 2.61413193e-01 -1.98368802e-01 -6.43234645e-02 5.75880430e-01 5.96646495e-01
250  1.45865363e-01 -5.16577220e-01 -2.09866367e-01 5.55699395e-01 5.73062990e-01 1.74778224e-01 8.77170506e-03
251  -3.24609748e-01 2.89664179e-01 -7.50613752e-01 -1.63060005e-01 1.70803377e-01 1.63784167e-01 4.05850414e-01
252  1.32570771e-03 3.80299846e-01 -5.08439810e-02 7.59141791e-01 -4.77844928e-01 2.12149087e-01 5.60882349e-02
253  1.62781208e-01 -5.75073150e-01 -3.60485665e-01 -1.70070331e-02 -5.72251258e-01 -3.50549638e-01 2.49394158e-01
254  UU.T=
255  1.00000000e+00 -2.77555756e-17 5.55111512e-17 -8.67361738e-18 -1.11022302e-16 -6.07153217e-17 6.93889390e-17
256  -2.77555756e-17 1.00000000e+00 5.55111512e-17 -1.84748050e-16 5.55111512e-17 -6.93889390e-18 8.32667268e-17
257  5.55111512e-17 5.55111512e-17 1.00000000e+00 1.45716772e-16 -5.55111512e-17 1.38777878e-16 5.55111512e-17
258  -8.67361738e-18 -1.84748050e-16 1.45716772e-16 1.00000000e+00 -1.95156391e-17 1.80411242e-16 8.02309608e-17
259  -1.11022302e-16 5.55111512e-17 -5.55111512e-17 -1.95156391e-17 1.00000000e+00 7.28583860e-17 1.11022302e-16
260  -6.07153217e-17 -6.93889390e-18 1.38777878e-16 1.80411242e-16 7.28583860e-17 1.00000000e+00 -1.12757026e-16
261  6.93889390e-17 8.32667268e-17 5.55111512e-17 8.02309608e-17 1.11022302e-16 -1.12757026e-16 1.00000000e+00
262 
263  NB: There's nothing special about this rotation. I purposefully made something 'ugly' to make a worst case scenario...
264  */
265  SPOSet::ValueMatrix rot_mat(orbitalsetsize, orbitalsetsize);
266  rot_mat[0][0] = 8.06061880e-01;
267  rot_mat[0][1] = 3.54921598e-01;
268  rot_mat[0][2] = -3.40706426e-01;
269  rot_mat[0][3] = -5.90163619e-02;
270  rot_mat[0][4] = 1.11650454e-01;
271  rot_mat[0][5] = -1.99768450e-01;
272  rot_mat[0][6] = -2.28818375e-01;
273  rot_mat[1][0] = 1.58069821e-02;
274  rot_mat[1][1] = 2.10363421e-01;
275  rot_mat[1][2] = 2.74922448e-01;
276  rot_mat[1][3] = 2.12581764e-01;
277  rot_mat[1][4] = 2.64602356e-01;
278  rot_mat[1][5] = -6.34971914e-01;
279  rot_mat[1][6] = 6.01265560e-01;
280  rot_mat[2][0] = 4.43696646e-01;
281  rot_mat[2][1] = -6.06912539e-02;
282  rot_mat[2][2] = 2.61413193e-01;
283  rot_mat[2][3] = -1.98368802e-01;
284  rot_mat[2][4] = -6.43234645e-02;
285  rot_mat[2][5] = 5.75880430e-01;
286  rot_mat[2][6] = 5.96646495e-01;
287  rot_mat[3][0] = 1.45865363e-01;
288  rot_mat[3][1] = -5.16577220e-01;
289  rot_mat[3][2] = -2.09866367e-01;
290  rot_mat[3][3] = 5.55699395e-01;
291  rot_mat[3][4] = 5.73062990e-01;
292  rot_mat[3][5] = 1.74778224e-01;
293  rot_mat[3][6] = 8.77170506e-03;
294  rot_mat[4][0] = -3.24609748e-01;
295  rot_mat[4][1] = 2.89664179e-01;
296  rot_mat[4][2] = -7.50613752e-01;
297  rot_mat[4][3] = -1.63060005e-01;
298  rot_mat[4][4] = 1.70803377e-01;
299  rot_mat[4][5] = 1.63784167e-01;
300  rot_mat[4][6] = 4.05850414e-01;
301  rot_mat[5][0] = 1.32570771e-03;
302  rot_mat[5][1] = 3.80299846e-01;
303  rot_mat[5][2] = -5.08439810e-02;
304  rot_mat[5][3] = 7.59141791e-01;
305  rot_mat[5][4] = -4.77844928e-01;
306  rot_mat[5][5] = 2.12149087e-01;
307  rot_mat[5][6] = 5.60882349e-02;
308  rot_mat[6][0] = 1.62781208e-01;
309  rot_mat[6][1] = -5.75073150e-01;
310  rot_mat[6][2] = -3.60485665e-01;
311  rot_mat[6][3] = -1.70070331e-02;
312  rot_mat[6][4] = -5.72251258e-01;
313  rot_mat[6][5] = -3.50549638e-01;
314  rot_mat[6][6] = 2.49394158e-01;
315  spo->storeParamsBeforeRotation(); // avoids coefs_copy_ nullptr segfault
316  spo->applyRotation(rot_mat, false);
317 
318  // Get the data for rotated orbitals
319  SPOSet::ValueMatrix psiM_rot(elec_.R.size(), orbitalsetsize);
320  SPOSet::GradMatrix dpsiM_rot(elec_.R.size(), orbitalsetsize);
321  SPOSet::ValueMatrix d2psiM_rot(elec_.R.size(), orbitalsetsize);
322  spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM_rot, dpsiM_rot, d2psiM_rot);
323 
324  // Compute the expected value, gradient, and laplacian by hand using the transformation above
325  // Check value
326  SPOSet::ValueMatrix psiM_rot_manual(elec_.R.size(), orbitalsetsize);
327  for (int i = 0; i < elec_.R.size(); i++)
328  for (int j = 0; j < orbitalsetsize; j++)
329  {
330  psiM_rot_manual[i][j] = 0.;
331  for (int k = 0; k < orbitalsetsize; k++)
332  psiM_rot_manual[i][j] += psiM_bare[i][k] * rot_mat[k][j];
333  }
334  auto check = checkMatrix(psiM_rot_manual, psiM_rot, true);
335  CHECKED_ELSE(check.result) { FAIL(check.result_message); }
336 
337  // Check grad
338  SPOSet::GradMatrix dpsiM_rot_manual(elec_.R.size(), orbitalsetsize);
339  for (int i = 0; i < elec_.R.size(); i++)
340  for (int j = 0; j < orbitalsetsize; j++)
341  {
342  dpsiM_rot_manual[i][j][0] = 0.;
343  dpsiM_rot_manual[i][j][1] = 0.;
344  dpsiM_rot_manual[i][j][2] = 0.;
345  for (int k = 0; k < orbitalsetsize; k++)
346  for (int l = 0; l < 3; l++)
347  dpsiM_rot_manual[i][j][l] += dpsiM_bare[i][k][l] * rot_mat[k][j];
348  }
349 
350  // No checkMatrix for tensors? Gotta do it by hand...
351  double res = 0.;
352  for (int i = 0; i < elec_.R.size(); i++)
353  for (int j = 0; j < orbitalsetsize; j++)
354  for (int k = 0; k < 3; k++)
355  res += std::sqrt(std::norm(dpsiM_rot_manual[i][j][k] - dpsiM_rot[i][j][k]));
356 
357  CHECK(res == Approx(0.).epsilon(2e-4)); // increase threshold to accomodate mixed precision.
358 
359  // Check laplacian
360  SPOSet::ValueMatrix d2psiM_rot_manual(elec_.R.size(), orbitalsetsize);
361  for (int i = 0; i < elec_.R.size(); i++)
362  for (int j = 0; j < orbitalsetsize; j++)
363  {
364  d2psiM_rot_manual[i][j] = 0.;
365  for (int k = 0; k < orbitalsetsize; k++)
366  d2psiM_rot_manual[i][j] += d2psiM_bare[i][k] * rot_mat[k][j];
367  }
368 
369  check = checkMatrix(d2psiM_rot_manual, d2psiM_rot, true, 2e-4);
370  CHECKED_ELSE(check.result) { FAIL(check.result_message); }
371 
372 } // TEST_CASE
373 
374 
375 TEST_CASE("Spline applyRotation two rotations", "[wavefunction]")
376 {
377  // How to get rid of all this annoying boilerplate?
378  Communicate* c = OHMMS::Controller;
379 
380  // diamondC_1x1x1
381  ParticleSet::ParticleLayout lattice;
382  lattice.R = {3.37316115, 3.37316115, 0.0, 0.0, 3.37316115, 3.37316115, 3.37316115, 0.0, 3.37316115};
383 
384  ParticleSetPool ptcl = ParticleSetPool(c);
385  ptcl.setSimulationCell(lattice);
386  auto ions_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
387  auto elec_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
388  ParticleSet& ions_(*ions_uptr);
389  ParticleSet& elec_(*elec_uptr);
390 
391  ions_.setName("ion");
392  ptcl.addParticleSet(std::move(ions_uptr));
393  ions_.create({2});
394  ions_.R[0] = {0.0, 0.0, 0.0};
395  ions_.R[1] = {1.68658058, 1.68658058, 1.68658058};
396 
397  elec_.setName("elec");
398  ptcl.addParticleSet(std::move(elec_uptr));
399  elec_.create({2});
400  elec_.R[0] = {0.0, 0.0, 0.0};
401  elec_.R[1] = {0.0, 1.0, 0.0};
402 
403  SpeciesSet& tspecies = elec_.getSpeciesSet();
404  int upIdx = tspecies.addSpecies("u");
405  int chargeIdx = tspecies.addAttribute("charge");
406  tspecies(chargeIdx, upIdx) = -1;
407 
408  // Load diamondC_1x1x1 wfn and explicitly construct a SplineC2C object with 7 orbitals
409  // This results in padding of the spline coefs table and thus is a more stringent test.
410  const char* particles = R"(<tmp>
411 <determinantset type="einspline" href="diamondC_1x1x1.pwscf.h5" tilematrix="1 0 0 0 1 0 0 0 1" twistnum="0" source="ion" gpu="no" precision="double" size="7"/>
412 </tmp>)";
413 
414  Libxml2Document doc;
415  bool okay = doc.parseFromString(particles);
416  REQUIRE(okay);
417  xmlNodePtr root = doc.getRoot();
418  xmlNodePtr ein1 = xmlFirstElementChild(root);
419  EinsplineSetBuilder einSet(elec_, ptcl.getPool(), c, ein1);
420  auto spo = einSet.createSPOSetFromXML(ein1);
421  REQUIRE(spo);
422 
423  const auto orbitalsetsize = spo->getOrbitalSetSize();
424  REQUIRE(orbitalsetsize == 7);
425  SPOSet::ValueMatrix psiM_bare(elec_.R.size(), orbitalsetsize);
426  SPOSet::GradMatrix dpsiM_bare(elec_.R.size(), orbitalsetsize);
427  SPOSet::ValueMatrix d2psiM_bare(elec_.R.size(), orbitalsetsize);
428  spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM_bare, dpsiM_bare, d2psiM_bare);
429 
430  // Check before rotation is applied. From 'test_einset.cpp'
431  CHECK(std::real(psiM_bare[1][0]) == Approx(-0.8886948824));
432  CHECK(std::real(psiM_bare[1][1]) == Approx(1.4194120169));
433  // grad
434  CHECK(std::real(dpsiM_bare[1][0][0]) == Approx(-0.0000183403));
435  CHECK(std::real(dpsiM_bare[1][0][1]) == Approx(0.1655139178));
436  CHECK(std::real(dpsiM_bare[1][0][2]) == Approx(-0.0000193077));
437  CHECK(std::real(dpsiM_bare[1][1][0]) == Approx(-1.3131694794));
438  CHECK(std::real(dpsiM_bare[1][1][1]) == Approx(-1.1174004078));
439  CHECK(std::real(dpsiM_bare[1][1][2]) == Approx(-0.8462534547));
440  // lapl
441  CHECK(std::real(d2psiM_bare[1][0]) == Approx(1.3313053846));
442  CHECK(std::real(d2psiM_bare[1][1]) == Approx(-4.712583065));
443 
444  /* Apply a rotation, U = exp(-K) as shown below, to the splines.
445  Because K = -K.T, U is unitary.
446  K=
447  0.00000000e+00 -2.33449314e-01 5.21735139e-01 1.67744276e-01 -1.68565994e-01 1.26632312e-02 2.29272040e-01
448  2.33449314e-01 0.00000000e+00 -3.81982621e-01 -6.76995896e-01 7.15727948e-02 9.10926961e-01 -6.94864205e-01
449  -5.21735139e-01 3.81982621e-01 0.00000000e+00 2.37433139e-01 -7.09878105e-01 -6.33841289e-01 -7.49009582e-01
450  -1.67744276e-01 6.76995896e-01 -2.37433139e-01 0.00000000e+00 -6.72002172e-01 4.27152921e-01 -4.82983743e-03
451  1.68565994e-01 -7.15727948e-02 7.09878105e-01 6.72002172e-01 0.00000000e+00 -7.93860012e-01 -7.76624731e-01
452  -1.26632312e-02 -9.10926961e-01 6.33841289e-01 -4.27152921e-01 7.93860012e-01 0.00000000e+00 -2.74931052e-01
453  -2.29272040e-01 6.94864205e-01 7.49009582e-01 4.82983743e-03 7.76624731e-01 2.74931052e-01 0.00000000e+00
454  U=
455  8.06061880e-01 3.54921598e-01 -3.40706426e-01 -5.90163619e-02 1.11650454e-01 -1.99768450e-01 -2.28818375e-01
456  1.58069821e-02 2.10363421e-01 2.74922448e-01 2.12581764e-01 2.64602356e-01 -6.34971914e-01 6.01265560e-01
457  4.43696646e-01 -6.06912539e-02 2.61413193e-01 -1.98368802e-01 -6.43234645e-02 5.75880430e-01 5.96646495e-01
458  1.45865363e-01 -5.16577220e-01 -2.09866367e-01 5.55699395e-01 5.73062990e-01 1.74778224e-01 8.77170506e-03
459  -3.24609748e-01 2.89664179e-01 -7.50613752e-01 -1.63060005e-01 1.70803377e-01 1.63784167e-01 4.05850414e-01
460  1.32570771e-03 3.80299846e-01 -5.08439810e-02 7.59141791e-01 -4.77844928e-01 2.12149087e-01 5.60882349e-02
461  1.62781208e-01 -5.75073150e-01 -3.60485665e-01 -1.70070331e-02 -5.72251258e-01 -3.50549638e-01 2.49394158e-01
462  UU.T=
463  1.00000000e+00 -2.77555756e-17 5.55111512e-17 -8.67361738e-18 -1.11022302e-16 -6.07153217e-17 6.93889390e-17
464  -2.77555756e-17 1.00000000e+00 5.55111512e-17 -1.84748050e-16 5.55111512e-17 -6.93889390e-18 8.32667268e-17
465  5.55111512e-17 5.55111512e-17 1.00000000e+00 1.45716772e-16 -5.55111512e-17 1.38777878e-16 5.55111512e-17
466  -8.67361738e-18 -1.84748050e-16 1.45716772e-16 1.00000000e+00 -1.95156391e-17 1.80411242e-16 8.02309608e-17
467  -1.11022302e-16 5.55111512e-17 -5.55111512e-17 -1.95156391e-17 1.00000000e+00 7.28583860e-17 1.11022302e-16
468  -6.07153217e-17 -6.93889390e-18 1.38777878e-16 1.80411242e-16 7.28583860e-17 1.00000000e+00 -1.12757026e-16
469  6.93889390e-17 8.32667268e-17 5.55111512e-17 8.02309608e-17 1.11022302e-16 -1.12757026e-16 1.00000000e+00
470 
471  NB: There's nothing special about this rotation. I purposefully made something 'ugly' to make a worst case scenario...
472  */
473  SPOSet::ValueMatrix rot_mat1(orbitalsetsize, orbitalsetsize);
474  rot_mat1[0][0] = 8.06061880e-01;
475  rot_mat1[0][1] = 3.54921598e-01;
476  rot_mat1[0][2] = -3.40706426e-01;
477  rot_mat1[0][3] = -5.90163619e-02;
478  rot_mat1[0][4] = 1.11650454e-01;
479  rot_mat1[0][5] = -1.99768450e-01;
480  rot_mat1[0][6] = -2.28818375e-01;
481  rot_mat1[1][0] = 1.58069821e-02;
482  rot_mat1[1][1] = 2.10363421e-01;
483  rot_mat1[1][2] = 2.74922448e-01;
484  rot_mat1[1][3] = 2.12581764e-01;
485  rot_mat1[1][4] = 2.64602356e-01;
486  rot_mat1[1][5] = -6.34971914e-01;
487  rot_mat1[1][6] = 6.01265560e-01;
488  rot_mat1[2][0] = 4.43696646e-01;
489  rot_mat1[2][1] = -6.06912539e-02;
490  rot_mat1[2][2] = 2.61413193e-01;
491  rot_mat1[2][3] = -1.98368802e-01;
492  rot_mat1[2][4] = -6.43234645e-02;
493  rot_mat1[2][5] = 5.75880430e-01;
494  rot_mat1[2][6] = 5.96646495e-01;
495  rot_mat1[3][0] = 1.45865363e-01;
496  rot_mat1[3][1] = -5.16577220e-01;
497  rot_mat1[3][2] = -2.09866367e-01;
498  rot_mat1[3][3] = 5.55699395e-01;
499  rot_mat1[3][4] = 5.73062990e-01;
500  rot_mat1[3][5] = 1.74778224e-01;
501  rot_mat1[3][6] = 8.77170506e-03;
502  rot_mat1[4][0] = -3.24609748e-01;
503  rot_mat1[4][1] = 2.89664179e-01;
504  rot_mat1[4][2] = -7.50613752e-01;
505  rot_mat1[4][3] = -1.63060005e-01;
506  rot_mat1[4][4] = 1.70803377e-01;
507  rot_mat1[4][5] = 1.63784167e-01;
508  rot_mat1[4][6] = 4.05850414e-01;
509  rot_mat1[5][0] = 1.32570771e-03;
510  rot_mat1[5][1] = 3.80299846e-01;
511  rot_mat1[5][2] = -5.08439810e-02;
512  rot_mat1[5][3] = 7.59141791e-01;
513  rot_mat1[5][4] = -4.77844928e-01;
514  rot_mat1[5][5] = 2.12149087e-01;
515  rot_mat1[5][6] = 5.60882349e-02;
516  rot_mat1[6][0] = 1.62781208e-01;
517  rot_mat1[6][1] = -5.75073150e-01;
518  rot_mat1[6][2] = -3.60485665e-01;
519  rot_mat1[6][3] = -1.70070331e-02;
520  rot_mat1[6][4] = -5.72251258e-01;
521  rot_mat1[6][5] = -3.50549638e-01;
522  rot_mat1[6][6] = 2.49394158e-01;
523  spo->storeParamsBeforeRotation(); // avoids coefs_copy_ nullptr segfault
524  spo->applyRotation(rot_mat1, false); // Apply rotation 1
525 
526 
527  /* Apply another rotation, U = exp(-K) as shown below, to the splines.
528  Because K = -K.T, U is unitary.
529  K=
530  0.00000000e+00 3.08898211e-01 9.70811450e-01 9.96582440e-01 2.59290113e-01 9.08544511e-01 7.47970513e-01
531  -3.08898211e-01 0.00000000e+00 4.90958419e-01 7.98394113e-01 9.02926177e-01 3.24156332e-01 1.83039207e-01
532  -9.70811450e-01 -4.90958419e-01 0.00000000e+00 7.73447329e-01 2.71156433e-01 2.18009012e-01 1.75304629e-01
533  -9.96582440e-01 -7.98394113e-01 -7.73447329e-01 0.00000000e+00 9.27679516e-01 8.09853231e-01 7.71485500e-01
534  -2.59290113e-01 -9.02926177e-01 -2.71156433e-01 -9.27679516e-01 0.00000000e+00 6.19032776e-01 6.05979744e-02
535  -9.08544511e-01 -3.24156332e-01 -2.18009012e-01 -8.09853231e-01 -6.19032776e-01 0.00000000e+00 7.91708106e-01
536  -7.47970513e-01 -1.83039207e-01 -1.75304629e-01 -7.71485500e-01 -6.05979744e-02 -7.91708106e-01 0.00000000e+00
537 
538  U=
539  -2.98194829e-02 -6.14346084e-01 -6.43923495e-01 -1.18552217e-01 3.77244445e-01 -2.06047353e-01 9.07122365e-02
540  -3.48818370e-01 4.38506143e-01 -5.57140467e-01 -3.70956624e-01 -2.23467853e-01 3.80270838e-01 2.08518583e-01
541  1.87644050e-01 -1.34182635e-01 3.74618322e-01 -8.74132032e-01 1.57504581e-01 4.78555353e-02 1.23455215e-01
542  -7.68247448e-02 -2.48006571e-01 -7.30866440e-02 -2.06817660e-01 -7.99320897e-01 -4.53458397e-01 -1.99842648e-01
543  2.59360916e-02 5.64417889e-01 -1.05874452e-01 -1.09701962e-01 2.76429977e-01 -7.59899812e-01 6.04531910e-02
544  3.50608068e-01 1.58922313e-01 -2.49914906e-01 -1.49783413e-01 1.03865829e-01 1.56180926e-01 -8.55420691e-01
545  8.44230723e-01 8.33975900e-02 -2.35816939e-01 8.35859456e-02 -2.38381031e-01 5.64243083e-02 3.97132056e-01
546 
547  UU.T=
548  1.00000000e+00 1.97376302e-16 2.32170398e-16 8.48209595e-17 1.52095225e-16 1.20549468e-16 -6.20012223e-17
549  1.97376302e-16 1.00000000e+00 2.49491048e-16 -1.32824448e-16 -2.95131454e-17 1.89706297e-17 -9.33948863e-17
550  2.32170398e-16 2.49491048e-16 1.00000000e+00 -1.53719614e-16 -1.07440039e-16 1.11845140e-17 1.25992152e-16
551  8.48209595e-17 -1.32824448e-16 -1.53719614e-16 1.00000000e+00 3.98945291e-17 8.97184517e-17 -1.23231760e-16
552  1.52095225e-16 -2.95131454e-17 -1.07440039e-16 3.98945291e-17 1.00000000e+00 2.40723889e-17 3.26140430e-17
553  1.20549468e-16 1.89706297e-17 1.11845140e-17 8.97184517e-17 2.40723889e-17 1.00000000e+00 2.75131978e-17
554  -6.20012223e-17 -9.33948863e-17 1.25992152e-16 -1.23231760e-16 3.26140430e-17 2.75131978e-17 1.00000000e+00
555 
556  NB: There's nothing special about this rotation. I purposefully made something 'ugly' to make a worst case scenario...
557  */
558  SPOSet::ValueMatrix rot_mat2(orbitalsetsize, orbitalsetsize);
559  rot_mat2[0][0] = -2.98194829e-02;
560  rot_mat2[0][1] = -6.14346084e-01;
561  rot_mat2[0][2] = -6.43923495e-01;
562  rot_mat2[0][3] = -1.18552217e-01;
563  rot_mat2[0][4] = 3.77244445e-01;
564  rot_mat2[0][5] = -2.06047353e-01;
565  rot_mat2[0][6] = 9.07122365e-02;
566  rot_mat2[1][0] = -3.48818370e-01;
567  rot_mat2[1][1] = 4.38506143e-01;
568  rot_mat2[1][2] = -5.57140467e-01;
569  rot_mat2[1][3] = -3.70956624e-01;
570  rot_mat2[1][4] = -2.23467853e-01;
571  rot_mat2[1][5] = 3.80270838e-01;
572  rot_mat2[1][6] = 2.08518583e-01;
573  rot_mat2[2][0] = 1.87644050e-01;
574  rot_mat2[2][1] = -1.34182635e-01;
575  rot_mat2[2][2] = 3.74618322e-01;
576  rot_mat2[2][3] = -8.74132032e-01;
577  rot_mat2[2][4] = 1.57504581e-01;
578  rot_mat2[2][5] = 4.78555353e-02;
579  rot_mat2[2][6] = 1.23455215e-01;
580  rot_mat2[3][0] = -7.68247448e-02;
581  rot_mat2[3][1] = -2.48006571e-01;
582  rot_mat2[3][2] = -7.30866440e-02;
583  rot_mat2[3][3] = -2.06817660e-01;
584  rot_mat2[3][4] = -7.99320897e-01;
585  rot_mat2[3][5] = -4.53458397e-01;
586  rot_mat2[3][6] = -1.99842648e-01;
587  rot_mat2[4][0] = 2.59360916e-02;
588  rot_mat2[4][1] = 5.64417889e-01;
589  rot_mat2[4][2] = -1.05874452e-01;
590  rot_mat2[4][3] = -1.09701962e-01;
591  rot_mat2[4][4] = 2.76429977e-01;
592  rot_mat2[4][5] = -7.59899812e-01;
593  rot_mat2[4][6] = 6.04531910e-02;
594  rot_mat2[5][0] = 3.50608068e-01;
595  rot_mat2[5][1] = 1.58922313e-01;
596  rot_mat2[5][2] = -2.49914906e-01;
597  rot_mat2[5][3] = -1.49783413e-01;
598  rot_mat2[5][4] = 1.03865829e-01;
599  rot_mat2[5][5] = 1.56180926e-01;
600  rot_mat2[5][6] = -8.55420691e-01;
601  rot_mat2[6][0] = 8.44230723e-01;
602  rot_mat2[6][1] = 8.33975900e-02;
603  rot_mat2[6][2] = -2.35816939e-01;
604  rot_mat2[6][3] = 8.35859456e-02;
605  rot_mat2[6][4] = -2.38381031e-01;
606  rot_mat2[6][5] = 5.64243083e-02;
607  rot_mat2[6][6] = 3.97132056e-01;
608  spo->storeParamsBeforeRotation(); // avoids coefs_copy_ nullptr segfault
609  spo->applyRotation(rot_mat2, false); // Apply rotation2
610 
611  // Total rotation is product of rot_mat1 and rot_mat2
612  SPOSet::ValueMatrix rot_mat_tot(orbitalsetsize, orbitalsetsize);
613  for (int i = 0; i < orbitalsetsize; i++)
614  for (int j = 0; j < orbitalsetsize; j++)
615  {
616  rot_mat_tot[i][j] = 0;
617  for (int k = 0; k < orbitalsetsize; k++)
618  rot_mat_tot[i][j] += rot_mat1[i][k] * rot_mat2[k][j];
619  }
620 
621  // Get the data for rotated orbitals
622  SPOSet::ValueMatrix psiM_rot(elec_.R.size(), orbitalsetsize);
623  SPOSet::GradMatrix dpsiM_rot(elec_.R.size(), orbitalsetsize);
624  SPOSet::ValueMatrix d2psiM_rot(elec_.R.size(), orbitalsetsize);
625  spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM_rot, dpsiM_rot, d2psiM_rot);
626 
627  // Compute the expected value, gradient, and laplacian by hand using the transformation above
628  // Check value
629  SPOSet::ValueMatrix psiM_rot_manual(elec_.R.size(), orbitalsetsize);
630  for (int i = 0; i < elec_.R.size(); i++)
631  for (int j = 0; j < orbitalsetsize; j++)
632  {
633  psiM_rot_manual[i][j] = 0.;
634  for (int k = 0; k < orbitalsetsize; k++)
635  psiM_rot_manual[i][j] += psiM_bare[i][k] * rot_mat_tot[k][j];
636  }
637  auto check = checkMatrix(psiM_rot_manual, psiM_rot, true);
638  CHECKED_ELSE(check.result) { FAIL(check.result_message); }
639 
640  // Check grad
641  SPOSet::GradMatrix dpsiM_rot_manual(elec_.R.size(), orbitalsetsize);
642  for (int i = 0; i < elec_.R.size(); i++)
643  for (int j = 0; j < orbitalsetsize; j++)
644  {
645  dpsiM_rot_manual[i][j][0] = 0.;
646  dpsiM_rot_manual[i][j][1] = 0.;
647  dpsiM_rot_manual[i][j][2] = 0.;
648  for (int k = 0; k < orbitalsetsize; k++)
649  for (int l = 0; l < 3; l++)
650  dpsiM_rot_manual[i][j][l] += dpsiM_bare[i][k][l] * rot_mat_tot[k][j];
651  }
652 
653  // No checkMatrix for tensors? Gotta do it by hand...
654  double res = 0.;
655  for (int i = 0; i < elec_.R.size(); i++)
656  for (int j = 0; j < orbitalsetsize; j++)
657  for (int k = 0; k < 3; k++)
658  res += std::sqrt(std::norm(dpsiM_rot_manual[i][j][k] - dpsiM_rot[i][j][k]));
659 
660  CHECK(res == Approx(0.).epsilon(2e-4)); // increase threshold to accomodate mixed precision.
661 
662  // Check laplacian
663  SPOSet::ValueMatrix d2psiM_rot_manual(elec_.R.size(), orbitalsetsize);
664  for (int i = 0; i < elec_.R.size(); i++)
665  for (int j = 0; j < orbitalsetsize; j++)
666  {
667  d2psiM_rot_manual[i][j] = 0.;
668  for (int k = 0; k < orbitalsetsize; k++)
669  d2psiM_rot_manual[i][j] += d2psiM_bare[i][k] * rot_mat_tot[k][j];
670  }
671  check = checkMatrix(d2psiM_rot_manual, d2psiM_rot, true, 2e-4);
672  CHECKED_ELSE(check.result) { FAIL(check.result_message); }
673 
674 } // TEST_CASE
675 
676 
677 #ifdef QMC_COMPLEX
678 TEST_CASE("Spline applyRotation complex rotation", "[wavefunction]")
679 {
680  // How to get rid of all this annoying boilerplate?
681  Communicate* c = OHMMS::Controller;
682 
683  // diamondC_1x1x1
684  ParticleSet::ParticleLayout lattice;
685  lattice.R = {3.37316115, 3.37316115, 0.0, 0.0, 3.37316115, 3.37316115, 3.37316115, 0.0, 3.37316115};
686 
687  ParticleSetPool ptcl = ParticleSetPool(c);
688  ptcl.setSimulationCell(lattice);
689  auto ions_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
690  auto elec_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
691  ParticleSet& ions_(*ions_uptr);
692  ParticleSet& elec_(*elec_uptr);
693 
694  ions_.setName("ion");
695  ptcl.addParticleSet(std::move(ions_uptr));
696  ions_.create({2});
697  ions_.R[0] = {0.0, 0.0, 0.0};
698  ions_.R[1] = {1.68658058, 1.68658058, 1.68658058};
699 
700  elec_.setName("elec");
701  ptcl.addParticleSet(std::move(elec_uptr));
702  elec_.create({2});
703  elec_.R[0] = {0.0, 0.0, 0.0};
704  elec_.R[1] = {0.0, 1.0, 0.0};
705 
706  SpeciesSet& tspecies = elec_.getSpeciesSet();
707  int upIdx = tspecies.addSpecies("u");
708  int chargeIdx = tspecies.addAttribute("charge");
709  tspecies(chargeIdx, upIdx) = -1;
710 
711  // Load diamondC_1x1x1 wfn and explicitly construct a SplineC2C object with 7 orbitals
712  // This results in padding of the spline coefs table and thus is a more stringent test.
713  const char* particles = R"(<tmp>
714 <determinantset type="einspline" href="diamondC_1x1x1.pwscf.h5" tilematrix="1 0 0 0 1 0 0 0 1" twistnum="0" source="ion" gpu="no" precision="double" size="7"/>
715 </tmp>)";
716 
717  Libxml2Document doc;
718  bool okay = doc.parseFromString(particles);
719  REQUIRE(okay);
720  xmlNodePtr root = doc.getRoot();
721  xmlNodePtr ein1 = xmlFirstElementChild(root);
722  EinsplineSetBuilder einSet(elec_, ptcl.getPool(), c, ein1);
723  auto spo = einSet.createSPOSetFromXML(ein1);
724  REQUIRE(spo);
725 
726  const auto orbitalsetsize = spo->getOrbitalSetSize();
727  REQUIRE(orbitalsetsize == 7);
728  SPOSet::ValueMatrix psiM_bare(elec_.R.size(), orbitalsetsize);
729  SPOSet::GradMatrix dpsiM_bare(elec_.R.size(), orbitalsetsize);
730  SPOSet::ValueMatrix d2psiM_bare(elec_.R.size(), orbitalsetsize);
731  spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM_bare, dpsiM_bare, d2psiM_bare);
732 
733  // Check before rotation is applied. From 'test_einset.cpp'
734  CHECK(std::real(psiM_bare[1][0]) == Approx(-0.8886948824));
735  CHECK(std::real(psiM_bare[1][1]) == Approx(1.4194120169));
736  // grad
737  CHECK(std::real(dpsiM_bare[1][0][0]) == Approx(-0.0000183403));
738  CHECK(std::real(dpsiM_bare[1][0][1]) == Approx(0.1655139178));
739  CHECK(std::real(dpsiM_bare[1][0][2]) == Approx(-0.0000193077));
740  CHECK(std::real(dpsiM_bare[1][1][0]) == Approx(-1.3131694794));
741  CHECK(std::real(dpsiM_bare[1][1][1]) == Approx(-1.1174004078));
742  CHECK(std::real(dpsiM_bare[1][1][2]) == Approx(-0.8462534547));
743  // lapl
744  CHECK(std::real(d2psiM_bare[1][0]) == Approx(1.3313053846));
745  CHECK(std::real(d2psiM_bare[1][1]) == Approx(-4.712583065));
746 
747  /* Apply a rotation, U = exp(-K) as shown below, to the splines.
748  Because K = -K.T, U is unitary.
749  K= (0+i)*eye(orbitalsetsize)
750  U= (5.40302306e-01 -8.41470985e-01i) * eye(orbitalsetsize)
751  U * conj(U.T) = eye(7)
752 
753  NB: There's nothing special about this rotation. I purposefully made something 'ugly' to make a worst case scenario...
754  */
755  SPOSet::ValueMatrix rot_mat(orbitalsetsize, orbitalsetsize);
756  const SPOSet::ValueType z(5.40302306e-01, -8.41470985e-01);
757  rot_mat = 0;
758  for (int i = 0; i < orbitalsetsize; i++)
759  rot_mat[i][i] = z;
760 
761  spo->storeParamsBeforeRotation(); // avoids coefs_copy_ nullptr segfault
762  spo->applyRotation(rot_mat, false);
763 
764  // Get the data for rotated orbitals
765  SPOSet::ValueMatrix psiM_rot(elec_.R.size(), orbitalsetsize);
766  SPOSet::GradMatrix dpsiM_rot(elec_.R.size(), orbitalsetsize);
767  SPOSet::ValueMatrix d2psiM_rot(elec_.R.size(), orbitalsetsize);
768  spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM_rot, dpsiM_rot, d2psiM_rot);
769 
770  // Compute the expected value, gradient, and laplacian by hand using the transformation above
771  // Check value
772  SPOSet::ValueMatrix psiM_rot_manual(elec_.R.size(), orbitalsetsize);
773  for (int i = 0; i < elec_.R.size(); i++)
774  for (int j = 0; j < orbitalsetsize; j++)
775  {
776  psiM_rot_manual[i][j] = 0.;
777  for (int k = 0; k < orbitalsetsize; k++)
778  psiM_rot_manual[i][j] += psiM_bare[i][k] * rot_mat[k][j];
779  }
780  auto check = checkMatrix(psiM_rot_manual, psiM_rot, true);
781  CHECKED_ELSE(check.result) { FAIL(check.result_message); }
782 
783  // Check grad
784  SPOSet::GradMatrix dpsiM_rot_manual(elec_.R.size(), orbitalsetsize);
785  for (int i = 0; i < elec_.R.size(); i++)
786  for (int j = 0; j < orbitalsetsize; j++)
787  {
788  dpsiM_rot_manual[i][j][0] = 0.;
789  dpsiM_rot_manual[i][j][1] = 0.;
790  dpsiM_rot_manual[i][j][2] = 0.;
791  for (int k = 0; k < orbitalsetsize; k++)
792  for (int l = 0; l < 3; l++)
793  dpsiM_rot_manual[i][j][l] += dpsiM_bare[i][k][l] * rot_mat[k][j];
794  }
795 
796  // No checkMatrix for tensors? Gotta do it by hand...
797  double res = 0.;
798  for (int i = 0; i < elec_.R.size(); i++)
799  for (int j = 0; j < orbitalsetsize; j++)
800  for (int k = 0; k < 3; k++)
801  res += std::sqrt(std::norm(dpsiM_rot_manual[i][j][k] - dpsiM_rot[i][j][k]));
802 
803  CHECK(res == Approx(0.).epsilon(2e-4)); // increase threshold to accomodate mixed precision.
804 
805  // Check laplacian
806  SPOSet::ValueMatrix d2psiM_rot_manual(elec_.R.size(), orbitalsetsize);
807  for (int i = 0; i < elec_.R.size(); i++)
808  for (int j = 0; j < orbitalsetsize; j++)
809  {
810  d2psiM_rot_manual[i][j] = 0.;
811  for (int k = 0; k < orbitalsetsize; k++)
812  d2psiM_rot_manual[i][j] += d2psiM_bare[i][k] * rot_mat[k][j];
813  }
814  check = checkMatrix(d2psiM_rot_manual, d2psiM_rot, true, 2e-4);
815  CHECKED_ELSE(check.result) { FAIL(check.result_message); }
816 } // TEST_CASE
817 #endif
818 } // namespace qmcplusplus
qmcplusplus::CrystalLattice
a class that defines a supercell in D-dimensional Euclean space.
Definition: CrystalLattice.h:55
qmcplusplus::ParticleSet::setName
void setName(const std::string &aname)
Definition: ParticleSet.h:237
qmcplusplus::ParticleSetPool::setSimulationCell
void setSimulationCell(const SimulationCell &simulation_cell)
set simulation cell
Definition: ParticleSetPool.h:98
Libxml2Document
class that handles xmlDoc
Definition: Libxml2Doc.h:76
qmcplusplus::SpeciesSet::addSpecies
int addSpecies(const std::string &aname)
When a name species does not exist, add a new species.
Definition: SpeciesSet.cpp:33
qmcplusplus::Units::real
QMCTraits::RealType real
Definition: unit_conversion.h:23
qmcplusplus
helper functions for EinsplineSetBuilder
Definition: Configuration.h:43
qmcplusplus::PtclOnLatticeTraits::ParticleLayout
CrystalLattice< OHMMS_PRECISION, OHMMS_DIM > ParticleLayout
Definition: Configuration.h:79
qmcplusplus::TEST_CASE
TEST_CASE("complex_helper", "[type_traits]")
Definition: test_complex_helper.cpp:38
qmcplusplus::ParticleSetPool::addParticleSet
void addParticleSet(std::unique_ptr< ParticleSet > &&p)
add a ParticleSet* to the pool with its ownership transferred ParticleSet built outside the ParticleS...
Definition: ParticleSetPool.cpp:67
Libxml2Document::getRoot
xmlNodePtr getRoot()
Definition: Libxml2Doc.h:88
EinsplineSetBuilder.h
Builder class for einspline-based SPOSet objects.
qmcplusplus::SPOSet::ValueMatrix
OrbitalSetTraits< ValueType >::ValueMatrix ValueMatrix
Definition: SPOSet.h:50
qmcplusplus::CHECKED_ELSE
CHECKED_ELSE(check_matrix_result.result)
Definition: test_cuBLAS_LU.cpp:455
EinsplineSpinorSetBuilder.h
Derives EinsplineSetBuilder.
OHMMS::Controller
Communicate * Controller
Global Communicator for a process.
Definition: Communicate.cpp:35
qmcplusplus::SpeciesSet::addAttribute
int addAttribute(const std::string &aname)
for a new attribute, allocate the data, !More often used to get the index of a species ...
Definition: SpeciesSet.cpp:45
qmcplusplus::SPOSet::GradMatrix
OrbitalSetTraits< ValueType >::GradMatrix GradMatrix
Definition: SPOSet.h:52
Communicate
Wrapping information on parallelism.
Definition: Communicate.h:68
qmcplusplus::lattice
CrystalLattice< OHMMS_PRECISION, OHMMS_DIM > lattice
Definition: test_SpinDensityNew.cpp:142
qmcplusplus::ParticleSet
Specialized paritlce class for atomistic simulations.
Definition: ParticleSet.h:55
qmcplusplus::Vector::size
size_type size() const
return the current size
Definition: OhmmsVector.h:162
norm
double norm(const zVec &c)
Definition: VectorOps.h:118
qmcplusplus::QMCTraits::ValueType
QTBase::ValueType ValueType
Definition: Configuration.h:60
qmcplusplus::REQUIRE
REQUIRE(std::filesystem::exists(filename))
qmcplusplus::ParticleSetPool
Manage a collection of ParticleSet objects.
Definition: ParticleSetPool.h:34
qmcplusplus::ParticleSet::R
ParticlePos R
Position.
Definition: ParticleSet.h:79
qmcplusplus::ParticleSet::getSpeciesSet
SpeciesSet & getSpeciesSet()
retrun the SpeciesSet of this particle set
Definition: ParticleSet.h:231
qmcplusplus::ParticleSet::create
void create(const std::vector< int > &agroup)
create grouped particles
Definition: ParticleSet.cpp:125
qmcplusplus::sqrt
MakeReturn< UnaryNode< FnSqrt, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t sqrt(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:136
Libxml2Document::parseFromString
bool parseFromString(const std::string_view data)
Definition: Libxml2Doc.cpp:204
qmcplusplus::checkMatrix
CheckMatrixResult checkMatrix(M1 &a_mat, M2 &b_mat, const bool check_all=false, std::optional< const double > eps=std::nullopt)
This function checks equality a_mat and b_mat elements M1, M2 need to have their element type declare...
Definition: checkMatrix.hpp:63
qmcplusplus::CHECK
CHECK(log_values[0]==ComplexApprox(std::complex< double >{ 5.603777579195571, -6.1586603331188225 }))
WaveFunctionComponent.h
Declaration of WaveFunctionComponent.
qmcplusplus::SpeciesSet
Custom container for set of attributes for a set of species.
Definition: SpeciesSet.h:33
qmcplusplus::ParticleSetPool::getPool
const PoolType & getPool() const
get the Pool object
Definition: ParticleSetPool.h:92
qmcplusplus::EinsplineSetBuilder::createSPOSetFromXML
std::unique_ptr< SPOSet > createSPOSetFromXML(xmlNodePtr cur) override
initialize the Antisymmetric wave function for electrons
Definition: EinsplineSetBuilder_createSPOs.cpp:106
SplineC2C.h
class to handle complex splines to complex orbitals with splines of arbitrary precision ...
qmcplusplus::ParticleSetPool::getSimulationCell
const auto & getSimulationCell() const
get simulation cell
Definition: ParticleSetPool.h:95
ParticleSetPool.h
Declaration of ParticleSetPool.