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 "OhmmsData/Libxml2Doc.h"
16 #include "OhmmsPETE/OhmmsMatrix.h"
17 #include "Particle/ParticleSet.h"
18 #include "QMCWaveFunctions/WaveFunctionComponent.h"
19 #include "QMCWaveFunctions/Jastrow/BsplineFunctor.h"
20 #include "QMCWaveFunctions/Jastrow/RadialJastrowBuilder.h"
21 #include "ParticleBase/ParticleAttribOps.h"
22 #include "QMCWaveFunctions/Jastrow/J2OrbitalSoA.h"
23 #include "QMCWaveFunctions/Jastrow/J1OrbitalSoA.h"
24 
25 #include <stdio.h>
26 #include <string>
27 
28 
29 // Uncomment to print information and values from the underlying functor
30 //#define PRINT_SPLINE_DATA
31 
32 using std::string;
33 
34 namespace qmcplusplus
35 {
36 using RealType     = WaveFunctionComponent::RealType;
37 using PsiValueType = WaveFunctionComponent::PsiValueType;
38 
39 TEST_CASE("BSpline functor zero", "[wavefunction]")
40 {
41   BsplineFunctor<double> bf;
42 
43   double r = 1.2;
44   double u = bf.evaluate(r);
45   REQUIRE(u == 0.0);
46 }
47 
48 TEST_CASE("BSpline functor one", "[wavefunction]")
49 {
50   BsplineFunctor<double> bf;
51 
52   bf.resize(1);
53 
54   double r = 1.2;
55   double u = bf.evaluate(r);
56   REQUIRE(u == 0.0);
57 }
58 
59 TEST_CASE("BSpline builder Jastrow J2", "[wavefunction]")
60 {
61   Communicate* c;
62   c = OHMMS::Controller;
63 
64   ParticleSet ions_;
65   ParticleSet elec_;
66 
67   ions_.setName("ion");
68   ions_.create(1);
69   ions_.R[0][0] = 2.0;
70   ions_.R[0][1] = 0.0;
71   ions_.R[0][2] = 0.0;
72 
73   elec_.setName("elec");
74   std::vector<int> ud(2);
75   ud[0] = ud[1] = 1;
76   elec_.create(ud);
77   elec_.R[0][0] = 1.00;
78   elec_.R[0][1] = 0.0;
79   elec_.R[0][2] = 0.0;
80   elec_.R[1][0] = 0.0;
81   elec_.R[1][1] = 0.0;
82   elec_.R[1][2] = 0.0;
83 
84   SpeciesSet& tspecies         = elec_.getSpeciesSet();
85   int upIdx                    = tspecies.addSpecies("u");
86   int downIdx                  = tspecies.addSpecies("d");
87   int chargeIdx                = tspecies.addAttribute("charge");
88   tspecies(chargeIdx, upIdx)   = -1;
89   tspecies(chargeIdx, downIdx) = -1;
90   elec_.resetGroups();
91 
92   const char* particles = "<tmp> \
93 <jastrow name=\"J2\" type=\"Two-Body\" function=\"Bspline\" print=\"yes\"> \
94    <correlation rcut=\"10\" size=\"10\" speciesA=\"u\" speciesB=\"d\"> \
95       <coefficients id=\"ud\" type=\"Array\"> 0.02904699284 -0.1004179 -0.1752703883 -0.2232576505 -0.2728029201 -0.3253286875 -0.3624525145 -0.3958223107 -0.4268582166 -0.4394531176</coefficients> \
96     </correlation> \
97 </jastrow> \
98 </tmp> \
99 ";
100   Libxml2Document doc;
101   bool okay = doc.parseFromString(particles);
102   REQUIRE(okay);
103 
104   xmlNodePtr root = doc.getRoot();
105 
106   xmlNodePtr jas1 = xmlFirstElementChild(root);
107 
108   RadialJastrowBuilder jastrow(c, elec_);
109 
110   typedef J2OrbitalSoA<BsplineFunctor<RealType>> J2Type;
111   std::unique_ptr<J2Type> j2(dynamic_cast<J2Type*>(jastrow.buildComponent(jas1)));
112   REQUIRE(j2);
113 
114   // update all distance tables
115   elec_.update();
116 
117   double logpsi_real = std::real(j2->evaluateLog(elec_, elec_.G, elec_.L));
118   REQUIRE(logpsi_real == Approx(0.1012632641)); // note: number not validated
119 
120   double KE = -0.5 * (Dot(elec_.G, elec_.G) + Sum(elec_.L));
121   REQUIRE(KE == Approx(-0.1616624771)); // note: number not validated
122 
123 
124   // now test evaluateHessian
125   WaveFunctionComponent::HessVector_t grad_grad_psi;
126   grad_grad_psi.resize(elec_.getTotalNum());
127   grad_grad_psi = 0.0;
128 
129   std::cout << "eval hess" << std::endl;
130   j2->evaluateHessian(elec_, grad_grad_psi);
131   std::vector<double> hess_values = {
132       -0.0627236, 0, 0, 0, 0.10652, 0, 0, 0, 0.10652, -0.0627236, 0, 0, 0, 0.10652, 0, 0, 0, 0.10652,
133   };
134 
135   int m = 0;
136   for (int n = 0; n < elec_.getTotalNum(); n++)
137     for (int i = 0; i < OHMMS_DIM; i++)
138       for (int j = 0; j < OHMMS_DIM; j++, m++)
139       {
140         REQUIRE(std::real(grad_grad_psi[n](i, j)) == Approx(hess_values[m]));
141       }
142 
143 
144   struct JValues
145   {
146     double r;
147     double u;
148     double du;
149     double ddu;
150   };
151 
152   // Cut and paste from output of gen_bspline_jastrow.py
153   const int N     = 20;
154   JValues Vals[N] = {{0.00, 0.1374071801, -0.5, 0.7866949593},
155                      {0.60, -0.04952403966, -0.1706645865, 0.3110897524},
156                      {1.20, -0.121361995, -0.09471371432, 0.055337302},
157                      {1.80, -0.1695590431, -0.06815900213, 0.0331784053},
158                      {2.40, -0.2058414025, -0.05505192964, 0.01049597156},
159                      {3.00, -0.2382237097, -0.05422744821, -0.002401552969},
160                      {3.60, -0.2712606182, -0.05600918024, -0.003537553803},
161                      {4.20, -0.3047843679, -0.05428535477, 0.0101841028},
162                      {4.80, -0.3347515004, -0.04506573714, 0.01469003611},
163                      {5.40, -0.3597048574, -0.03904232165, 0.005388015505},
164                      {6.00, -0.3823503292, -0.03657502025, 0.003511355265},
165                      {6.60, -0.4036800017, -0.03415678101, 0.007891305516},
166                      {7.20, -0.4219818468, -0.02556305518, 0.02075444724},
167                      {7.80, -0.4192355508, 0.06799438701, 0.3266190181},
168                      {8.40, -0.3019238309, 0.32586994, 0.2880861726},
169                      {9.00, -0.09726352421, 0.2851358014, -0.4238666348},
170                      {9.60, -0.006239062395, 0.04679296796, -0.2339648398},
171                      {10.20, 0, 0, 0},
172                      {10.80, 0, 0, 0},
173                      {11.40, 0, 0, 0}};
174 
175 
176   BsplineFunctor<RealType>* bf = j2->F[0];
177 
178   for (int i = 0; i < N; i++)
179   {
180     RealType dv  = 0.0;
181     RealType ddv = 0.0;
182     RealType val = bf->evaluate(Vals[i].r, dv, ddv);
183     REQUIRE(Vals[i].u == Approx(val));
184     REQUIRE(Vals[i].du == Approx(dv));
185     REQUIRE(Vals[i].ddu == Approx(ddv));
186   }
187 
188 #ifdef PRINT_SPLINE_DATA
189   // write out values of the Bspline functor
190   //BsplineFunctor<double> *bf = j2->F[0];
191   printf("NumParams = %d\n", bf->NumParams);
192   printf("CuspValue = %g\n", bf->CuspValue);
193   printf("DeltaR = %g\n", bf->DeltaR);
194   printf("SplineCoeffs size = %d\n", bf->SplineCoefs.size());
195   for (int j = 0; j < bf->SplineCoefs.size(); j++)
196   {
197     printf("%d %g\n", j, bf->SplineCoefs[j]);
198   }
199   printf("\n");
200 
201   for (int i = 0; i < 20; i++)
202   {
203     double r      = 0.6 * i;
204     elec_.R[0][0] = r;
205     elec_.update();
206     double logpsi_real = std::real(j2->evaluateLog(elec_, elec_.G, elec_.L));
207     //double alt_val = bf->evaluate(r);
208     double dv      = 0.0;
209     double ddv     = 0.0;
210     double alt_val = bf->evaluate(r, dv, ddv);
211     printf("%g %g %g %g %g\n", r, logpsi_real, alt_val, dv, ddv);
212   }
213 #endif
214 
215   typedef QMCTraits::ValueType ValueType;
216   typedef QMCTraits::PosType PosType;
217 
218   // set virtutal particle position
219   PosType newpos(0.3, 0.2, 0.5);
220 
221   elec_.makeVirtualMoves(newpos);
222   std::vector<ValueType> ratios(elec_.getTotalNum());
223   j2->evaluateRatiosAlltoOne(elec_, ratios);
224 
225   REQUIRE(std::real(ratios[0]) == Approx(0.9522052017));
226   REQUIRE(std::real(ratios[1]) == Approx(0.9871985577));
227 
228   elec_.makeMove(0, newpos - elec_.R[0]);
229   PsiValueType ratio_0 = j2->ratio(elec_, 0);
230   elec_.rejectMove(0);
231 
232   REQUIRE(std::real(ratio_0) == Approx(0.9522052017));
233 
234   VirtualParticleSet VP(elec_, 2);
235   std::vector<PosType> newpos2(2);
236   std::vector<ValueType> ratios2(2);
237   newpos2[0] = newpos - elec_.R[1];
238   newpos2[1] = PosType(0.2, 0.5, 0.3) - elec_.R[1];
239   VP.makeMoves(1, elec_.R[1], newpos2);
240   j2->evaluateRatios(VP, ratios);
241 
242   REQUIRE(std::real(ratios[0]) == Approx(0.9871985577));
243   REQUIRE(std::real(ratios[1]) == Approx(0.9989268241));
244 
245   //test acceptMove
246   elec_.makeMove(1, newpos - elec_.R[1]);
247   PsiValueType ratio_1 = j2->ratio(elec_, 1);
248   j2->acceptMove(elec_, 1);
249   elec_.acceptMove(1);
250 
251   REQUIRE(std::real(ratio_1) == Approx(0.9871985577));
252   REQUIRE(std::real(j2->LogValue) == Approx(0.0883791773));
253 }
254 
255 TEST_CASE("BSpline builder Jastrow J1", "[wavefunction]")
256 {
257   Communicate* c;
258   c = OHMMS::Controller;
259 
260   ParticleSet ions_;
261   ParticleSet elec_;
262 
263   ions_.setName("ion");
264   ions_.create(1);
265   ions_.R[0][0] = 2.0;
266   ions_.R[0][1] = 0.0;
267   ions_.R[0][2] = 0.0;
268 
269   SpeciesSet& ispecies       = ions_.getSpeciesSet();
270   int CIdx                   = ispecies.addSpecies("C");
271   int ichargeIdx             = ispecies.addAttribute("charge");
272   ispecies(ichargeIdx, CIdx) = 4;
273   ions_.resetGroups();
274   ions_.update();
275 
276   elec_.setName("elec");
277   std::vector<int> ud(2);
278   ud[0] = ud[1] = 1;
279   elec_.create(ud);
280   elec_.R[0][0] = 1.00;
281   elec_.R[0][1] = 0.0;
282   elec_.R[0][2] = 0.0;
283   elec_.R[1][0] = 0.0;
284   elec_.R[1][1] = 0.0;
285   elec_.R[1][2] = 0.0;
286 
287   SpeciesSet& tspecies         = elec_.getSpeciesSet();
288   int upIdx                    = tspecies.addSpecies("u");
289   int downIdx                  = tspecies.addSpecies("d");
290   int chargeIdx                = tspecies.addAttribute("charge");
291   tspecies(chargeIdx, upIdx)   = -1;
292   tspecies(chargeIdx, downIdx) = -1;
293   elec_.resetGroups();
294 
295   const char* particles = "<tmp> \
296    <jastrow type=\"One-Body\" name=\"J1\" function=\"bspline\" source=\"ion\" print=\"yes\"> \
297        <correlation elementType=\"C\" rcut=\"10\" size=\"8\" cusp=\"0.0\"> \
298                <coefficients id=\"eC\" type=\"Array\"> \
299 -0.2032153051 -0.1625595974 -0.143124599 -0.1216434956 -0.09919771951 -0.07111729038 \
300 -0.04445345869 -0.02135082917 \
301                </coefficients> \
302             </correlation> \
303          </jastrow> \
304 </tmp> \
305 ";
306   Libxml2Document doc;
307   bool okay = doc.parseFromString(particles);
308   REQUIRE(okay);
309 
310   xmlNodePtr root = doc.getRoot();
311 
312   xmlNodePtr jas1 = xmlFirstElementChild(root);
313 
314   RadialJastrowBuilder jastrow(c, elec_, ions_);
315 
316   typedef J1OrbitalSoA<BsplineFunctor<RealType>> J1Type;
317   std::unique_ptr<J1Type> j1(dynamic_cast<J1Type*>(jastrow.buildComponent(jas1)));
318   REQUIRE(j1);
319 
320   // update all distance tables
321   elec_.update();
322 
323   double logpsi_real = std::real(j1->evaluateLog(elec_, elec_.G, elec_.L));
324   REQUIRE(logpsi_real == Approx(0.3160552244)); // note: number not validated
325 
326   //Ionic Derivative Test.
327   QMCTraits::GradType gsource(0.0);
328   TinyVector<ParticleSet::ParticleGradient_t, OHMMS_DIM> grad_grad_source;
329   TinyVector<ParticleSet::ParticleLaplacian_t, OHMMS_DIM> lapl_grad_source;
330   int nelecs = elec_.getTotalNum();
331 
332   for (int dim = 0; dim < OHMMS_DIM; dim++)
333   {
334     grad_grad_source[dim].resize(nelecs);
335     lapl_grad_source[dim].resize(nelecs);
336   }
337 
338   /////////////////////////////////////////
339   //Testing the ion gradient w.r.t. ion # 0
340   /////////////////////////////////////////
341   //NOTE:  All test values in this section are validated against finite differences
342   //       for this configuration.
343 
344   //First we test evalGradSource(P,ions,ionid);
345 
346   gsource = j1->evalGradSource(elec_, ions_, 0);
347 
348   //Gradient comparison
349   REQUIRE(std::real(gsource[0]) == Approx(-0.04695203659));
350   REQUIRE(std::real(gsource[1]) == Approx(0.00000000000));
351   REQUIRE(std::real(gsource[2]) == Approx(0.00000000000));
352 
353   //Now we test evalGradSource that returns higher order derivatives.
354   gsource = j1->evalGradSource(elec_, ions_, 0, grad_grad_source, lapl_grad_source);
355 
356   //Gradient comparison
357   REQUIRE(std::real(gsource[0]) == Approx(-0.04695203659));
358   REQUIRE(std::real(gsource[1]) == Approx(0.00000000000));
359   REQUIRE(std::real(gsource[2]) == Approx(0.00000000000));
360 
361   //Ion gradient of electron gradient comparison.
362   REQUIRE(std::real(grad_grad_source[0][0][0]) == Approx(-0.008883672));
363   REQUIRE(std::real(grad_grad_source[0][1][0]) == Approx(-0.002111879));
364   REQUIRE(std::real(grad_grad_source[1][0][1]) == Approx(0.028489287));
365   REQUIRE(std::real(grad_grad_source[1][1][1]) == Approx(0.009231375));
366   REQUIRE(std::real(grad_grad_source[2][0][2]) == Approx(0.028489287));
367   REQUIRE(std::real(grad_grad_source[2][1][2]) == Approx(0.009231375));
368 
369   //Ion gradient of electron laplacians.
370   REQUIRE(std::real(lapl_grad_source[0][0]) == Approx(0.1494918378));
371   REQUIRE(std::real(lapl_grad_source[0][1]) == Approx(-0.0056182539));
372   REQUIRE(std::real(lapl_grad_source[1][0]) == Approx(0.0000000000));
373   REQUIRE(std::real(lapl_grad_source[1][1]) == Approx(0.0000000000));
374   REQUIRE(std::real(lapl_grad_source[2][0]) == Approx(0.0000000000));
375   REQUIRE(std::real(lapl_grad_source[2][1]) == Approx(0.0000000000));
376 
377 
378   // now test evaluateHessian
379   WaveFunctionComponent::HessVector_t grad_grad_psi;
380   grad_grad_psi.resize(elec_.getTotalNum());
381   grad_grad_psi = 0.0;
382 
383   j1->evaluateHessian(elec_, grad_grad_psi);
384 
385   std::vector<double> hess_values = {
386       0.00888367, 0, 0, 0, -0.0284893, 0, 0, 0, -0.0284893, 0.00211188, 0, 0, 0, -0.00923137, 0, 0, 0, -0.00923137,
387   };
388 
389   int m = 0;
390   for (int n = 0; n < elec_.getTotalNum(); n++)
391     for (int i = 0; i < OHMMS_DIM; i++)
392       for (int j = 0; j < OHMMS_DIM; j++, m++)
393       {
394         REQUIRE(std::real(grad_grad_psi[n](i, j)) == Approx(hess_values[m]));
395       }
396 
397   j1->evaluateLog(elec_, elec_.G, elec_.L); // evaluateHessian has side effects
398 
399 
400   struct JValues
401   {
402     double r;
403     double u;
404     double du;
405     double ddu;
406   };
407 
408   // Cut and paste from output of gen_bspline_jastrow.py
409   const int N     = 20;
410   JValues Vals[N] = {{0.00, -0.1896634025, 0, 0.06586224647},
411                      {0.60, -0.1804990512, 0.02606308248, 0.02101469513},
412                      {1.20, -0.1637586749, 0.0255799351, -0.01568108497},
413                      {1.80, -0.1506226948, 0.01922435549, -0.005504180392},
414                      {2.40, -0.1394848415, 0.01869442683, 0.001517191423},
415                      {3.00, -0.128023472, 0.01946283614, 0.00104417293},
416                      {3.60, -0.1161729491, 0.02009651096, 0.001689229059},
417                      {4.20, -0.1036884223, 0.02172284322, 0.003731878464},
418                      {4.80, -0.08992443283, 0.0240346508, 0.002736384838},
419                      {5.40, -0.07519614609, 0.02475121662, -0.000347832122},
420                      {6.00, -0.06054074137, 0.02397053075, -0.001842295859},
421                      {6.60, -0.04654631918, 0.0225837382, -0.002780345968},
422                      {7.20, -0.03347994129, 0.02104406699, -0.00218107833},
423                      {7.80, -0.0211986378, 0.01996899618, -0.00173646255},
424                      {8.40, -0.01004416026, 0.01635533409, -0.01030907776},
425                      {9.00, -0.002594125744, 0.007782377232, -0.01556475446},
426                      {9.60, -0.0001660240476, 0.001245180357, -0.006225901786},
427                      {10.20, 0, 0, 0},
428                      {10.80, 0, 0, 0},
429                      {11.40, 0, 0, 0}};
430 
431   BsplineFunctor<RealType>* bf = j1->F[0];
432 
433   for (int i = 0; i < N; i++)
434   {
435     RealType dv  = 0.0;
436     RealType ddv = 0.0;
437     RealType val = bf->evaluate(Vals[i].r, dv, ddv);
438     REQUIRE(Vals[i].u == Approx(val));
439     REQUIRE(Vals[i].du == Approx(dv));
440     REQUIRE(Vals[i].ddu == Approx(ddv));
441   }
442 
443 #ifdef PRINT_SPLINE_DATA
444   // write out values of the Bspline functor
445   //BsplineFunctor<double> *bf = j1->F[0];
446   printf("NumParams = %d\n", bf->NumParams);
447   printf("CuspValue = %g\n", bf->CuspValue);
448   printf("DeltaR = %g\n", bf->DeltaR);
449   printf("SplineCoeffs size = %d\n", bf->SplineCoefs.size());
450   for (int j = 0; j < bf->SplineCoefs.size(); j++)
451   {
452     printf("%d %g\n", j, bf->SplineCoefs[j]);
453   }
454   printf("\n");
455 
456   for (int i = 0; i < 20; i++)
457   {
458     double r      = 0.6 * i;
459     elec_.R[0][0] = r;
460     elec_.update();
461     double logpsi_real = std::real(j1->evaluateLog(elec_, elec_.G, elec_.L));
462     //double alt_val = bf->evaluate(r);
463     double dv      = 0.0;
464     double ddv     = 0.0;
465     double alt_val = bf->evaluate(r, dv, ddv);
466     printf("%g %g %g %g %g\n", r, logpsi_real, alt_val, dv, ddv);
467   }
468 #endif
469 
470   typedef QMCTraits::ValueType ValueType;
471   typedef QMCTraits::PosType PosType;
472 
473   // set virtutal particle position
474   PosType newpos(0.3, 0.2, 0.5);
475 
476   elec_.makeVirtualMoves(newpos);
477   std::vector<ValueType> ratios(elec_.getTotalNum());
478   j1->evaluateRatiosAlltoOne(elec_, ratios);
479 
480   REQUIRE(std::real(ratios[0]) == Approx(0.9819208747));
481   REQUIRE(std::real(ratios[1]) == Approx(1.0040884258));
482 
483   elec_.makeMove(0, newpos - elec_.R[0]);
484   PsiValueType ratio_0 = j1->ratio(elec_, 0);
485   elec_.rejectMove(0);
486 
487   REQUIRE(std::real(ratio_0) == Approx(0.9819208747));
488 
489   // test acceptMove results
490   elec_.makeMove(1, newpos - elec_.R[1]);
491   PsiValueType ratio_1 = j1->ratio(elec_, 1);
492   j1->acceptMove(elec_, 1);
493   elec_.acceptMove(1);
494 
495   REQUIRE(std::real(ratio_1) == Approx(1.0040884258));
496   REQUIRE(std::real(j1->LogValue) == Approx(0.32013531536));
497 
498   // test to make sure that setting cusp for J1 works properly
499   const char* particles2 = "<tmp> \
500    <jastrow type=\"One-Body\" name=\"J1\" function=\"bspline\" source=\"ion\" print=\"yes\"> \
501        <correlation elementType=\"C\" rcut=\"10\" size=\"8\" cusp=\"2.0\"> \
502                <coefficients id=\"eC\" type=\"Array\"> \
503 -0.2032153051 -0.1625595974 -0.143124599 -0.1216434956 -0.09919771951 -0.07111729038 \
504 -0.04445345869 -0.02135082917 \
505                </coefficients> \
506             </correlation> \
507          </jastrow> \
508 </tmp> \
509 ";
510 
511   Libxml2Document doc2;
512   bool okay2 = doc2.parseFromString(particles2);
513   REQUIRE(okay2);
514 
515   xmlNodePtr root2 = doc2.getRoot();
516 
517   xmlNodePtr jas2 = xmlFirstElementChild(root2);
518 
519   RadialJastrowBuilder jastrow2(c, elec_, ions_);
520 
521   std::unique_ptr<J1Type> j12(dynamic_cast<J1Type*>(jastrow2.buildComponent(jas2)));
522   REQUIRE(j12);
523 
524   // Cut and paste from output of gen_bspline_jastrow.py
525   // note only the first two rows should change from above
526   const int N2      = 20;
527   JValues Vals2[N2] = {{0.00, -0.9304041433, 2, -3.534137754},
528                        {0.60, -0.252599792, 0.4492630825, -1.634985305},
529                        {1.20, -0.1637586749, 0.0255799351, -0.01568108497},
530                        {1.80, -0.1506226948, 0.01922435549, -0.005504180392},
531                        {2.40, -0.1394848415, 0.01869442683, 0.001517191423},
532                        {3.00, -0.128023472, 0.01946283614, 0.00104417293},
533                        {3.60, -0.1161729491, 0.02009651096, 0.001689229059},
534                        {4.20, -0.1036884223, 0.02172284322, 0.003731878464},
535                        {4.80, -0.08992443283, 0.0240346508, 0.002736384838},
536                        {5.40, -0.07519614609, 0.02475121662, -0.000347832122},
537                        {6.00, -0.06054074137, 0.02397053075, -0.001842295859},
538                        {6.60, -0.04654631918, 0.0225837382, -0.002780345968},
539                        {7.20, -0.03347994129, 0.02104406699, -0.00218107833},
540                        {7.80, -0.0211986378, 0.01996899618, -0.00173646255},
541                        {8.40, -0.01004416026, 0.01635533409, -0.01030907776},
542                        {9.00, -0.002594125744, 0.007782377232, -0.01556475446},
543                        {9.60, -0.0001660240476, 0.001245180357, -0.006225901786},
544                        {10.20, 0, 0, 0},
545                        {10.80, 0, 0, 0},
546                        {11.40, 0, 0, 0}};
547 
548   BsplineFunctor<RealType>* bf2 = j12->F[0];
549 
550   for (int i = 0; i < N2; i++)
551   {
552     RealType dv  = 0.0;
553     RealType ddv = 0.0;
554     RealType val = bf2->evaluate(Vals2[i].r, dv, ddv);
555     REQUIRE(Vals2[i].du == Approx(dv));
556     REQUIRE(Vals2[i].u == Approx(val));
557     REQUIRE(Vals2[i].ddu == Approx(ddv));
558   }
559 }
560 } // namespace qmcplusplus
561