QMCPACK
test_soa_cusp_corr.cpp
Go to the documentation of this file.
1 
2 //////////////////////////////////////////////////////////////////////////////////////
3 // This file is distributed under the University of Illinois/NCSA Open Source License.
4 // See LICENSE file in top directory for details.
5 //
6 // Copyright (c) 2018 Jeongnim Kim and QMCPACK developers.
7 //
8 // File developed by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory
9 //
10 // File created by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory
11 //////////////////////////////////////////////////////////////////////////////////////
12 
13 
14 #include "catch.hpp"
15 
16 #include "Configuration.h"
17 
18 #include "Message/Communicate.h"
19 #include "Numerics/OneDimGridBase.h"
20 #include "ParticleIO/XMLParticleIO.h"
21 #include "Numerics/GaussianBasisSet.h"
22 
23 #include "QMCWaveFunctions/LCAO/LCAOrbitalSet.h"
24 #include "QMCWaveFunctions/LCAO/CuspCorrectionConstruction.h"
25 
26 #include "QMCWaveFunctions/SPOSetBuilderFactory.h"
27 
28 namespace qmcplusplus
29 {
30 TEST_CASE("readCuspInfo", "[wavefunction]")
31 {
32  Communicate* c = OHMMS::Controller;
33 
34  using GridType = OneDimGridBase<double>;
35 
36  Matrix<CuspCorrectionParameters> info;
37  int num_center = 3;
38  int orbital_set_size = 7;
39  info.resize(num_center, orbital_set_size);
40 
41  bool okay = readCuspInfo("hcn_downdet.cuspInfo.xml", "downdet", orbital_set_size, info);
42  REQUIRE(okay);
43 
44  // N
45  CHECK(info(0, 0).redo == Approx(0.0)); // redo
46  CHECK(info(0, 0).C == Approx(0.0)); // C
47  CHECK(info(0, 0).sg == Approx(1.0)); // sg
48  CHECK(info(0, 0).Rc == Approx(0.0769130700800000)); // rc
49  CHECK(info(0, 0).alpha[0] == Approx(2.29508580995773)); // a1
50  CHECK(info(0, 0).alpha[1] == Approx(-7.00028778782666)); // a2
51  CHECK(info(0, 0).alpha[2] == Approx(0.834942828252775)); // a3
52  CHECK(info(0, 0).alpha[3] == Approx(-4.61597420905980)); // a4
53  CHECK(info(0, 0).alpha[4] == Approx(31.6558091872316)); // a5
54 
55  // Spot check a few values from these centers
56  // C
57  CHECK(info(0, 6).C == Approx(0.0)); // C
58  CHECK(info(0, 6).alpha[4] == Approx(0.0)); // a5
59 
60  // H
61  CHECK(info(2, 4).alpha[4] == Approx(-404.733151049101)); // a5
62 }
63 
64 
65 TEST_CASE("applyCuspInfo", "[wavefunction]")
66 {
67  Communicate* c = OHMMS::Controller;
68 
69  Libxml2Document doc;
70  bool okay = doc.parse("hcn.structure.xml");
71  REQUIRE(okay);
72  xmlNodePtr root = doc.getRoot();
73 
74  const SimulationCell simulation_cell;
75  auto ions_ptr = std::make_unique<ParticleSet>(simulation_cell);
76  auto& ions(*ions_ptr);
77  XMLParticleParser parse_ions(ions);
78  OhmmsXPathObject particleset_ion("//particleset[@name='ion0']", doc.getXPathContext());
79  REQUIRE(particleset_ion.size() == 1);
80  parse_ions.readXML(particleset_ion[0]);
81 
82  REQUIRE(ions.groups() == 3);
83  REQUIRE(ions.R.size() == 3);
84  ions.update();
85 
86  auto elec_ptr = std::make_unique<ParticleSet>(simulation_cell);
87  auto& elec(*elec_ptr);
88  XMLParticleParser parse_elec(elec);
89  OhmmsXPathObject particleset_elec("//particleset[@name='e']", doc.getXPathContext());
90  REQUIRE(particleset_elec.size() == 1);
91  parse_elec.readXML(particleset_elec[0]);
92 
93  REQUIRE(elec.groups() == 2);
94  REQUIRE(elec.R.size() == 14);
95 
96  elec.R = 0.0;
97 
98  elec.addTable(ions);
99  elec.update();
100 
101  Libxml2Document doc2;
102  okay = doc2.parse("hcn.wfnoj.xml");
103  REQUIRE(okay);
104  xmlNodePtr root2 = doc2.getRoot();
105 
106  WaveFunctionComponentBuilder::PSetMap particle_set_map;
107  particle_set_map.emplace(elec_ptr->getName(), std::move(elec_ptr));
108  particle_set_map.emplace(ions_ptr->getName(), std::move(ions_ptr));
109 
110  SPOSetBuilderFactory bf(c, elec, particle_set_map);
111 
112  OhmmsXPathObject MO_base("//determinantset", doc2.getXPathContext());
113  REQUIRE(MO_base.size() == 1);
114 
115  const auto bb_ptr = bf.createSPOSetBuilder(MO_base[0]);
116  auto& bb(*bb_ptr);
117 
118  OhmmsXPathObject slater_base("//determinant", doc2.getXPathContext());
119  auto sposet = bb.createSPOSet(slater_base[0]);
120 
121  LCAOrbitalSet& lcob = dynamic_cast<LCAOrbitalSet&>(*sposet);
122 
123  LCAOrbitalSet phi("phi", std::unique_ptr<LCAOrbitalSet::basis_type>(lcob.myBasisSet->makeClone()),
124  lcob.getOrbitalSetSize(), lcob.isIdentity(), lcob.isOMPoffload());
125 
126  LCAOrbitalSet eta("eta", std::unique_ptr<LCAOrbitalSet::basis_type>(lcob.myBasisSet->makeClone()),
127  lcob.getOrbitalSetSize(), lcob.isIdentity(), lcob.isOMPoffload());
128 
129  *(eta.C) = *(lcob.C);
130  *(phi.C) = *(lcob.C);
131 
132 
133  int num_center = 3;
134  std::vector<bool> corrCenter(num_center, "true");
135 
136  // N is first atom
137  int center_idx = 0;
138 
139  using RealType = QMCTraits::RealType;
140 
141  splitPhiEta(center_idx, corrCenter, phi, eta);
142 
143  // 1S orbital on N
144  CHECK((*phi.C)(0, 0) == Approx(1.00180500));
145  CHECK((*eta.C)(0, 0) == Approx(0.0));
146 
147  int orbital_set_size = 7;
148  Matrix<CuspCorrectionParameters> info;
149  info.resize(num_center, orbital_set_size);
150  okay = readCuspInfo("hcn_downdet.cuspInfo.xml", "downdet", orbital_set_size, info);
151 
152  REQUIRE(okay);
153  Vector<RealType> xgrid;
154  Vector<RealType> rad_orb;
155  int ngrid = 10;
156  xgrid.resize(ngrid);
157  for (int i = 0; i < ngrid; i++)
158  {
159  xgrid[i] = 0.012 * (i + 1);
160  }
161 
162  rad_orb.resize(ngrid);
163 
164  int mo_idx = 0;
165  computeRadialPhiBar(&elec, &ions, mo_idx, center_idx, &phi, xgrid, rad_orb, info(center_idx, mo_idx));
166 
167  // Comparisons generated from gen_cusp_corr.py
168  // Center 0 MO 0 rc = 0.07691307008
169  CHECK(rad_orb[0] == Approx(9.1266186340)); // x = 0.012
170  CHECK(rad_orb[1] == Approx(8.3939106599)); // x = 0.024
171  CHECK(rad_orb[2] == Approx(7.7213972780)); // x = 0.036
172  CHECK(rad_orb[3] == Approx(7.1039662640)); // x = 0.048
173  CHECK(rad_orb[4] == Approx(6.5370601478)); // x = 0.06
174  CHECK(rad_orb[5] == Approx(6.0165935481)); // x = 0.072
175  CHECK(rad_orb[6] == Approx(5.5390213984)); // x = 0.084
176  CHECK(rad_orb[7] == Approx(5.1023814795)); // x = 0.096
177  CHECK(rad_orb[8] == Approx(4.7033287383)); // x = 0.108
178  CHECK(rad_orb[9] == Approx(4.3370522377)); // x = 0.12
179 
180 
181  mo_idx = 1;
182  computeRadialPhiBar(&elec, &ions, mo_idx, center_idx, &phi, xgrid, rad_orb, info(center_idx, mo_idx));
183 
184  // Center 0 MO 1 rc = 0.060909477888
185  CHECK(rad_orb[0] == Approx(-0.0099816961)); // x = 0.012
186  CHECK(rad_orb[1] == Approx(-0.0092950723)); // x = 0.024
187  CHECK(rad_orb[2] == Approx(-0.0086498844)); // x = 0.036
188  CHECK(rad_orb[3] == Approx(-0.0080440071)); // x = 0.048
189  CHECK(rad_orb[4] == Approx(-0.0074778482)); // x = 0.06
190  CHECK(rad_orb[5] == Approx(-0.0069529708)); // x = 0.072
191  CHECK(rad_orb[6] == Approx(-0.0064707256)); // x = 0.084
192  CHECK(rad_orb[7] == Approx(-0.0060313791)); // x = 0.096
193  CHECK(rad_orb[8] == Approx(-0.0056312867)); // x = 0.108
194  CHECK(rad_orb[9] == Approx(-0.0052652668)); // x = 0.12
195 
196 
197  // Reset the MO matrices for another center
198 
199  *(eta.C) = *(lcob.C);
200  *(phi.C) = *(lcob.C);
201 
202 
203  // C is second atom
204  center_idx = 1;
205  splitPhiEta(center_idx, corrCenter, phi, eta);
206 
207  // 1S orbital on N
208  CHECK((*phi.C)(0, 0) == Approx(0.0));
209  CHECK((*eta.C)(0, 0) == Approx(1.00180500));
210 
211  mo_idx = 0;
212  computeRadialPhiBar(&elec, &ions, mo_idx, center_idx, &phi, xgrid, rad_orb, info(center_idx, mo_idx));
213 
214  // Center 1 MO 0 rc = 0.105
215  CHECK(rad_orb[0] == Approx(0.0017535517)); // x = 0.012
216  CHECK(rad_orb[1] == Approx(0.0016496533)); // x = 0.024
217  CHECK(rad_orb[2] == Approx(0.0015544835)); // x = 0.036
218  CHECK(rad_orb[3] == Approx(0.0014678130)); // x = 0.048
219  CHECK(rad_orb[4] == Approx(0.0013891000)); // x = 0.06
220  CHECK(rad_orb[5] == Approx(0.0013175785)); // x = 0.072
221  CHECK(rad_orb[6] == Approx(0.0012523246)); // x = 0.084
222  CHECK(rad_orb[7] == Approx(0.0011923038)); // x = 0.096
223  CHECK(rad_orb[8] == Approx(0.0011364095)); // x = 0.108
224  CHECK(rad_orb[9] == Approx(0.0010837868)); // x = 0.12
225 
226 
227  removeSTypeOrbitals(corrCenter, lcob);
228 
229  CHECK((*lcob.C)(0, 0) == Approx(0.0));
230  CHECK((*lcob.C)(0, 1) == Approx(0.0));
231  CHECK((*lcob.C)(0, 2) == Approx(0.0));
232  CHECK((*lcob.C)(0, 3) != 0.0);
233 }
234 
235 TEST_CASE("HCN MO with cusp", "[wavefunction]")
236 {
237  Communicate* c = OHMMS::Controller;
238 
239  Libxml2Document doc;
240  bool okay = doc.parse("hcn.structure.xml");
241  REQUIRE(okay);
242  xmlNodePtr root = doc.getRoot();
243 
244  const SimulationCell simulation_cell;
245  auto ions_ptr = std::make_unique<ParticleSet>(simulation_cell);
246  auto& ions(*ions_ptr);
247  XMLParticleParser parse_ions(ions);
248  OhmmsXPathObject particleset_ion("//particleset[@name='ion0']", doc.getXPathContext());
249  REQUIRE(particleset_ion.size() == 1);
250  parse_ions.readXML(particleset_ion[0]);
251 
252  REQUIRE(ions.groups() == 3);
253  REQUIRE(ions.R.size() == 3);
254  ions.update();
255 
256  auto elec_ptr = std::make_unique<ParticleSet>(simulation_cell);
257  auto& elec(*elec_ptr);
258  XMLParticleParser parse_elec(elec);
259  OhmmsXPathObject particleset_elec("//particleset[@name='e']", doc.getXPathContext());
260  REQUIRE(particleset_elec.size() == 1);
261  parse_elec.readXML(particleset_elec[0]);
262 
263  REQUIRE(elec.groups() == 2);
264  REQUIRE(elec.R.size() == 14);
265 
266  elec.R = 0.0;
267 
268  elec.addTable(ions);
269  elec.update();
270 
271  Libxml2Document doc2;
272  okay = doc2.parse("hcn.wfnoj.xml");
273  REQUIRE(okay);
274  xmlNodePtr root2 = doc2.getRoot();
275 
276  WaveFunctionComponentBuilder::PSetMap particle_set_map;
277  particle_set_map.emplace(elec_ptr->getName(), std::move(elec_ptr));
278  particle_set_map.emplace(ions_ptr->getName(), std::move(ions_ptr));
279 
280  SPOSetBuilderFactory bf(c, elec, particle_set_map);
281 
282  OhmmsXPathObject MO_base("//determinantset", doc2.getXPathContext());
283  REQUIRE(MO_base.size() == 1);
284 
285  xmlSetProp(MO_base[0], castCharToXMLChar("cuspCorrection"), castCharToXMLChar("yes"));
286 
287  const auto bb_ptr = bf.createSPOSetBuilder(MO_base[0]);
288  auto& bb(*bb_ptr);
289 
290  OhmmsXPathObject slater_base("//determinant", doc2.getXPathContext());
291  auto sposet = bb.createSPOSet(slater_base[0]);
292 
293  SPOSet::ValueVector values;
294  SPOSet::GradVector dpsi;
295  SPOSet::ValueVector d2psi;
296  values.resize(7);
297  dpsi.resize(7);
298  d2psi.resize(7);
299 
300  elec.R = 0.0;
301  elec.update();
302  ParticleSet::SingleParticlePos newpos;
303  elec.makeMove(0, newpos);
304 
305  sposet->evaluateValue(elec, 0, values);
306 
307  // Values from gen_cusp_corr.py
308  CHECK(values[0] == Approx(0.00945227));
309  CHECK(values[1] == Approx(0.0200836));
310  CHECK(values[2] == Approx(0.416375));
311  CHECK(values[3] == Approx(-0.0885443));
312  CHECK(values[4] == Approx(0.273159));
313  CHECK(values[5] == Approx(0));
314  CHECK(values[6] == Approx(0));
315 
316  // Put electron near N atom
317  elec.R[0][0] = -1.09;
318  elec.update();
319  elec.makeMove(0, newpos);
320 
321  values = 0.0;
322  sposet->evaluateValue(elec, 0, values);
323  //std::cout << "values = " << values << std::endl;
324  // Values from gen_cusp_corr.py
325  CHECK(values[0] == Approx(9.5150713253));
326  CHECK(values[1] == Approx(-0.0086731542));
327  CHECK(values[2] == Approx(-1.6426151116));
328  CHECK(values[3] == Approx(0.6569242017));
329  CHECK(values[4] == Approx(0.9775522176));
330  CHECK(values[5] == Approx(0.0000000000));
331  CHECK(values[6] == Approx(0.0000000000));
332 
333 
334  values = 0.0;
335  sposet->evaluateVGL(elec, 0, values, dpsi, d2psi);
336 
337  //std::cout << "values = " << values << std::endl;
338  //std::cout << "dpsi = " << dpsi << std::endl;
339  //std::cout << "d2psi = " << d2psi << std::endl;
340 
341  // Values from gen_cusp_corr.py
342  CHECK(values[0] == Approx(9.5150713253));
343  CHECK(values[1] == Approx(-0.0086731542));
344  CHECK(values[2] == Approx(-1.6426151116));
345  CHECK(values[3] == Approx(0.6569242017));
346  CHECK(values[4] == Approx(0.9775522176));
347  CHECK(values[5] == Approx(0.0000000000));
348  CHECK(values[6] == Approx(0.0000000000));
349 
350  CHECK(dpsi[0][0] == Approx(-66.5007223213));
351  CHECK(dpsi[0][1] == Approx(0.0000000000));
352  CHECK(dpsi[0][2] == Approx(0.0000000000));
353  CHECK(d2psi[0] == Approx(-21540.9990552510));
354 
355  CHECK(values[1] == Approx(-0.0086731542));
356  CHECK(dpsi[1][0] == Approx(0.0616909346));
357  CHECK(dpsi[1][1] == Approx(0.0000000000));
358  CHECK(dpsi[1][2] == Approx(0.0000000000));
359  CHECK(d2psi[1] == Approx(19.8720529007));
360 
361 
362  SPOSet::ValueMatrix all_values;
363  SPOSet::GradMatrix all_grad;
364  SPOSet::ValueMatrix all_lap;
365  all_values.resize(7, 7);
366  all_grad.resize(7, 7);
367  all_lap.resize(7, 7);
368 
369 
370  sposet->evaluate_notranspose(elec, 0, 7, all_values, all_grad, all_lap);
371 
372  // Values from gen_cusp_corr.py
373  CHECK(values[0] == Approx(9.5150713253));
374 
375  CHECK(all_values[0][0] == Approx(9.5150713253));
376  CHECK(all_grad[0][0][0] == Approx(-66.5007223213));
377  CHECK(all_grad[0][0][1] == Approx(0.0000000000));
378  CHECK(all_grad[0][0][2] == Approx(0.0000000000));
379  CHECK(all_lap[0][0] == Approx(-21540.9990552510));
380 
381  CHECK(all_values[0][1] == Approx(-0.0086731542));
382  CHECK(all_grad[0][1][0] == Approx(0.0616909346));
383  CHECK(all_grad[0][1][1] == Approx(0.0000000000));
384  CHECK(all_grad[0][1][2] == Approx(0.0000000000));
385  CHECK(all_lap[0][1] == Approx(19.8720529007));
386 
387 
388  // Test the makeClone method
389  std::unique_ptr<SPOSet> sposet_clone(sposet->makeClone());
390 
391  sposet_clone->evaluate_notranspose(elec, 0, 7, all_values, all_grad, all_lap);
392 
393  // Values from gen_cusp_corr.py
394  CHECK(values[0] == Approx(9.5150713253));
395 
396  CHECK(all_values[0][0] == Approx(9.5150713253));
397  CHECK(all_grad[0][0][0] == Approx(-66.5007223213));
398  CHECK(all_grad[0][0][1] == Approx(0.0000000000));
399  CHECK(all_grad[0][0][2] == Approx(0.0000000000));
400  CHECK(all_lap[0][0] == Approx(-21540.9990552510));
401 
402  CHECK(all_values[0][1] == Approx(-0.0086731542));
403  CHECK(all_grad[0][1][0] == Approx(0.0616909346));
404  CHECK(all_grad[0][1][1] == Approx(0.0000000000));
405  CHECK(all_grad[0][1][2] == Approx(0.0000000000));
406  CHECK(all_lap[0][1] == Approx(19.8720529007));
407 }
408 
409 // Test case with multiple atoms of the same type
410 TEST_CASE("Ethanol MO with cusp", "[wavefunction]")
411 {
412  Communicate* c = OHMMS::Controller;
413 
414  Libxml2Document doc;
415  bool okay = doc.parse("ethanol.structure.xml");
416  REQUIRE(okay);
417  xmlNodePtr root = doc.getRoot();
418 
419  const SimulationCell simulation_cell;
420  auto ions_ptr = std::make_unique<ParticleSet>(simulation_cell);
421  auto& ions(*ions_ptr);
422  XMLParticleParser parse_ions(ions);
423  OhmmsXPathObject particleset_ion("//particleset[@name='ion0']", doc.getXPathContext());
424  REQUIRE(particleset_ion.size() == 1);
425  parse_ions.readXML(particleset_ion[0]);
426 
427  REQUIRE(ions.groups() == 3);
428  REQUIRE(ions.R.size() == 9);
429  ions.update();
430 
431  auto elec_ptr = std::make_unique<ParticleSet>(simulation_cell);
432  auto& elec(*elec_ptr);
433  XMLParticleParser parse_elec(elec);
434  OhmmsXPathObject particleset_elec("//particleset[@name='e']", doc.getXPathContext());
435  REQUIRE(particleset_elec.size() == 1);
436  parse_elec.readXML(particleset_elec[0]);
437 
438  REQUIRE(elec.groups() == 2);
439  REQUIRE(elec.R.size() == 26);
440 
441  elec.R = 0.0;
442 
443  elec.addTable(ions);
444  elec.update();
445 
446  Libxml2Document doc2;
447  okay = doc2.parse("ethanol.wfnoj.xml");
448  REQUIRE(okay);
449  xmlNodePtr root2 = doc2.getRoot();
450 
451  WaveFunctionComponentBuilder::PSetMap particle_set_map;
452  particle_set_map.emplace(elec_ptr->getName(), std::move(elec_ptr));
453  particle_set_map.emplace(ions_ptr->getName(), std::move(ions_ptr));
454 
455  SPOSetBuilderFactory bf(c, elec, particle_set_map);
456 
457  OhmmsXPathObject MO_base("//determinantset", doc2.getXPathContext());
458  REQUIRE(MO_base.size() == 1);
459 
460  xmlSetProp(MO_base[0], castCharToXMLChar("cuspCorrection"), castCharToXMLChar("yes"));
461 
462  const auto bb_ptr = bf.createSPOSetBuilder(MO_base[0]);
463  auto& bb(*bb_ptr);
464 
465  OhmmsXPathObject slater_base("//determinant", doc2.getXPathContext());
466  auto sposet = bb.createSPOSet(slater_base[0]);
467 
468  SPOSet::ValueVector values;
469  SPOSet::GradVector dpsi;
470  SPOSet::ValueVector d2psi;
471  values.resize(13);
472  dpsi.resize(13);
473  d2psi.resize(13);
474 
475  elec.R = 0.0;
476  // Put electron near O atom
477  elec.R[0][0] = -2.10;
478  elec.R[0][1] = 0.50;
479 
480  elec.update();
481  ParticleSet::SingleParticlePos newpos;
482  elec.makeMove(0, newpos);
483 
484  sposet->evaluateValue(elec, 0, values);
485 
486  // Values from gen_cusp_corr.py
487  CHECK(values[0] == Approx(4.3617329704));
488  CHECK(values[1] == Approx(0.0014119853));
489  CHECK(values[2] == Approx(0.0001156461));
490  CHECK(values[3] == Approx(-0.6722670611));
491  CHECK(values[4] == Approx(0.2762949842));
492  CHECK(values[5] == Approx(0.2198735778));
493  CHECK(values[6] == Approx(0.0659454461));
494  CHECK(values[7] == Approx(0.2952071056));
495  CHECK(values[8] == Approx(0.0322071389));
496  CHECK(values[9] == Approx(0.0877981239));
497  CHECK(values[10] == Approx(-0.2151873873));
498  CHECK(values[11] == Approx(0.4250074750));
499  CHECK(values[12] == Approx(0.0767950823));
500 
501  sposet->evaluateVGL(elec, 0, values, dpsi, d2psi);
502 
503  CHECK(values[0] == Approx(4.3617329704));
504  CHECK(values[1] == Approx(0.0014119853));
505  CHECK(values[2] == Approx(0.0001156461));
506  CHECK(values[3] == Approx(-0.6722670611));
507  CHECK(values[4] == Approx(0.2762949842));
508  CHECK(values[5] == Approx(0.2198735778));
509  CHECK(values[6] == Approx(0.0659454461));
510  CHECK(values[7] == Approx(0.2952071056));
511  CHECK(values[8] == Approx(0.0322071389));
512  CHECK(values[9] == Approx(0.0877981239));
513  CHECK(values[10] == Approx(-0.2151873873));
514  CHECK(values[11] == Approx(0.4250074750));
515  CHECK(values[12] == Approx(0.0767950823));
516 
517  CHECK(dpsi[0][0] == Approx(-27.2844138432));
518  CHECK(dpsi[0][1] == Approx(15.9958208598));
519  CHECK(dpsi[0][2] == Approx(0.0195317131));
520  CHECK(d2psi[0] == Approx(-293.2869628790));
521 
522  CHECK(dpsi[12][0] == Approx(1.7548511775));
523  CHECK(dpsi[12][1] == Approx(2.2759333828));
524  CHECK(dpsi[12][2] == Approx(-1.4878277937));
525  CHECK(d2psi[12] == Approx(-4.3399821309));
526 
527 
528  SPOSet::ValueMatrix all_values;
529  SPOSet::GradMatrix all_grad;
530  SPOSet::ValueMatrix all_lap;
531  all_values.resize(10, 13);
532  all_grad.resize(10, 13);
533  all_lap.resize(10, 13);
534 
535  sposet->evaluate_notranspose(elec, 0, 7, all_values, all_grad, all_lap);
536 
537  CHECK(all_values[0][0] == Approx(4.3617329704));
538  CHECK(all_grad[0][0][0] == Approx(-27.2844138432));
539  CHECK(all_grad[0][0][1] == Approx(15.9958208598));
540  CHECK(all_grad[0][0][2] == Approx(0.0195317131));
541  CHECK(all_lap[0][0] == Approx(-293.2869628790));
542 
543  CHECK(all_values[0][11] == Approx(0.4250074750));
544  CHECK(all_grad[0][11][0] == Approx(-0.3947036210));
545  CHECK(all_grad[0][11][1] == Approx(0.9883840215));
546  CHECK(all_grad[0][11][2] == Approx(1.7863218842));
547  CHECK(all_lap[0][11] == Approx(-33.5202249813));
548 }
549 
550 
551 TEST_CASE("broadcastCuspInfo", "[wavefunction]")
552 {
553  Communicate* c = OHMMS::Controller;
554  CuspCorrectionParameters cp;
555  int root = 0;
556  if (c->rank() == root)
557  {
558  cp.Rc = 2.0;
559  cp.C = 3.0;
560  cp.sg = -1.0;
561  cp.alpha[0] = 1.1;
562  cp.alpha[1] = 1.2;
563  cp.alpha[2] = 1.3;
564  cp.alpha[3] = 1.4;
565  cp.alpha[4] = 1.5;
566  cp.redo = 1;
567  }
568 
569  broadcastCuspInfo(cp, *c, root);
570 
571  CHECK(cp.Rc == Approx(2.0));
572  CHECK(cp.C == Approx(3.0));
573  CHECK(cp.sg == Approx(-1.0));
574  CHECK(cp.alpha[0] == Approx(1.1));
575  CHECK(cp.alpha[1] == Approx(1.2));
576  CHECK(cp.alpha[2] == Approx(1.3));
577  CHECK(cp.alpha[3] == Approx(1.4));
578  CHECK(cp.alpha[4] == Approx(1.5));
579  REQUIRE(cp.redo == 1);
580 }
581 
582 
583 } // namespace qmcplusplus
qmcplusplus::Vector::resize
void resize(size_type n, Type_t val=Type_t())
Resize the container.
Definition: OhmmsVector.h:166
qmcplusplus::TinyVector
Fixed-size array.
Definition: OhmmsTinyMeta.h:30
Libxml2Document
class that handles xmlDoc
Definition: Libxml2Doc.h:76
qmcplusplus::CuspCorrectionParameters::alpha
TinyVector< ValueType, 5 > alpha
The coefficients of the polynomial in Eq 8.
Definition: CuspCorrection.h:52
qmcplusplus::LinearGrid
One-Dimensional linear-grid.
Definition: OneDimGridBase.h:134
qmcplusplus
helper functions for EinsplineSetBuilder
Definition: Configuration.h:43
Communicate::rank
int rank() const
return the rank
Definition: Communicate.h:116
qmcplusplus::QMCTraits::RealType
QTBase::RealType RealType
Definition: Configuration.h:58
qmcplusplus::CuspCorrectionParameters
Cusp correction parameters.
Definition: CuspCorrection.h:37
Libxml2Document::getXPathContext
xmlXPathContextPtr getXPathContext()
Definition: Libxml2Doc.cpp:100
qmcplusplus::CuspCorrectionParameters::sg
RealType sg
The sign of the wavefunction at the nucleus.
Definition: CuspCorrection.h:49
qmcplusplus::CuspCorrectionParameters::Rc
RealType Rc
The cutoff radius.
Definition: CuspCorrection.h:43
qmcplusplus::TEST_CASE
TEST_CASE("complex_helper", "[type_traits]")
Definition: test_complex_helper.cpp:38
Libxml2Document::getRoot
xmlNodePtr getRoot()
Definition: Libxml2Doc.h:88
qmcplusplus::LCAOrbitalSet::myBasisSet
std::unique_ptr< basis_type > myBasisSet
pointer to the basis set
Definition: LCAOrbitalSet.h:40
qmcplusplus::SPOSet::ValueMatrix
OrbitalSetTraits< ValueType >::ValueMatrix ValueMatrix
Definition: SPOSet.h:50
qmcplusplus::Matrix::resize
void resize(size_type n, size_type m)
Resize the container.
Definition: OhmmsMatrix.h:99
OHMMS::Controller
Communicate * Controller
Global Communicator for a process.
Definition: Communicate.cpp:35
qmcplusplus::broadcastCuspInfo
void broadcastCuspInfo(CuspCorrectionParameters &param, Communicate &Comm, int root)
Broadcast cusp correction parameters.
Definition: CuspCorrectionConstruction.cpp:198
qmcplusplus::WaveFunctionComponentBuilder::PSetMap
std::map< std::string, const std::unique_ptr< ParticleSet > > PSetMap
Definition: WaveFunctionComponentBuilder.h:43
OhmmsXPathObject
class to handle xmlXPathObject
Definition: Libxml2Doc.h:26
qmcplusplus::CuspCorrectionParameters::redo
int redo
Flag to indicate the correction should be recalculated.
Definition: CuspCorrection.h:55
qmcplusplus::OneDimGridBase
An abstract base class to implement a One-Dimensional grid.
Definition: OneDimGridBase.h:35
qmcplusplus::readCuspInfo
bool readCuspInfo(const std::string &cuspInfoFile, const std::string &objectName, int OrbitalSetSize, Matrix< CuspCorrectionParameters > &info)
Read cusp correction parameters from XML file.
Definition: CuspCorrectionConstruction.cpp:30
qmcplusplus::SPOSet::GradMatrix
OrbitalSetTraits< ValueType >::GradMatrix GradMatrix
Definition: SPOSet.h:52
Communicate
Wrapping information on parallelism.
Definition: Communicate.h:68
OneDimGridBase.h
Decalaration of One-Dimesional grids.
qmcplusplus::removeSTypeOrbitals
void removeSTypeOrbitals(const std::vector< bool > &corrCenter, LCAOrbitalSet &Phi)
Remove S atomic orbitals from all molecular orbitals on all centers.
Definition: CuspCorrectionConstruction.cpp:256
qmcplusplus::REQUIRE
REQUIRE(std::filesystem::exists(filename))
qmcplusplus::SPOSet::ValueVector
OrbitalSetTraits< ValueType >::ValueVector ValueVector
Definition: SPOSet.h:49
qmcplusplus::LCAOrbitalSet::C
std::shared_ptr< OffloadValueMatrix > C
pointer to matrix containing the coefficients
Definition: LCAOrbitalSet.h:42
qmcplusplus::SPOSet::getOrbitalSetSize
int getOrbitalSetSize() const
return the size of the orbitals
Definition: SPOSet.h:88
qmcplusplus::XMLParticleParser::readXML
bool readXML(xmlNodePtr cur)
process xmlnode <particleset/> which contains everything about the particle set to initialize ...
Definition: XMLParticleIO.cpp:90
qmcplusplus::computeRadialPhiBar
void computeRadialPhiBar(ParticleSet *targetP, ParticleSet *sourceP, int curOrb_, int curCenter_, SPOSet *Phi, Vector< QMCTraits::RealType > &xgrid, Vector< QMCTraits::RealType > &rad_orb, const CuspCorrectionParameters &data)
Compute the radial part of the corrected wavefunction.
Definition: CuspCorrectionConstruction.cpp:280
qmcplusplus::SPOSet::GradVector
OrbitalSetTraits< ValueType >::GradVector GradVector
Definition: SPOSet.h:51
qmcplusplus::CHECK
CHECK(log_values[0]==ComplexApprox(std::complex< double >{ 5.603777579195571, -6.1586603331188225 }))
Libxml2Document::parse
bool parse(const std::string &fname)
Definition: Libxml2Doc.cpp:180
qmcplusplus::LCAOrbitalSet
class to handle linear combinations of basis orbitals used to evaluate the Dirac determinants.
Definition: LCAOrbitalSet.h:30
qmcplusplus::SPOSetBuilderFactory::createSPOSetBuilder
std::unique_ptr< SPOSetBuilder > createSPOSetBuilder(xmlNodePtr rootNode)
Definition: SPOSetBuilderFactory.cpp:71
qmcplusplus::CuspCorrectionParameters::C
RealType C
A shift to keep correction to a single sign.
Definition: CuspCorrection.h:46
qmcplusplus::LCAOrbitalSet::isOMPoffload
bool isOMPoffload() const override
Query if this SPOSet uses OpenMP offload.
Definition: LCAOrbitalSet.h:58
castCharToXMLChar
xmlChar * castCharToXMLChar(char *c)
unsafe char* to xmlChar* cast
Definition: libxmldefs.h:70
qmcplusplus::splitPhiEta
void splitPhiEta(int center, const std::vector< bool > &corrCenter, LCAOrbitalSet &Phi, LCAOrbitalSet &Eta)
Divide molecular orbital into atomic S-orbitals on this center (phi), and everything else (eta)...
Definition: CuspCorrectionConstruction.cpp:226