QMCPACK
test_MO.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) 2018 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory
8 //
9 // File created by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory
10 //////////////////////////////////////////////////////////////////////////////////////
11 
12 
13 #include "catch.hpp"
14 
15 #include "Configuration.h"
16 #include "Message/Communicate.h"
17 #include "Numerics/OneDimGridBase.h"
18 #include "ParticleIO/XMLParticleIO.h"
19 #include "Numerics/GaussianBasisSet.h"
20 #include "QMCWaveFunctions/LCAO/LCAOrbitalBuilder.h"
21 #include "QMCWaveFunctions/SPOSetBuilderFactory.h"
22 #include <ResourceCollection.h>
23 
24 namespace qmcplusplus
25 {
26 void test_He(bool transform)
27 {
28  std::ostringstream section_name;
29  section_name << "He, transform orbitals to grid: " << (transform ? "T" : "F");
30 
31  SECTION(section_name.str())
32  {
33  Communicate* c = OHMMS::Controller;
34 
35  const SimulationCell simulation_cell;
36  auto elec_ptr = std::make_unique<ParticleSet>(simulation_cell);
37  auto& elec(*elec_ptr);
38  std::vector<int> agroup(2);
39  agroup[0] = 1;
40  agroup[1] = 1;
41  elec.setName("e");
42  elec.create(agroup);
43  elec.R[0] = 0.0;
44 
45  SpeciesSet& tspecies = elec.getSpeciesSet();
46  int upIdx = tspecies.addSpecies("u");
47  int downIdx = tspecies.addSpecies("d");
48  int massIdx = tspecies.addAttribute("mass");
49  tspecies(massIdx, upIdx) = 1.0;
50  tspecies(massIdx, downIdx) = 1.0;
51 
52  auto ions_ptr = std::make_unique<ParticleSet>(simulation_cell);
53  auto& ions(*ions_ptr);
54  ions.setName("ion0");
55  ions.create({1});
56  ions.R[0] = 0.0;
57  SpeciesSet& ispecies = ions.getSpeciesSet();
58  int heIdx = ispecies.addSpecies("He");
59  ions.update();
60 
61  elec.addTable(ions);
62  elec.update();
63 
64  Libxml2Document doc;
65  bool okay = doc.parse("he_sto3g.wfj.xml");
66  REQUIRE(okay);
67  xmlNodePtr root = doc.getRoot();
68 
69  WaveFunctionComponentBuilder::PSetMap particle_set_map;
70  particle_set_map.emplace(elec_ptr->getName(), std::move(elec_ptr));
71  particle_set_map.emplace(ions_ptr->getName(), std::move(ions_ptr));
72 
73  SPOSetBuilderFactory bf(c, elec, particle_set_map);
74 
75  OhmmsXPathObject MO_base("//determinantset", doc.getXPathContext());
76  REQUIRE(MO_base.size() == 1);
77  if (!transform)
78  {
79  // input file is set to transform GTO's to numerical orbitals by default
80  // use direct evaluation of GTO's
81  xmlSetProp(MO_base[0], castCharToXMLChar("transform"), castCharToXMLChar("no"));
82  xmlSetProp(MO_base[0], castCharToXMLChar("key"), castCharToXMLChar("GTO"));
83  }
84 
85  const auto bb_ptr = bf.createSPOSetBuilder(MO_base[0]);
86  auto& bb(*bb_ptr);
87 
88  OhmmsXPathObject slater_base("//determinant", doc.getXPathContext());
89  auto sposet = bb.createSPOSet(slater_base[0]);
90 
91  //std::cout << "basis set size = " << sposet->getBasisSetSize() << std::endl;
92 
93  SPOSet::ValueVector values;
94  SPOSet::GradVector dpsi;
95  SPOSet::ValueVector d2psi;
96  values.resize(1);
97  dpsi.resize(1);
98  d2psi.resize(1);
99 
100  // Call makeMove to compute the distances
101  ParticleSet::SingleParticlePos newpos(0.0001, 0.0, 0.0);
102  elec.makeMove(0, newpos);
103 
104  sposet->evaluateValue(elec, 0, values);
105 
106  // Generated from gen_mo.py for position [0.0001, 0.0, 0.0]
107  CHECK(values[0] == Approx(0.9996037001));
108 
109  sposet->evaluateVGL(elec, 0, values, dpsi, d2psi);
110 
111  // Generated from gen_mo.py for position [0.0001, 0.0, 0.0]
112  CHECK(values[0] == Approx(0.9996037001));
113  CHECK(dpsi[0][0] == Approx(-0.0006678035459));
114  CHECK(dpsi[0][1] == Approx(0));
115  CHECK(dpsi[0][2] == Approx(0));
116  CHECK(d2psi[0] == Approx(-20.03410564));
117 
118 
119  ParticleSet::SingleParticlePos disp(1.0, 0.0, 0.0);
120  elec.makeMove(0, disp);
121 
122  sposet->evaluateVGL(elec, 0, values, dpsi, d2psi);
123  // Generated from gen_mo.py for position [1.0, 0.0, 0.0]
124  CHECK(values[0] == Approx(0.2315567641));
125  CHECK(dpsi[0][0] == Approx(-0.3805431885));
126  CHECK(dpsi[0][1] == Approx(0));
127  CHECK(dpsi[0][2] == Approx(0));
128  CHECK(d2psi[0] == Approx(-0.2618497452));
129  }
130 }
131 
132 TEST_CASE("ReadMolecularOrbital GTO He", "[wavefunction]") { test_He(false); }
133 
134 TEST_CASE("ReadMolecularOrbital Numerical He", "[wavefunction]") { test_He(true); }
135 
136 void test_He_mw(bool transform)
137 {
138  // set up ion particle set as normal
139  Communicate* c = OHMMS::Controller;
140 
141  const SimulationCell simulation_cell;
142  auto elec_ptr = std::make_unique<ParticleSet>(simulation_cell);
143  auto& elec(*elec_ptr);
144  std::vector<int> agroup(2);
145  agroup[0] = 1;
146  agroup[1] = 1;
147  elec.setName("e");
148  elec.create(agroup);
149  elec.R[0] = 0.0;
150 
151  SpeciesSet& tspecies = elec.getSpeciesSet();
152  int upIdx = tspecies.addSpecies("u");
153  int downIdx = tspecies.addSpecies("d");
154  int massIdx = tspecies.addAttribute("mass");
155  tspecies(massIdx, upIdx) = 1.0;
156  tspecies(massIdx, downIdx) = 1.0;
157 
158  auto ions_ptr = std::make_unique<ParticleSet>(simulation_cell);
159  auto& ions(*ions_ptr);
160  ions.setName("ion0");
161  ions.create({1});
162  ions.R[0] = 0.0;
163  SpeciesSet& ispecies = ions.getSpeciesSet();
164  int heIdx = ispecies.addSpecies("He");
165  ions.update();
166 
167  elec.addTable(ions);
168  elec.update();
169 
170  Libxml2Document doc;
171  bool okay = doc.parse("he_sto3g.wfj.xml");
172  REQUIRE(okay);
173  xmlNodePtr root = doc.getRoot();
174 
175  WaveFunctionComponentBuilder::PSetMap particle_set_map;
176  particle_set_map.emplace(elec_ptr->getName(), std::move(elec_ptr));
177  particle_set_map.emplace(ions_ptr->getName(), std::move(ions_ptr));
178 
179  SPOSetBuilderFactory bf(c, elec, particle_set_map);
180 
181  OhmmsXPathObject MO_base("//determinantset", doc.getXPathContext());
182  REQUIRE(MO_base.size() == 1);
183  if (!transform)
184  {
185  // input file is set to transform GTO's to numerical orbitals by default
186  // use direct evaluation of GTO's
187  xmlSetProp(MO_base[0], castCharToXMLChar("transform"), castCharToXMLChar("no"));
188  xmlSetProp(MO_base[0], castCharToXMLChar("key"), castCharToXMLChar("GTO"));
189  }
190 
191  const auto bb_ptr = bf.createSPOSetBuilder(MO_base[0]);
192  auto& bb(*bb_ptr);
193 
194  OhmmsXPathObject slater_base("//determinant", doc.getXPathContext());
195  auto sposet = bb.createSPOSet(slater_base[0]);
196 
197  //std::cout << "basis set size = " << sposet->getBasisSetSize() << std::endl;
198 
199  SPOSet::ValueVector psi;
200  SPOSet::GradVector dpsi;
201  SPOSet::ValueVector d2psi;
202  psi.resize(1);
203  dpsi.resize(1);
204  d2psi.resize(1);
205 
206  // Call makeMove to compute the distances
207  ParticleSet::SingleParticlePos newpos(0.0001, 0.0, 0.0);
208  elec.makeMove(0, newpos);
209  // set up second walkers
210  // auto elec2 = elec.makeClone();
211 
212  sposet->evaluateVGL(elec, 0, psi, dpsi, d2psi);
213  CHECK(std::real(psi[0]) == Approx(0.9996037001));
214  CHECK(std::real(dpsi[0][0]) == Approx(-0.000667803579));
215  CHECK(std::real(dpsi[0][1]) == Approx(0));
216  CHECK(std::real(dpsi[0][2]) == Approx(0));
217  CHECK(std::real(d2psi[0]) == Approx(-20.0342132));
218 
219 
220  // vectors of SPOSets, ParticleSets, V/G/L (leading dim of each == nwalkers)
221  RefVectorWithLeader<SPOSet> spo_list(*sposet);
222  spo_list.push_back(*sposet);
223 
224  RefVectorWithLeader<ParticleSet> P_list(elec);
225  P_list.push_back(elec);
226 
227  RefVector<SPOSet::ValueVector> psi_list;
228  RefVector<SPOSet::GradVector> dpsi_list;
229  RefVector<SPOSet::ValueVector> d2psi_list;
230 
231  // create V,G,L arrays for walker 1
232  SPOSet::ValueVector psi_1(sposet->getOrbitalSetSize());
233  SPOSet::GradVector dpsi_1(sposet->getOrbitalSetSize());
234  SPOSet::ValueVector d2psi_1(sposet->getOrbitalSetSize());
235 
236  psi_list.push_back(psi_1);
237  dpsi_list.push_back(dpsi_1);
238  d2psi_list.push_back(d2psi_1);
239 
240  //second walker
241  // interchange positions and shift y instead of x
242  ParticleSet::SingleParticlePos dy(0.0, 0.0001, 0.0);
243  ParticleSet elec_2(elec);
244  elec_2.R[0] = elec.R[1];
245  elec_2.R[1] = elec.R[0];
246  elec_2.update();
247  elec_2.makeMove(0, dy);
248 
249  std::unique_ptr<SPOSet> sposet_2(sposet->makeClone());
250  SPOSet::ValueVector psi_2(sposet->getOrbitalSetSize());
251  SPOSet::GradVector dpsi_2(sposet->getOrbitalSetSize());
252  SPOSet::ValueVector d2psi_2(sposet->getOrbitalSetSize());
253  spo_list.push_back(*sposet_2);
254  P_list.push_back(elec_2);
255  psi_list.push_back(psi_2);
256  dpsi_list.push_back(dpsi_2);
257  d2psi_list.push_back(d2psi_2);
258 
259  ResourceCollection pset_res("test_pset_res");
260  ResourceCollection spo_res("test_spo_res");
261 
262  elec.createResource(pset_res);
263  sposet->createResource(spo_res);
264 
265  ResourceCollectionTeamLock<ParticleSet> mw_pset_lock(pset_res, P_list);
266  ResourceCollectionTeamLock<SPOSet> mw_sposet_lock(spo_res, spo_list);
267 
268  sposet->mw_evaluateVGL(spo_list, P_list, 0, psi_list, dpsi_list, d2psi_list);
269 
270  CHECK(std::real(psi_list[0].get()[0]) == Approx(psi[0]));
271  CHECK(std::real(dpsi_list[0].get()[0][0]) == Approx(dpsi[0][0]));
272  CHECK(std::real(dpsi_list[0].get()[0][1]) == Approx(dpsi[0][1]));
273  CHECK(std::real(dpsi_list[0].get()[0][2]) == Approx(dpsi[0][2]));
274  CHECK(std::real(d2psi_list[0].get()[0]) == Approx(d2psi[0]));
275 
276  CHECK(std::real(psi_list[1].get()[0]) == Approx(psi[0]));
277  CHECK(std::real(dpsi_list[1].get()[0][0]) == Approx(dpsi[0][1])); // x, y switched here
278  CHECK(std::real(dpsi_list[1].get()[0][1]) == Approx(dpsi[0][0]));
279  CHECK(std::real(dpsi_list[1].get()[0][2]) == Approx(dpsi[0][2]));
280  CHECK(std::real(d2psi_list[1].get()[0]) == Approx(d2psi[0]));
281 }
282 
283 TEST_CASE("mw_evaluate Numerical He", "[wavefunction]") { test_He_mw(true); }
284 
285 void test_EtOH_mw(bool transform)
286 {
287  // set up ion particle set as normal
288  Communicate* c = OHMMS::Controller;
289 
290  Libxml2Document doc;
291  bool okay = doc.parse("ethanol.structure.xml");
292  REQUIRE(okay);
293  xmlNodePtr root = doc.getRoot();
294 
295  const SimulationCell simulation_cell;
296  auto ions_ptr = std::make_unique<ParticleSet>(simulation_cell);
297  auto& ions(*ions_ptr);
298  XMLParticleParser parse_ions(ions);
299  OhmmsXPathObject particleset_ion("//particleset[@name='ion0']", doc.getXPathContext());
300  REQUIRE(particleset_ion.size() == 1);
301  parse_ions.readXML(particleset_ion[0]);
302 
303  REQUIRE(ions.groups() == 3);
304  REQUIRE(ions.R.size() == 9);
305  ions.update();
306 
307  auto elec_ptr = std::make_unique<ParticleSet>(simulation_cell);
308  auto& elec(*elec_ptr);
309  XMLParticleParser parse_elec(elec);
310  OhmmsXPathObject particleset_elec("//particleset[@name='e']", doc.getXPathContext());
311  REQUIRE(particleset_elec.size() == 1);
312  parse_elec.readXML(particleset_elec[0]);
313 
314  REQUIRE(elec.groups() == 2);
315  REQUIRE(elec.R.size() == 26);
316 
317  elec.R = 0.0;
318 
319  elec.addTable(ions);
320  elec.update();
321 
322  Libxml2Document doc2;
323  okay = doc2.parse("ethanol.wfnoj.xml");
324  REQUIRE(okay);
325  xmlNodePtr root2 = doc2.getRoot();
326 
327  WaveFunctionComponentBuilder::PSetMap particle_set_map;
328  particle_set_map.emplace(elec_ptr->getName(), std::move(elec_ptr));
329  particle_set_map.emplace(ions_ptr->getName(), std::move(ions_ptr));
330 
331  SPOSetBuilderFactory bf(c, elec, particle_set_map);
332 
333  OhmmsXPathObject MO_base("//determinantset", doc2.getXPathContext());
334  REQUIRE(MO_base.size() == 1);
335  if (!transform)
336  {
337  // input file is set to transform GTO's to numerical orbitals by default
338  // use direct evaluation of GTO's
339  xmlSetProp(MO_base[0], castCharToXMLChar("transform"), castCharToXMLChar("no"));
340  xmlSetProp(MO_base[0], castCharToXMLChar("key"), castCharToXMLChar("GTO"));
341  }
342 
343  xmlSetProp(MO_base[0], castCharToXMLChar("cuspCorrection"), castCharToXMLChar("no"));
344 
345  const auto bb_ptr = bf.createSPOSetBuilder(MO_base[0]);
346  auto& bb(*bb_ptr);
347 
348  OhmmsXPathObject slater_base("//determinant", doc2.getXPathContext());
349  auto sposet = bb.createSPOSet(slater_base[0]);
350 
351 
352  //std::cout << "basis set size = " << sposet->getBasisSetSize() << std::endl;
353  size_t n_mo = sposet->getOrbitalSetSize();
354  SPOSet::ValueVector psiref_0(n_mo);
355  SPOSet::GradVector dpsiref_0(n_mo);
356  SPOSet::ValueVector d2psiref_0(n_mo);
357  // Call makeMove to compute the distances
358  //ParticleSet::SingleParticlePos newpos(0.0001, 0.0, 0.0);
359  //elec.makeMove(0, newpos);
360  //elec.update();
361  elec.R[0] = {0.0001, 0.0, 0.0};
362  elec.update();
363  // set up second walkers
364  // auto elec2 = elec.makeClone();
365  sposet->evaluateVGL(elec, 0, psiref_0, dpsiref_0, d2psiref_0);
366 
367  CHECK(std::real(psiref_0[0]) == Approx(-0.001664403313));
368  CHECK(std::real(psiref_0[1]) == Approx(0.01579976715));
369  CHECK(std::real(dpsiref_0[0][0]) == Approx(-0.0001961749098));
370  CHECK(std::real(dpsiref_0[0][1]) == Approx(0.003340392074));
371  CHECK(std::real(dpsiref_0[0][2]) == Approx(9.461877818e-05));
372  CHECK(std::real(dpsiref_0[1][0]) == Approx(-0.005476152264));
373  CHECK(std::real(dpsiref_0[1][1]) == Approx(-0.06648077046));
374  CHECK(std::real(dpsiref_0[1][2]) == Approx(2.086541402e-05));
375  CHECK(std::real(d2psiref_0[0]) == Approx(0.01071299243));
376  CHECK(std::real(d2psiref_0[1]) == Approx(0.4121970776));
377 
378  //ParticleSet::SingleParticlePos newpos2(0.0, 0.04, 0.02);
379  //elec.makeMove(1, newpos2);
380  elec.R[1] = {0.0, 0.04, 0.02};
381  elec.update();
382  SPOSet::ValueVector psiref_1(n_mo);
383  SPOSet::GradVector dpsiref_1(n_mo);
384  SPOSet::ValueVector d2psiref_1(n_mo);
385  sposet->evaluateVGL(elec, 1, psiref_1, dpsiref_1, d2psiref_1);
386 
387  CHECK(std::real(psiref_1[0]) == Approx(-0.001528135727));
388  CHECK(std::real(psiref_1[0]) == Approx(-0.001528135727));
389  CHECK(std::real(psiref_1[1]) == Approx(0.01351541907));
390  CHECK(std::real(d2psiref_1[0]) == Approx(0.01001796854));
391  CHECK(std::real(d2psiref_1[1]) == Approx(0.2912963205));
392  CHECK(std::real(dpsiref_1[0][0]) == Approx(-0.0004235196101));
393  CHECK(std::real(dpsiref_1[0][1]) == Approx(0.003351193375));
394  CHECK(std::real(dpsiref_1[0][2]) == Approx(0.0001374796409));
395  CHECK(std::real(dpsiref_1[1][0]) == Approx(-0.003873067027));
396  CHECK(std::real(dpsiref_1[1][1]) == Approx(-0.0483167767));
397  CHECK(std::real(dpsiref_1[1][2]) == Approx(-0.0008320732335));
398 
399  // vectors of SPOSets, ParticleSets, V/G/L (leading dim of each == nwalkers)
400  RefVectorWithLeader<SPOSet> spo_list(*sposet);
401  std::unique_ptr<SPOSet> sposet_2(sposet->makeClone());
402  spo_list.push_back(*sposet_2);
403  spo_list.push_back(*sposet);
404  // second walker
405  // exchange positions
406  ParticleSet elec_2(elec);
407  elec_2.R[0] = elec.R[1];
408  elec_2.R[1] = elec.R[0];
409  elec_2.update();
410  RefVectorWithLeader<ParticleSet> P_list(elec);
411  P_list.push_back(elec);
412  P_list.push_back(elec_2);
413  // create V,G,L arrays for walker 1
414  SPOSet::ValueVector psi_1(n_mo);
415  SPOSet::GradVector dpsi_1(n_mo);
416  SPOSet::ValueVector d2psi_1(n_mo);
417  SPOSet::ValueVector psi_2(n_mo);
418  SPOSet::GradVector dpsi_2(n_mo);
419  SPOSet::ValueVector d2psi_2(n_mo);
420  RefVector<SPOSet::ValueVector> psi_list = {psi_1, psi_2};
421  RefVector<SPOSet::GradVector> dpsi_list = {dpsi_1, dpsi_2};
422  RefVector<SPOSet::ValueVector> d2psi_list = {d2psi_1, d2psi_2};
423 
424  size_t nw = psi_list.size();
425  SPOSet::ValueVector psi_v_1(n_mo);
426  SPOSet::ValueVector psi_v_2(n_mo);
427  RefVector<SPOSet::ValueVector> psi_v_list{psi_v_1, psi_v_2};
428 
429  ResourceCollection pset_res("test_pset_res");
430  ResourceCollection spo_res("test_spo_res");
431 
432  elec.createResource(pset_res);
433  sposet->createResource(spo_res);
434 
435  ResourceCollectionTeamLock<ParticleSet> mw_pset_lock(pset_res, P_list);
436  ResourceCollectionTeamLock<SPOSet> mw_sposet_lock(spo_res, spo_list);
437 
438  sposet->mw_evaluateVGL(spo_list, P_list, 0, psi_list, dpsi_list, d2psi_list);
439  sposet->mw_evaluateValue(spo_list, P_list, 0, psi_v_list);
440 
441  for (size_t iorb = 0; iorb < n_mo; iorb++)
442  {
443  for (size_t iw = 0; iw < nw; iw++)
444  {
445  // test values from OffloadMWVArray impl.
446  CHECK(std::real(psi_v_list[iw].get()[iorb]) == Approx(psi_list[iw].get()[iorb]));
447  }
448  CHECK(std::real(psi_v_list[0].get()[iorb]) == Approx(psiref_0[iorb]));
449  CHECK(std::real(psi_v_list[1].get()[iorb]) == Approx(psiref_1[iorb]));
450  CHECK(std::real(psi_list[0].get()[iorb]) == Approx(psiref_0[iorb]));
451  CHECK(std::real(psi_list[1].get()[iorb]) == Approx(psiref_1[iorb]));
452  CHECK(std::real(d2psi_list[0].get()[iorb]) == Approx(d2psiref_0[iorb]));
453  CHECK(std::real(d2psi_list[1].get()[iorb]) == Approx(d2psiref_1[iorb]));
454  for (size_t idim = 0; idim < SPOSet::DIM; idim++)
455  {
456  CHECK(std::real(dpsi_list[0].get()[iorb][idim]) == Approx(dpsiref_0[iorb][idim]));
457  CHECK(std::real(dpsi_list[1].get()[iorb][idim]) == Approx(dpsiref_1[iorb][idim]));
458  }
459  }
460 }
461 
462 TEST_CASE("mw_evaluate Numerical EtOH", "[wavefunction]") { test_EtOH_mw(true); }
463 TEST_CASE("mw_evaluate GTO EtOH", "[wavefunction]") { test_EtOH_mw(false); }
464 
465 void test_Ne(bool transform)
466 {
467  std::ostringstream section_name;
468  section_name << "Neon, transform orbitals to grid: " << (transform ? "T" : "F");
469 
470  SECTION(section_name.str())
471  {
472  Communicate* c = OHMMS::Controller;
473 
474  const SimulationCell simulation_cell;
475  auto elec_ptr = std::make_unique<ParticleSet>(simulation_cell);
476  auto& elec(*elec_ptr);
477 
478  std::vector<int> agroup(2);
479  agroup[0] = 1;
480  agroup[1] = 1;
481  elec.setName("e");
482  elec.create(agroup);
483  elec.R[0] = 0.0;
484 
485  SpeciesSet& tspecies = elec.getSpeciesSet();
486  int upIdx = tspecies.addSpecies("u");
487  int downIdx = tspecies.addSpecies("d");
488  int massIdx = tspecies.addAttribute("mass");
489  tspecies(massIdx, upIdx) = 1.0;
490  tspecies(massIdx, downIdx) = 1.0;
491 
492  auto ions_ptr = std::make_unique<ParticleSet>(simulation_cell);
493  auto& ions(*ions_ptr);
494  ions.setName("ion0");
495  ions.create({1});
496  ions.R[0] = {0.0, 0.0, 0.0};
497  SpeciesSet& ispecies = ions.getSpeciesSet();
498  int heIdx = ispecies.addSpecies("Ne");
499  ions.update();
500 
501 
502  elec.addTable(ions);
503  elec.update();
504 
505  Libxml2Document doc;
506  bool okay = doc.parse("ne_def2_svp.wfnoj.xml");
507  REQUIRE(okay);
508  xmlNodePtr root = doc.getRoot();
509 
510  WaveFunctionComponentBuilder::PSetMap particle_set_map;
511  particle_set_map.emplace(elec_ptr->getName(), std::move(elec_ptr));
512  particle_set_map.emplace(ions_ptr->getName(), std::move(ions_ptr));
513 
514  SPOSetBuilderFactory bf(c, elec, particle_set_map);
515 
516  OhmmsXPathObject MO_base("//determinantset", doc.getXPathContext());
517  REQUIRE(MO_base.size() == 1);
518 
519  if (!transform)
520  {
521  // input file is set to transform GTO's to numerical orbitals by default
522  // use direct evaluation of GTO's
523  xmlSetProp(MO_base[0], castCharToXMLChar("transform"), castCharToXMLChar("no"));
524  xmlSetProp(MO_base[0], castCharToXMLChar("key"), castCharToXMLChar("GTO"));
525  }
526 
527  const auto bb_ptr = bf.createSPOSetBuilder(MO_base[0]);
528  auto& bb(*bb_ptr);
529 
530  OhmmsXPathObject slater_base("//determinant", doc.getXPathContext());
531  auto sposet = bb.createSPOSet(slater_base[0]);
532 
533  //std::cout << "basis set size = " << sposet->getBasisSetSize() << std::endl;
534 
535  const int norbs = 5;
536  SPOSet::ValueVector values;
537  SPOSet::GradVector dpsi;
538  SPOSet::ValueVector d2psi;
539  values.resize(norbs);
540  dpsi.resize(norbs);
541  d2psi.resize(norbs);
542 
543  ParticleSet::SingleParticlePos newpos(0.00001, 0.0, 0.0);
544  elec.makeMove(0, newpos);
545 
546  sposet->evaluateValue(elec, 0, values);
547 
548  // Generated from gen_mo.py for position [1e-05, 0.0, 0.0]
549  CHECK(values[0] == Approx(-16.11819042));
550 
551  sposet->evaluateVGL(elec, 0, values, dpsi, d2psi);
552  // Generated from gen_mo.py for position [1e-05, 0.0, 0.0]
553  CHECK(values[0] == Approx(-16.11819042));
554  CHECK(dpsi[0][0] == Approx(0.1747261458));
555  CHECK(dpsi[0][1] == Approx(0));
556  CHECK(dpsi[0][2] == Approx(0));
557 
558 
559  ParticleSet::SingleParticlePos disp(1.0, 0.0, 0.0);
560  elec.makeMove(0, disp);
561  sposet->evaluateValue(elec, 0, values);
562  // Generated from gen_mo.py for position [1.0, 0.0, 0.0]
563  CHECK(values[0] == Approx(-0.005041631374));
564 
565 
566  sposet->evaluateVGL(elec, 0, values, dpsi, d2psi);
567  // Generated from gen_mo.py for position [1.0, 0.0, 0.0]
568  CHECK(values[0] == Approx(-0.005041631374));
569  CHECK(dpsi[0][0] == Approx(0.01862216578));
570  CHECK(dpsi[0][1] == Approx(0));
571  CHECK(dpsi[0][2] == Approx(0));
572  CHECK(d2psi[0] == Approx(-0.01551755818));
573 
574  // when a determinant only contains a single particle.
575  SPOSet::ValueVector phi(1), phiinv(1);
576  phiinv[0] = 100;
577  VirtualParticleSet VP(elec, 2);
578  std::vector<ParticleSet::SingleParticlePos> newpos2(2);
579  std::vector<SPOSet::ValueType> ratios2(2);
580  newpos2[0] = disp;
581  newpos2[1] = -disp;
582  VP.makeMoves(elec, 0, newpos2);
583  sposet->evaluateDetRatios(VP, phi, phiinv, ratios2);
584  CHECK(ratios2[0] == Approx(-0.504163137)); // values[0] * phiinv[0]
585  CHECK(ratios2[1] == Approx(-0.504163137)); // symmetric move
586  }
587 }
588 
589 TEST_CASE("ReadMolecularOrbital GTO Ne", "[wavefunction]") { test_Ne(false); }
590 
591 TEST_CASE("ReadMolecularOrbital Numerical Ne", "[wavefunction]") { test_Ne(true); }
592 
593 TEST_CASE("ReadMolecularOrbital HCN", "[wavefunction]") {}
594 
595 void test_HCN(bool transform)
596 {
597  std::ostringstream section_name;
598  section_name << "HCN, transform orbitals to grid: " << (transform ? "T" : "F");
599 
600  SECTION(section_name.str())
601  {
602  Communicate* c = OHMMS::Controller;
603 
604  Libxml2Document doc;
605  bool okay = doc.parse("hcn.structure.xml");
606  REQUIRE(okay);
607  xmlNodePtr root = doc.getRoot();
608 
609  const SimulationCell simulation_cell;
610  auto ions_ptr = std::make_unique<ParticleSet>(simulation_cell);
611  auto& ions(*ions_ptr);
612  XMLParticleParser parse_ions(ions);
613  OhmmsXPathObject particleset_ion("//particleset[@name='ion0']", doc.getXPathContext());
614  REQUIRE(particleset_ion.size() == 1);
615  parse_ions.readXML(particleset_ion[0]);
616 
617  REQUIRE(ions.groups() == 3);
618  REQUIRE(ions.R.size() == 3);
619  ions.update();
620 
621  auto elec_ptr = std::make_unique<ParticleSet>(simulation_cell);
622  auto& elec(*elec_ptr);
623  XMLParticleParser parse_elec(elec);
624  OhmmsXPathObject particleset_elec("//particleset[@name='e']", doc.getXPathContext());
625  REQUIRE(particleset_elec.size() == 1);
626  parse_elec.readXML(particleset_elec[0]);
627 
628  REQUIRE(elec.groups() == 2);
629  REQUIRE(elec.R.size() == 14);
630 
631  elec.R = 0.0;
632 
633  elec.addTable(ions);
634  elec.update();
635 
636  Libxml2Document doc2;
637  okay = doc2.parse("hcn.wfnoj.xml");
638  REQUIRE(okay);
639  xmlNodePtr root2 = doc2.getRoot();
640 
641  WaveFunctionComponentBuilder::PSetMap particle_set_map;
642  particle_set_map.emplace(elec_ptr->getName(), std::move(elec_ptr));
643  particle_set_map.emplace(ions_ptr->getName(), std::move(ions_ptr));
644 
645  SPOSetBuilderFactory bf(c, elec, particle_set_map);
646 
647  OhmmsXPathObject MO_base("//determinantset", doc2.getXPathContext());
648  REQUIRE(MO_base.size() == 1);
649 
650  if (!transform)
651  {
652  // input file is set to transform GTO's to numerical orbitals by default
653  // use direct evaluation of GTO's
654  xmlSetProp(MO_base[0], castCharToXMLChar("transform"), castCharToXMLChar("no"));
655  xmlSetProp(MO_base[0], castCharToXMLChar("key"), castCharToXMLChar("GTO"));
656  }
657 
658  const auto bb_ptr = bf.createSPOSetBuilder(MO_base[0]);
659  auto& bb(*bb_ptr);
660 
661  OhmmsXPathObject slater_base("//determinant", doc2.getXPathContext());
662  auto sposet = bb.createSPOSet(slater_base[0]);
663 
664 
665  SPOSet::ValueVector values;
666  SPOSet::GradVector dpsi;
667  SPOSet::ValueVector d2psi;
668  values.resize(7);
669  dpsi.resize(7);
670  d2psi.resize(7);
671 
672  ParticleSet::SingleParticlePos newpos;
673  elec.makeMove(0, newpos);
674 
675  sposet->evaluateValue(elec, 0, values);
676 
677  CHECK(values[0] == Approx(0.009452265234));
678  CHECK(values[1] == Approx(0.02008357407));
679  CHECK(values[2] == Approx(0.4163749594));
680  CHECK(values[3] == Approx(-0.08854428343));
681  CHECK(values[4] == Approx(0.273158705));
682  CHECK(values[5] == Approx(0));
683  CHECK(values[6] == Approx(0));
684 
685  sposet->evaluateVGL(elec, 0, values, dpsi, d2psi);
686 
687 
688  // Generated form gen_mo.py for position [0.0, 0.0, 0.0]
689  CHECK(values[0] == Approx(0.009452265234));
690  CHECK(dpsi[0][0] == Approx(-0.05400764372));
691  CHECK(dpsi[0][1] == Approx(0));
692  CHECK(dpsi[0][2] == Approx(0));
693  CHECK(d2psi[0] == Approx(0.2532157143));
694 
695  CHECK(values[1] == Approx(0.02008357407));
696  CHECK(dpsi[1][0] == Approx(0.1009262252));
697  CHECK(dpsi[1][1] == Approx(0));
698  CHECK(dpsi[1][2] == Approx(0));
699  CHECK(d2psi[1] == Approx(0.3423520138));
700 
701  CHECK(values[2] == Approx(0.4163749594));
702  CHECK(dpsi[2][0] == Approx(-0.1202256419));
703  CHECK(dpsi[2][1] == Approx(0));
704  CHECK(dpsi[2][2] == Approx(0));
705  CHECK(d2psi[2] == Approx(-1.178149899));
706 
707  CHECK(values[3] == Approx(-0.08854428343));
708  CHECK(dpsi[3][0] == Approx(-0.004505552544));
709  CHECK(dpsi[3][1] == Approx(0));
710  CHECK(dpsi[3][2] == Approx(0));
711  CHECK(d2psi[3] == Approx(0.2838238091));
712 
713  CHECK(values[4] == Approx(0.273158705));
714  CHECK(dpsi[4][0] == Approx(-0.01125044248));
715  CHECK(dpsi[4][1] == Approx(0));
716  CHECK(dpsi[4][2] == Approx(0));
717  CHECK(d2psi[4] == Approx(-0.9173261582));
718 
719  CHECK(values[5] == Approx(0));
720  CHECK(dpsi[5][0] == Approx(0));
721  CHECK(dpsi[5][1] == Approx(0.4221165864));
722  CHECK(dpsi[5][2] == Approx(-0.08191634629));
723  CHECK(d2psi[5] == Approx(0));
724 
725  CHECK(values[6] == Approx(0));
726  CHECK(dpsi[6][0] == Approx(0));
727  CHECK(dpsi[6][1] == Approx(0.08191634629));
728  CHECK(dpsi[6][2] == Approx(0.4221165864));
729  CHECK(d2psi[6] == Approx(0));
730 
731  //==========Hessian and Grad Hessian Test==========
732  SPOSet::HessVector dhpsi;
733  dhpsi.resize(7);
734 
735  sposet->evaluateVGH(elec, 0, values, dpsi, dhpsi);
736 
737  // Generated from gen_mo.py for position [0.0, 0.0, 0.0]
738  CHECK(values[0] == Approx(0.009452265234));
739  CHECK(dpsi[0][0] == Approx(-0.05400764372));
740  CHECK(dpsi[0][1] == Approx(0));
741  CHECK(dpsi[0][2] == Approx(0));
742  //Hessian (xx,xy,xz,yy,yz,zz)
743  CHECK(dhpsi[0](0, 0) == Approx(0.3523924743));
744  CHECK(dhpsi[0](0, 1) == Approx(0));
745  CHECK(dhpsi[0](0, 2) == Approx(0));
746  CHECK(dhpsi[0](1, 1) == Approx(-0.04958838002));
747  CHECK(dhpsi[0](1, 2) == Approx(0));
748  CHECK(dhpsi[0](2, 2) == Approx(-0.04958838002));
749 
750  /////////////////////////////////////////////////////////////////////////////
751  //Now we're going to test the API's called by SPOSet for higher derivatives//
752  /////////////////////////////////////////////////////////////////////////////
753  SPOSet::ValueMatrix psiM(elec.R.size(), sposet->getOrbitalSetSize());
754  SPOSet::GradMatrix dpsiM(elec.R.size(), sposet->getOrbitalSetSize());
755  SPOSet::HessMatrix hesspsiV(elec.R.size(), sposet->getOrbitalSetSize());
756  SPOSet::GGGMatrix d3psiV(elec.R.size(), sposet->getOrbitalSetSize());
757 
758  sposet->evaluate_notranspose(elec, 0, elec.R.size(), psiM, dpsiM, hesspsiV, d3psiV);
759 
760 
761  // Generated from gen_mo.py for position [0.0, 0.0, 0.0]
762  CHECK(psiM[0][0] == Approx(0.009452265234));
763  CHECK(dpsiM[0][0][0] == Approx(-0.05400764372));
764  CHECK(dpsiM[0][0][1] == Approx(0));
765  CHECK(dpsiM[0][0][2] == Approx(0));
766  //Hessian (xx,xy,xz,yy,yz,zz)
767  CHECK(hesspsiV[0][0](0, 0) == Approx(0.3523924743));
768  CHECK(hesspsiV[0][0](0, 1) == Approx(0));
769  CHECK(hesspsiV[0][0](0, 2) == Approx(0));
770  CHECK(hesspsiV[0][0](1, 1) == Approx(-0.04958838002));
771  CHECK(hesspsiV[0][0](1, 2) == Approx(0));
772  CHECK(hesspsiV[0][0](2, 2) == Approx(-0.04958838002));
773 
774  //GradHessian (xxx,xxy,xxz,xyy,xyz,xzz,yyy,yyz,yzz,zzz)
775  CHECK(d3psiV[0][0][0](0, 0) == Approx(-2.241965465));
776  CHECK(d3psiV[0][0][0](0, 1) == Approx(0));
777  CHECK(d3psiV[0][0][0](0, 2) == Approx(0));
778  CHECK(d3psiV[0][0][0](1, 1) == Approx(0.3714481861));
779  CHECK(d3psiV[0][0][0](1, 2) == Approx(0));
780  CHECK(d3psiV[0][0][0](2, 2) == Approx(0.3714481861));
781  CHECK(d3psiV[0][0][1](1, 1) == Approx(0));
782  CHECK(d3psiV[0][0][1](1, 2) == Approx(0));
783  CHECK(d3psiV[0][0][1](2, 2) == Approx(0));
784  CHECK(d3psiV[0][0][2](2, 2) == Approx(0));
785 
786  // Generated from gen_mo.py for position [0.0, 0.0, 0.0]
787  CHECK(psiM[0][1] == Approx(0.02008357407));
788  CHECK(dpsiM[0][1][0] == Approx(0.1009262252));
789  CHECK(dpsiM[0][1][1] == Approx(0));
790  CHECK(dpsiM[0][1][2] == Approx(0));
791  //Hessian (xx,xy,xz,yy,yz,zz)
792  CHECK(hesspsiV[0][1](0, 0) == Approx(0.5298289497));
793  CHECK(hesspsiV[0][1](0, 1) == Approx(0));
794  CHECK(hesspsiV[0][1](0, 2) == Approx(0));
795  CHECK(hesspsiV[0][1](1, 1) == Approx(-0.09373846794));
796  CHECK(hesspsiV[0][1](1, 2) == Approx(0));
797  CHECK(hesspsiV[0][1](2, 2) == Approx(-0.09373846794));
798 
799  //GradHessian (xxx,xxy,xxz,xyy,xyz,xzz,yyy,yyz,yzz,zzz)
800  CHECK(d3psiV[0][1][0](0, 0) == Approx(2.594787656));
801  CHECK(d3psiV[0][1][0](0, 1) == Approx(0));
802  CHECK(d3psiV[0][1][0](0, 2) == Approx(0));
803  CHECK(d3psiV[0][1][0](1, 1) == Approx(-0.5720485625));
804  CHECK(d3psiV[0][1][0](1, 2) == Approx(0));
805  CHECK(d3psiV[0][1][0](2, 2) == Approx(-0.5720485625));
806  CHECK(d3psiV[0][1][1](1, 1) == Approx(0));
807  CHECK(d3psiV[0][1][1](1, 2) == Approx(0));
808  CHECK(d3psiV[0][1][1](2, 2) == Approx(0));
809  CHECK(d3psiV[0][1][2](2, 2) == Approx(0));
810  // Generated from gen_mo.py for position [0.0, 0.0, 0.0]
811  CHECK(psiM[0][2] == Approx(0.4163749594));
812  CHECK(dpsiM[0][2][0] == Approx(-0.1202256419));
813  CHECK(dpsiM[0][2][1] == Approx(0));
814  CHECK(dpsiM[0][2][2] == Approx(0));
815  //Hessian (xx,xy,xz,yy,yz,zz)
816  CHECK(hesspsiV[0][2](0, 0) == Approx(-0.02607695984));
817  CHECK(hesspsiV[0][2](0, 1) == Approx(0));
818  CHECK(hesspsiV[0][2](0, 2) == Approx(0));
819  CHECK(hesspsiV[0][2](1, 1) == Approx(-0.5760364698));
820  CHECK(hesspsiV[0][2](1, 2) == Approx(0));
821  CHECK(hesspsiV[0][2](2, 2) == Approx(-0.5760364698));
822  //GradHessian (xxx,xxy,xxz,xyy,xyz,xzz,yyy,yyz,yzz,zzz)
823  CHECK(d3psiV[0][2][0](0, 0) == Approx(-0.227147312));
824  CHECK(d3psiV[0][2][0](0, 1) == Approx(0));
825  CHECK(d3psiV[0][2][0](0, 2) == Approx(0));
826  CHECK(d3psiV[0][2][0](1, 1) == Approx(0.2992015499));
827  CHECK(d3psiV[0][2][0](1, 2) == Approx(0));
828  CHECK(d3psiV[0][2][0](2, 2) == Approx(0.2992015499));
829  CHECK(d3psiV[0][2][1](1, 1) == Approx(0));
830  CHECK(d3psiV[0][2][1](1, 2) == Approx(0));
831  CHECK(d3psiV[0][2][1](2, 2) == Approx(0));
832  CHECK(d3psiV[0][2][2](2, 2) == Approx(0));
833 
834 
835  //Move electron 0 to some nontrivial position:
836  ParticleSet::SingleParticlePos disp(0.02, -0.1, 0.05);
837  elec.makeMove(0, disp);
838 
839 
840  SPOSet::GradMatrix dionpsi(elec.R.size(), sposet->getOrbitalSetSize());
841  SPOSet::HessMatrix diongradpsi(elec.R.size(), sposet->getOrbitalSetSize());
842  SPOSet::GradMatrix dionlaplpsi(elec.R.size(), sposet->getOrbitalSetSize());
843 
844  sposet->evaluateGradSource(elec, 0, elec.R.size(), ions, 0, dionpsi, diongradpsi, dionlaplpsi);
845 
846  //============== Ion 0 Component 0 ===================
847  CHECK(dionpsi[0][0][0] == Approx(0.0453112082));
848  CHECK(diongradpsi[0][0](0, 0) == Approx(-0.2943513994));
849  CHECK(diongradpsi[0][0](0, 1) == Approx(0.030468047));
850  CHECK(diongradpsi[0][0](0, 2) == Approx(-0.0152340235));
851  CHECK(dionlaplpsi[0][0][0] == Approx(1.333755581));
852  CHECK(dionpsi[0][1][0] == Approx(-0.0006473819623));
853  CHECK(diongradpsi[0][1](0, 0) == Approx(0.0004713407512));
854  CHECK(diongradpsi[0][1](0, 1) == Approx(-0.0001254975603));
855  CHECK(diongradpsi[0][1](0, 2) == Approx(6.274878013e-05));
856  CHECK(dionlaplpsi[0][1][0] == Approx(0.001057940846));
857  CHECK(dionpsi[0][2][0] == Approx(0.265578336));
858  CHECK(diongradpsi[0][2](0, 0) == Approx(-0.08685804115));
859  CHECK(diongradpsi[0][2](0, 1) == Approx(0.05438178417));
860  CHECK(diongradpsi[0][2](0, 2) == Approx(-0.02719089209));
861  CHECK(dionlaplpsi[0][2][0] == Approx(-1.882489819));
862  CHECK(dionpsi[0][3][0] == Approx(-0.06444305979));
863  CHECK(diongradpsi[0][3](0, 0) == Approx(-0.002013151923));
864  CHECK(diongradpsi[0][3](0, 1) == Approx(-0.002535923431));
865  CHECK(diongradpsi[0][3](0, 2) == Approx(0.001267961716));
866  CHECK(dionlaplpsi[0][3][0] == Approx(0.4547401581));
867  CHECK(dionpsi[0][4][0] == Approx(0.1454357726));
868  CHECK(diongradpsi[0][4](0, 0) == Approx(-0.2330499431));
869  CHECK(diongradpsi[0][4](0, 1) == Approx(0.09667641762));
870  CHECK(diongradpsi[0][4](0, 2) == Approx(-0.04833820881));
871  CHECK(dionlaplpsi[0][4][0] == Approx(-0.9197558839));
872  CHECK(dionpsi[0][5][0] == Approx(-0.04329985085));
873  CHECK(diongradpsi[0][5](0, 0) == Approx(0.09051993304));
874  CHECK(diongradpsi[0][5](0, 1) == Approx(0.382375474));
875  CHECK(diongradpsi[0][5](0, 2) == Approx(-0.07043361927));
876  CHECK(dionlaplpsi[0][5][0] == Approx(0.2201672051));
877  CHECK(dionpsi[0][6][0] == Approx(0.01207541177));
878  CHECK(diongradpsi[0][6](0, 0) == Approx(-0.02524405435));
879  CHECK(diongradpsi[0][6](0, 1) == Approx(0.0800332842));
880  CHECK(diongradpsi[0][6](0, 2) == Approx(0.3929818664));
881  CHECK(dionlaplpsi[0][6][0] == Approx(-0.0614000824));
882 
883 
884  sposet->evaluateGradSource(elec, 0, elec.R.size(), ions, 1, dionpsi, diongradpsi, dionlaplpsi);
885 
886  //============== Ion 1 Component 1 ===================
887  CHECK(dionpsi[0][0][1] == Approx(0.0001412373768));
888  CHECK(diongradpsi[0][0](1, 0) == Approx(0.0001124265646));
889  CHECK(diongradpsi[0][0](1, 1) == Approx(-0.001383378615));
890  CHECK(diongradpsi[0][0](1, 2) == Approx(-1.449757545e-05));
891  CHECK(dionlaplpsi[0][0][1] == Approx(-0.001252043663));
892  CHECK(dionpsi[0][1][1] == Approx(-0.01029290716));
893  CHECK(diongradpsi[0][1](1, 0) == Approx(-0.06178485148));
894  CHECK(diongradpsi[0][1](1, 1) == Approx(0.0971577216));
895  CHECK(diongradpsi[0][1](1, 2) == Approx(0.002885675005));
896  CHECK(dionlaplpsi[0][1][1] == Approx(-0.1403103458));
897  CHECK(dionpsi[0][2][1] == Approx(-0.0230872583));
898  CHECK(diongradpsi[0][2](1, 0) == Approx(-0.02537847709));
899  CHECK(diongradpsi[0][2](1, 1) == Approx(0.2268946564));
900  CHECK(diongradpsi[0][2](1, 2) == Approx(0.001988963201));
901  CHECK(dionlaplpsi[0][2][1] == Approx(0.2028851421));
902  CHECK(dionpsi[0][3][1] == Approx(0.01850231814));
903  CHECK(diongradpsi[0][3](1, 0) == Approx(0.05709948475));
904  CHECK(diongradpsi[0][3](1, 1) == Approx(-0.1776515965));
905  CHECK(diongradpsi[0][3](1, 2) == Approx(-0.003685792479));
906  CHECK(dionlaplpsi[0][3][1] == Approx(-0.1280699725));
907  CHECK(dionpsi[0][4][1] == Approx(-0.02136209962));
908  CHECK(diongradpsi[0][4](1, 0) == Approx(-0.03836586276));
909  CHECK(diongradpsi[0][4](1, 1) == Approx(0.2084578148));
910  CHECK(diongradpsi[0][4](1, 2) == Approx(0.002581590766));
911  CHECK(dionlaplpsi[0][4][1] == Approx(0.1792683544));
912  CHECK(dionpsi[0][5][1] == Approx(-0.1942343714));
913  CHECK(diongradpsi[0][5](1, 0) == Approx(-0.3037357197));
914  CHECK(diongradpsi[0][5](1, 1) == Approx(-0.09561978734));
915  CHECK(diongradpsi[0][5](1, 2) == Approx(0.02118492506));
916  CHECK(dionlaplpsi[0][5][1] == Approx(0.6410434658));
917  CHECK(dionpsi[0][6][1] == Approx(-0.03930992259));
918  CHECK(diongradpsi[0][6](1, 0) == Approx(-0.06331544695));
919  CHECK(diongradpsi[0][6](1, 1) == Approx(-0.002807368817));
920  CHECK(diongradpsi[0][6](1, 2) == Approx(-0.02801340823));
921  CHECK(dionlaplpsi[0][6][1] == Approx(0.1369061053));
922 
923  sposet->evaluateGradSource(elec, 0, elec.R.size(), ions, 2, dionpsi, diongradpsi, dionlaplpsi);
924 
925  //============== Ion 2 Component 2 ===================
926  CHECK(dionpsi[0][0][2] == Approx(1.302648961e-06));
927  CHECK(diongradpsi[0][0](2, 0) == Approx(1.865129579e-06));
928  CHECK(diongradpsi[0][0](2, 1) == Approx(6.142092043e-08));
929  CHECK(diongradpsi[0][0](2, 2) == Approx(2.602225618e-05));
930  CHECK(dionlaplpsi[0][0][2] == Approx(1.234692903e-06));
931  CHECK(dionpsi[0][1][2] == Approx(3.248738084e-07));
932  CHECK(diongradpsi[0][1](2, 0) == Approx(-2.044420189e-06));
933  CHECK(diongradpsi[0][1](2, 1) == Approx(-7.011145137e-08));
934  CHECK(diongradpsi[0][1](2, 2) == Approx(6.532522353e-06));
935  CHECK(dionlaplpsi[0][1][2] == Approx(-6.10958506e-06));
936  CHECK(dionpsi[0][2][2] == Approx(3.264249981e-06));
937  CHECK(diongradpsi[0][2](2, 0) == Approx(2.820971234e-05));
938  CHECK(diongradpsi[0][2](2, 1) == Approx(9.405184964e-07));
939  CHECK(diongradpsi[0][2](2, 2) == Approx(6.481420782e-05));
940  CHECK(dionlaplpsi[0][2][2] == Approx(5.73961989e-05));
941  CHECK(dionpsi[0][3][2] == Approx(0.0001288974413));
942  CHECK(diongradpsi[0][3](2, 0) == Approx(0.0002840756879));
943  CHECK(diongradpsi[0][3](2, 1) == Approx(9.281700408e-06));
944  CHECK(diongradpsi[0][3](2, 2) == Approx(0.002573308008));
945  CHECK(dionlaplpsi[0][3][2] == Approx(0.0003025314443));
946  CHECK(dionpsi[0][4][2] == Approx(-7.300043903e-05));
947  CHECK(diongradpsi[0][4](2, 0) == Approx(-0.0001000016834));
948  CHECK(diongradpsi[0][4](2, 1) == Approx(-3.233243534e-06));
949  CHECK(diongradpsi[0][4](2, 2) == Approx(-0.001458391774));
950  CHECK(dionlaplpsi[0][4][2] == Approx(-3.546690719e-05));
951  CHECK(dionpsi[0][5][2] == Approx(2.910525987e-06));
952  CHECK(diongradpsi[0][5](2, 0) == Approx(1.307065133e-05));
953  CHECK(diongradpsi[0][5](2, 1) == Approx(1.560390706e-06));
954  CHECK(diongradpsi[0][5](2, 2) == Approx(-2.92731811e-06));
955  CHECK(dionlaplpsi[0][5][2] == Approx(3.797816228e-05));
956  CHECK(dionpsi[0][6][2] == Approx(-1.56074936e-05));
957  CHECK(diongradpsi[0][6](2, 0) == Approx(-7.009049656e-05));
958  CHECK(diongradpsi[0][6](2, 1) == Approx(-2.048666792e-06));
959  CHECK(diongradpsi[0][6](2, 2) == Approx(2.967709412e-06));
960  CHECK(dionlaplpsi[0][6][2] == Approx(-0.0002018111858));
961 
962  //Same tests as before, but for the gradient only.
963 
964  sposet->evaluateGradSource(elec, 0, elec.R.size(), ions, 0, dionpsi);
965  //============== Ion 0 Component 0 ===================
966  CHECK(dionpsi[0][0][0] == Approx(0.0453112082));
967  CHECK(dionpsi[0][1][0] == Approx(-0.0006473819623));
968  CHECK(dionpsi[0][2][0] == Approx(0.265578336));
969  CHECK(dionpsi[0][3][0] == Approx(-0.06444305979));
970  CHECK(dionpsi[0][4][0] == Approx(0.1454357726));
971  CHECK(dionpsi[0][5][0] == Approx(-0.04329985085));
972  CHECK(dionpsi[0][6][0] == Approx(0.01207541177));
973 
974  sposet->evaluateGradSource(elec, 0, elec.R.size(), ions, 1, dionpsi);
975  //============== Ion 1 Component 1 ===================
976  CHECK(dionpsi[0][0][1] == Approx(0.0001412373768));
977  CHECK(dionpsi[0][1][1] == Approx(-0.01029290716));
978  CHECK(dionpsi[0][2][1] == Approx(-0.0230872583));
979  CHECK(dionpsi[0][3][1] == Approx(0.01850231814));
980  CHECK(dionpsi[0][4][1] == Approx(-0.02136209962));
981  CHECK(dionpsi[0][5][1] == Approx(-0.1942343714));
982  CHECK(dionpsi[0][6][1] == Approx(-0.03930992259));
983 
984  sposet->evaluateGradSource(elec, 0, elec.R.size(), ions, 2, dionpsi);
985  //============== Ion 2 Component 2 ===================
986  CHECK(dionpsi[0][0][2] == Approx(1.302648961e-06));
987  CHECK(dionpsi[0][1][2] == Approx(3.248738084e-07));
988  CHECK(dionpsi[0][2][2] == Approx(3.264249981e-06));
989  CHECK(dionpsi[0][3][2] == Approx(0.0001288974413));
990  CHECK(dionpsi[0][4][2] == Approx(-7.300043903e-05));
991  CHECK(dionpsi[0][5][2] == Approx(2.910525987e-06));
992  CHECK(dionpsi[0][6][2] == Approx(-1.56074936e-05));
993 
994 
995  //Finally, going to test evaluateGradSourceRow. Same template and reference
996  //values as above.
997  SPOSet::GradVector dionpsivec;
998  dionpsivec.resize(7);
999 
1000  sposet->evaluateGradSourceRow(elec, 0, ions, 0, dionpsivec);
1001  //============== Ion 0 Component 0 ===================
1002  CHECK(dionpsivec[0][0] == Approx(0.0453112082));
1003  CHECK(dionpsivec[1][0] == Approx(-0.0006473819623));
1004  CHECK(dionpsivec[2][0] == Approx(0.265578336));
1005  CHECK(dionpsivec[3][0] == Approx(-0.06444305979));
1006  CHECK(dionpsivec[4][0] == Approx(0.1454357726));
1007  CHECK(dionpsivec[5][0] == Approx(-0.04329985085));
1008  CHECK(dionpsivec[6][0] == Approx(0.01207541177));
1009 
1010  sposet->evaluateGradSourceRow(elec, 0, ions, 1, dionpsivec);
1011  //============== Ion 1 Component 1 ===================
1012  CHECK(dionpsivec[0][1] == Approx(0.0001412373768));
1013  CHECK(dionpsivec[1][1] == Approx(-0.01029290716));
1014  CHECK(dionpsivec[2][1] == Approx(-0.0230872583));
1015  CHECK(dionpsivec[3][1] == Approx(0.01850231814));
1016  CHECK(dionpsivec[4][1] == Approx(-0.02136209962));
1017  CHECK(dionpsivec[5][1] == Approx(-0.1942343714));
1018  CHECK(dionpsivec[6][1] == Approx(-0.03930992259));
1019 
1020  sposet->evaluateGradSourceRow(elec, 0, ions, 2, dionpsivec);
1021  //============== Ion 2 Component 2 ===================
1022  CHECK(dionpsivec[0][2] == Approx(1.302648961e-06));
1023  CHECK(dionpsivec[1][2] == Approx(3.248738084e-07));
1024  CHECK(dionpsivec[2][2] == Approx(3.264249981e-06));
1025  CHECK(dionpsivec[3][2] == Approx(0.0001288974413));
1026  CHECK(dionpsivec[4][2] == Approx(-7.300043903e-05));
1027  CHECK(dionpsivec[5][2] == Approx(2.910525987e-06));
1028  CHECK(dionpsivec[6][2] == Approx(-1.56074936e-05));
1029  }
1030 }
1031 
1032 TEST_CASE("ReadMolecularOrbital GTO HCN", "[wavefunction]") { test_HCN(false); }
1033 
1034 TEST_CASE("ReadMolecularOrbital Numerical HCN", "[wavefunction]") { test_HCN(true); }
1035 
1036 } // namespace qmcplusplus
qmcplusplus::SPOSet::HessVector
OrbitalSetTraits< ValueType >::HessVector HessVector
Definition: SPOSet.h:53
qmcplusplus::TinyVector
Fixed-size array.
Definition: OhmmsTinyMeta.h:30
qmcplusplus::test_HCN
void test_HCN(bool transform)
Definition: test_MO.cpp:595
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
Libxml2Document::getXPathContext
xmlXPathContextPtr getXPathContext()
Definition: Libxml2Doc.cpp:100
qmcplusplus::TEST_CASE
TEST_CASE("complex_helper", "[type_traits]")
Definition: test_complex_helper.cpp:38
qmcplusplus::VirtualParticleSet
A ParticleSet that handles virtual moves of a selected particle of a given physical ParticleSet Virtu...
Definition: VirtualParticleSet.h:39
qmcplusplus::VirtualParticleSet::makeMoves
void makeMoves(const ParticleSet &refp, int jel, const std::vector< PosType > &deltaV, bool sphere=false, int iat=-1)
move virtual particles to new postions and update distance tables
Definition: VirtualParticleSet.cpp:129
Libxml2Document::getRoot
xmlNodePtr getRoot()
Definition: Libxml2Doc.h:88
qmcplusplus::SPOSet::ValueMatrix
OrbitalSetTraits< ValueType >::ValueMatrix ValueMatrix
Definition: SPOSet.h:50
qmcplusplus::test_EtOH_mw
void test_EtOH_mw(bool transform)
Definition: test_MO.cpp:285
qmcplusplus::ParticleSet::update
void update(bool skipSK=false)
update the internal data
Definition: ParticleSet.cpp:349
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::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::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::ParticleSet
Specialized paritlce class for atomistic simulations.
Definition: ParticleSet.h:55
qmcplusplus::REQUIRE
REQUIRE(std::filesystem::exists(filename))
qmcplusplus::SPOSet::ValueVector
OrbitalSetTraits< ValueType >::ValueVector ValueVector
Definition: SPOSet.h:49
qmcplusplus::XMLParticleParser::readXML
bool readXML(xmlNodePtr cur)
process xmlnode <particleset/> which contains everything about the particle set to initialize ...
Definition: XMLParticleIO.cpp:90
qmcplusplus::ParticleSet::R
ParticlePos R
Position.
Definition: ParticleSet.h:79
qmcplusplus::test_He_mw
void test_He_mw(bool transform)
Definition: test_MO.cpp:136
qmcplusplus::test_Ne
void test_Ne(bool transform)
Definition: test_MO.cpp:465
qmcplusplus::RefVector
std::vector< std::reference_wrapper< T > > RefVector
Definition: template_types.hpp:32
qmcplusplus::SPOSet::GradVector
OrbitalSetTraits< ValueType >::GradVector GradVector
Definition: SPOSet.h:51
qmcplusplus::ParticleSet::makeMove
void makeMove(Index_t iat, const SingleParticlePos &displ, bool maybe_accept=true)
move the iat-th particle to active_pos_
Definition: ParticleSet.cpp:382
qmcplusplus::test_He
void test_He(bool transform)
Definition: test_MO.cpp:26
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::ResourceCollectionTeamLock
handles acquire/release resource by the consumer (RefVectorWithLeader type).
Definition: ResourceCollection.h:58
qmcplusplus::SPOSetBuilderFactory::createSPOSetBuilder
std::unique_ptr< SPOSetBuilder > createSPOSetBuilder(xmlNodePtr rootNode)
Definition: SPOSetBuilderFactory.cpp:71
qmcplusplus::SpeciesSet
Custom container for set of attributes for a set of species.
Definition: SpeciesSet.h:33
qmcplusplus::SPOSet::GGGMatrix
OrbitalSetTraits< ValueType >::GradHessMatrix GGGMatrix
Definition: SPOSet.h:56
castCharToXMLChar
xmlChar * castCharToXMLChar(char *c)
unsafe char* to xmlChar* cast
Definition: libxmldefs.h:70
qmcplusplus::SPOSet::HessMatrix
OrbitalSetTraits< ValueType >::HessMatrix HessMatrix
Definition: SPOSet.h:54