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) 2017 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Mark Dewing, markdewing@gmail.com, University of Illinois at Urbana-Champaign
8 // Ye Luo, yeluo@anl.gov, Argonne National Laboratory
9 //
10 // File created by: Mark Dewing, markdewing@gmail.com, University of Illinois at Urbana-Champaign
11 //////////////////////////////////////////////////////////////////////////////////////
12
13
14 #include "catch.hpp"
15
16 #include "Configuration.h"
17 #include "Numerics/Quadrature.h"
18 #include "QMCHamiltonians/ECPComponentBuilder.h"
19 #include "QMCHamiltonians/NonLocalECPComponent.h"
20 #include "QMCHamiltonians/SOECPComponent.h"
21
22 //for wavefunction
23 #include "OhmmsData/Libxml2Doc.h"
24 #include "QMCWaveFunctions/WaveFunctionComponent.h"
25 #include "QMCWaveFunctions/TrialWaveFunction.h"
26 #include "QMCWaveFunctions/Jastrow/BsplineFunctor.h"
27 #include "QMCWaveFunctions/Jastrow/RadialJastrowBuilder.h"
28 #include "QMCWaveFunctions/Fermion/DiracDeterminant.h"
29 #include "QMCWaveFunctions/SpinorSet.h"
30 //for nonlocal moves
31 #include "QMCHamiltonians/NonLocalTOperator.h"
32
33
34 //for Hamiltonian manipulations.
35 #include "Particle/ParticleSet.h"
36 #include "Particle/ParticleSetPool.h"
37 #include "LongRange/EwaldHandler3D.h"
38
39 #ifdef QMC_COMPLEX //This is for the spinor test.
40 #include "QMCWaveFunctions/ElectronGas/ElectronGasComplexOrbitalBuilder.h"
41 #endif
42
43 namespace qmcplusplus
44 {
getSplinedSOPot(SOECPComponent * so_comp,int l,double r)45 QMCTraits::RealType getSplinedSOPot(SOECPComponent* so_comp, int l, double r) { return so_comp->sopp_m[l]->splint(r); }
46
47 TEST_CASE("CheckSphericalIntegration", "[hamiltonian]")
48 {
49 // Use the built-in quadrature rule check
50 for (int quadrature_index = 1; quadrature_index < 8; quadrature_index++)
51 {
52 Quadrature3D<QMCTraits::RealType> myRule(quadrature_index, false);
53 REQUIRE(myRule.quad_ok);
54 }
55 }
56
57 TEST_CASE("ReadFileBuffer_no_file", "[hamiltonian]")
58 {
59 ReadFileBuffer buf(NULL);
60 bool open_okay = buf.open_file("does_not_exist");
61 REQUIRE(open_okay == false);
62 }
63
64 TEST_CASE("ReadFileBuffer_simple_serial", "[hamiltonian]")
65 {
66 // Initializing with no Communicate pointer under MPI,
67 // this will read the file on every node. Should be okay
68 // for testing purposes.
69 ReadFileBuffer buf(NULL);
70 bool open_okay = buf.open_file("simple.txt");
71 REQUIRE(open_okay == true);
72
73 bool read_okay = buf.read_contents();
74 REQUIRE(read_okay);
75 REQUIRE(buf.length == 14);
76 REQUIRE(std::string("File contents\n") == buf.contents());
77 }
78
79 TEST_CASE("ReadFileBuffer_simple_mpi", "[hamiltonian]")
80 {
81 Communicate* c = OHMMS::Controller;
82
83 ReadFileBuffer buf(c);
84 bool open_okay = buf.open_file("simple.txt");
85 REQUIRE(open_okay == true);
86
87 bool read_okay = buf.read_contents();
88 REQUIRE(read_okay);
89 REQUIRE(buf.length == 14);
90 REQUIRE(std::string("File contents\n") == buf.contents());
91 }
92
93 TEST_CASE("ReadFileBuffer_ecp", "[hamiltonian]")
94 {
95 Communicate* c = OHMMS::Controller;
96
97 ECPComponentBuilder ecp("test_read_ecp", c);
98
99 bool okay = ecp.read_pp_file("C.BFD.xml");
100 REQUIRE(okay);
101
102 REQUIRE(ecp.Zeff == 4);
103
104 // TODO: add more checks that pseudopotential file was read correctly
105 }
106
107 TEST_CASE("ReadFileBuffer_sorep", "[hamiltonian]")
108 {
109 Communicate* c = OHMMS::Controller;
110
111 ECPComponentBuilder ecp("test_read_sorep", c);
112
113 bool okay = ecp.read_pp_file("so_ecp_test.xml");
114 REQUIRE(okay);
115
116 REQUIRE(ecp.Zeff == 13);
117
118 double rlist[5] = {0.001, 0.500, 1.000, 2.000, 10.000};
119 double so_p[5] = {0.999999000005, 0.778800783071, 0.3678794411714, 0.01831563888873418, 0.000};
120 double so_d[5] = {9.99998000e-01, 6.06530660e-01, 1.35335283e-01, 3.35462628e-04, 0.000};
121 double so_f[5] = {9.99997000e-01, 4.72366553e-01, 4.97870684e-02, 6.14421235e-06, 0.000};
122
123 for (int i = 0; i < 5; i++)
124 {
125 double r = rlist[i];
126 double so_p_ref = so_p[i];
127 double so_d_ref = so_d[i];
128 double so_f_ref = so_f[i];
129
130 double so_p_val = getSplinedSOPot(ecp.pp_so.get(), 0, r);
131 double so_d_val = getSplinedSOPot(ecp.pp_so.get(), 1, r);
132 double so_f_val = getSplinedSOPot(ecp.pp_so.get(), 2, r);
133
134 REQUIRE(so_p_val == Approx(so_p_ref));
135 REQUIRE(so_d_val == Approx(so_d_ref));
136 REQUIRE(so_f_val == Approx(so_f_ref));
137 }
138
139 // TODO: add more checks that pseudopotential file was read correctly
140 }
141
142
143 TEST_CASE("ReadFileBuffer_reopen", "[hamiltonian]")
144 {
145 // Initializing with no Communicate pointer under MPI,
146 // this will read the file on every node. Should be okay
147 // for testing purposes.
148 ReadFileBuffer buf(NULL);
149 bool open_okay = buf.open_file("simple.txt");
150 REQUIRE(open_okay == true);
151
152 bool read_okay = buf.read_contents();
153 REQUIRE(read_okay);
154 REQUIRE(buf.length == 14);
155
156 open_okay = buf.open_file("C.BFD.xml");
157 REQUIRE(open_okay == true);
158
159 read_okay = buf.read_contents();
160 REQUIRE(read_okay);
161 REQUIRE(buf.length > 14);
162 }
163
copyGridUnrotatedForTest(NonLocalECPComponent & nlpp)164 void copyGridUnrotatedForTest(NonLocalECPComponent& nlpp) { nlpp.rrotsgrid_m = nlpp.sgridxyz_m; }
copyGridUnrotatedForTest(SOECPComponent & sopp)165 void copyGridUnrotatedForTest(SOECPComponent& sopp) { sopp.rrotsgrid_m = sopp.sgridxyz_m; }
166
167 TEST_CASE("Evaluate_ecp", "[hamiltonian]")
168 {
169 typedef QMCTraits::RealType RealType;
170 typedef QMCTraits::ValueType ValueType;
171 typedef QMCTraits::PosType PosType;
172
173 Communicate* c = OHMMS::Controller;
174
175 //Cell definition:
176
177 CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> Lattice;
178 Lattice.BoxBConds = true; // periodic
179 Lattice.R.diagonal(20);
180 Lattice.LR_dim_cutoff = 15;
181 Lattice.reset();
182
183
184 ParticleSet ions;
185 ParticleSet elec;
186
187 ions.setName("ion0");
188 ions.create(2);
189 ions.R[0][0] = 0.0;
190 ions.R[0][1] = 0.0;
191 ions.R[0][2] = 0.0;
192 ions.R[1][0] = 6.0;
193 ions.R[1][1] = 0.0;
194 ions.R[1][2] = 0.0;
195
196
197 SpeciesSet& ion_species = ions.getSpeciesSet();
198 int pIdx = ion_species.addSpecies("Na");
199 int pChargeIdx = ion_species.addAttribute("charge");
200 int iatnumber = ion_species.addAttribute("atomic_number");
201 ion_species(pChargeIdx, pIdx) = 1;
202 ion_species(iatnumber, pIdx) = 11;
203 ions.Lattice = Lattice;
204 ions.createSK();
205
206
207 elec.Lattice = Lattice;
208 elec.setName("e");
209 std::vector<int> agroup(2, 1);
210 elec.create(agroup);
211 elec.R[0][0] = 2.0;
212 elec.R[0][1] = 0.0;
213 elec.R[0][2] = 0.0;
214 elec.R[1][0] = 3.0;
215 elec.R[1][1] = 0.0;
216 elec.R[1][2] = 0.0;
217
218
219 SpeciesSet& tspecies = elec.getSpeciesSet();
220 int upIdx = tspecies.addSpecies("u");
221 int downIdx = tspecies.addSpecies("d");
222 int chargeIdx = tspecies.addAttribute("charge");
223 int massIdx = tspecies.addAttribute("mass");
224 tspecies(chargeIdx, upIdx) = -1;
225 tspecies(chargeIdx, downIdx) = -1;
226 tspecies(massIdx, upIdx) = 1.0;
227 tspecies(massIdx, downIdx) = 1.0;
228
229 elec.createSK();
230
231 ParticleSetPool ptcl = ParticleSetPool(c);
232
233 ions.resetGroups();
234
235 // The call to resetGroups is needed transfer the SpeciesSet
236 // settings to the ParticleSet
237 elec.resetGroups();
238
239 //Cool. Now to construct a wavefunction with 1 and 2 body jastrow (no determinant)
240 TrialWaveFunction psi;
241
242 //Add the two body jastrow
243 const char* particles = "<tmp> \
244 <jastrow name=\"J2\" type=\"Two-Body\" function=\"Bspline\" print=\"yes\"> \
245 <correlation speciesA=\"u\" speciesB=\"d\" rcut=\"10\" size=\"8\"> \
246 <coefficients id=\"ud\" type=\"Array\"> 2.015599059 1.548994099 1.17959447 0.8769687661 0.6245736507 0.4133517767 0.2333851935 0.1035636904</coefficients> \
247 </correlation> \
248 </jastrow> \
249 </tmp> \
250 ";
251 Libxml2Document doc;
252 bool okay = doc.parseFromString(particles);
253 REQUIRE(okay);
254
255 xmlNodePtr root = doc.getRoot();
256
257 xmlNodePtr jas2 = xmlFirstElementChild(root);
258
259 RadialJastrowBuilder jastrow(c, elec);
260 psi.addComponent(jastrow.buildComponent(jas2));
261 // Done with two body jastrow.
262
263 //Add the one body jastrow.
264 const char* particles2 = "<tmp> \
265 <jastrow name=\"J1\" type=\"One-Body\" function=\"Bspline\" source=\"ion0\" print=\"yes\"> \
266 <correlation elementType=\"Na\" rcut=\"10\" size=\"10\" cusp=\"0\"> \
267 <coefficients id=\"eNa\" type=\"Array\"> 1.244201343 -1.188935609 -1.840397253 -1.803849126 -1.612058635 -1.35993202 -1.083353212 -0.8066295188 -0.5319252448 -0.3158819772</coefficients> \
268 </correlation> \
269 </jastrow> \
270 </tmp> \
271 ";
272 bool okay3 = doc.parseFromString(particles2);
273 REQUIRE(okay3);
274
275 root = doc.getRoot();
276
277 xmlNodePtr jas1 = xmlFirstElementChild(root);
278
279 RadialJastrowBuilder jastrow1bdy(c, elec, ions);
280 psi.addComponent(jastrow1bdy.buildComponent(jas1));
281
282 //Now we set up the nonlocal ECP component.
283 ECPComponentBuilder ecp("test_read_ecp", c);
284
285 bool okay2 = ecp.read_pp_file("Na.BFD.xml");
286
287 NonLocalECPComponent* nlpp = ecp.pp_nonloc.get();
288
289 //This line is required because the randomized quadrature Lattice is set by
290 //random number generator in NonLocalECPotential. We take the unrotated
291 //quadrature Lattice instead...
292 copyGridUnrotatedForTest(*nlpp);
293
294 const int myTableIndex = elec.addTable(ions);
295
296 const auto& myTable = elec.getDistTable(myTableIndex);
297
298 // update all distance tables
299 ions.update();
300 elec.update();
301
302 //Need to set up temporary data for this configuration in trial wavefunction. Needed for ratios.
303 double logpsi = psi.evaluateLog(elec);
304 REQUIRE(logpsi == Approx(5.1497823982));
305
306 double Value1(0.0);
307 //Using SoA distance tables, hence the guard.
308 for (int jel = 0; jel < elec.getTotalNum(); jel++)
309 {
310 const auto& dist = myTable.getDistRow(jel);
311 const auto& displ = myTable.getDisplRow(jel);
312 for (int iat = 0; iat < ions.getTotalNum(); iat++)
313 if (nlpp != nullptr && dist[iat] < nlpp->getRmax())
314 Value1 += nlpp->evaluateOne(elec, iat, psi, jel, dist[iat], -displ[iat], false);
315 }
316 //These numbers are validated against an alternate code path via wavefunction tester.
317 REQUIRE(Value1 == Approx(6.9015710211e-02));
318
319 opt_variables_type optvars;
320 std::vector<ValueType> dlogpsi;
321 std::vector<ValueType> dhpsioverpsi;
322
323 psi.checkInVariables(optvars);
324 optvars.resetIndex();
325 const int NumOptimizables(optvars.size());
326 psi.checkOutVariables(optvars);
327 dlogpsi.resize(NumOptimizables, ValueType(0));
328 dhpsioverpsi.resize(NumOptimizables, ValueType(0));
329 psi.evaluateDerivatives(elec, optvars, dlogpsi, dhpsioverpsi);
330 REQUIRE(std::real(dlogpsi[0]) == Approx(-0.2211666667));
331 REQUIRE(std::real(dlogpsi[2]) == Approx(-0.1215));
332 REQUIRE(std::real(dlogpsi[3]) == Approx(0.0));
333 REQUIRE(std::real(dlogpsi[9]) == Approx(-0.0853333333));
334 REQUIRE(std::real(dlogpsi[10]) == Approx(-0.745));
335
336 REQUIRE(std::real(dhpsioverpsi[0]) == Approx(-0.6463306581));
337 REQUIRE(std::real(dhpsioverpsi[2]) == Approx(1.5689981479));
338 REQUIRE(std::real(dhpsioverpsi[3]) == Approx(0.0));
339 REQUIRE(std::real(dhpsioverpsi[9]) == Approx(0.279561213));
340 REQUIRE(std::real(dhpsioverpsi[10]) == Approx(-0.3968828778));
341
342 Value1 = 0.0;
343 //Using SoA distance tables, hence the guard.
344 for (int jel = 0; jel < elec.getTotalNum(); jel++)
345 {
346 const auto& dist = myTable.getDistRow(jel);
347 const auto& displ = myTable.getDisplRow(jel);
348 for (int iat = 0; iat < ions.getTotalNum(); iat++)
349 if (nlpp != nullptr && dist[iat] < nlpp->getRmax())
350 Value1 += nlpp->evaluateValueAndDerivatives(elec, iat, psi, jel, dist[iat], -displ[iat], optvars, dlogpsi,
351 dhpsioverpsi);
352 }
353 REQUIRE(Value1 == Approx(6.9015710211e-02));
354
355 REQUIRE(std::real(dhpsioverpsi[0]) == Approx(-0.6379341942));
356 REQUIRE(std::real(dhpsioverpsi[2]) == Approx(1.5269279991));
357 REQUIRE(std::real(dhpsioverpsi[3]) == Approx(-0.0355730676));
358 REQUIRE(std::real(dhpsioverpsi[9]) == Approx(0.279561213));
359 REQUIRE(std::real(dhpsioverpsi[10]) == Approx(-0.3968763604));
360
361 double Value2(0.0);
362 double Value3(0.0);
363 ParticleSet::ParticlePos_t PulayTerm, HFTerm, HFTerm2;
364 HFTerm.resize(ions.getTotalNum());
365 HFTerm2.resize(ions.getTotalNum());
366 PulayTerm.resize(ions.getTotalNum());
367 HFTerm = 0;
368 HFTerm2 = 0;
369 PulayTerm = 0;
370
371 for (int jel = 0; jel < elec.getTotalNum(); jel++)
372 {
373 const auto& dist = myTable.getDistRow(jel);
374 const auto& displ = myTable.getDisplRow(jel);
375 for (int iat = 0; iat < ions.getTotalNum(); iat++)
376 if (nlpp != nullptr && dist[iat] < nlpp->getRmax())
377 {
378 Value2 += nlpp->evaluateOneWithForces(elec, iat, psi, jel, dist[iat], -displ[iat], HFTerm[iat]);
379 Value3 +=
380 nlpp->evaluateOneWithForces(elec, ions, iat, psi, jel, dist[iat], -displ[iat], HFTerm2[iat], PulayTerm);
381 }
382 }
383 //These values are validated against print statements.
384 //Two-body jastrow-only wave functions agree with finite difference of NLPP to machine precision.
385 // These numbers assume the Hellman Feynmann implementation is correct, and dump the values
386 // when a one body term is added in.
387
388 REQUIRE(Value2 == Approx(6.9015710211e-02));
389 REQUIRE(Value3 == Approx(6.9015710211e-02));
390
391 //The total force (HFTerm+PulayTerm) is validated against finite difference of nonlocal PP w.r.t
392 //ion coordinates. delta=1e-6. Should be good up to 7th or 8th sig fig. These are:
393 // F[0][0]= 0.3474359
394 // F[0][1]= 0
395 // F[0][2]= 0
396 // F[1][0]=-0.002734064
397 // F[1][1]= 0
398 // F[1][2]= 0
399
400 REQUIRE(HFTerm[0][0] == Approx(-0.3557369031));
401 REQUIRE(HFTerm[0][1] == Approx(0.0));
402 REQUIRE(HFTerm[0][2] == Approx(0.0));
403 REQUIRE(HFTerm[1][0] == Approx(0.001068673105));
404 REQUIRE(HFTerm[1][1] == Approx(0.0));
405 REQUIRE(HFTerm[1][2] == Approx(0.0));
406
407 REQUIRE(HFTerm2[0][0] == Approx(-0.3557369031));
408 REQUIRE(HFTerm2[0][1] == Approx(0.0));
409 REQUIRE(HFTerm2[0][2] == Approx(0.0));
410 REQUIRE(HFTerm2[1][0] == Approx(0.001068673105));
411 REQUIRE(HFTerm2[1][1] == Approx(0.0));
412 REQUIRE(HFTerm2[1][2] == Approx(0.0));
413
414 REQUIRE(PulayTerm[0][0] == Approx(0.008300993315));
415 REQUIRE(PulayTerm[0][1] == Approx(0.0));
416 REQUIRE(PulayTerm[0][2] == Approx(0.0));
417 REQUIRE(PulayTerm[1][0] == Approx(0.001665391103));
418 REQUIRE(PulayTerm[1][1] == Approx(0.0));
419 REQUIRE(PulayTerm[1][2] == Approx(0.0));
420
421 //Comparing against finite difference results above, here's what we get.
422 //HFTerm[0][0]+PulayTerm[0][0] = −0.34743591
423 //HFTerm[0][1]+PulayTerm[0][1] = 0.0
424 //HFTerm[0][2]+PulayTerm[0][2] = 0.0
425 //HFTerm[1][0]+PulayTerm[1][0] = 0.002734064
426 //HFTerm[1][1]+PulayTerm[1][1] = 0.0
427 //HFTerm[1][2]+PulayTerm[1][2] = 0.0
428 }
429
430 #ifdef QMC_COMPLEX
431 TEST_CASE("Evaluate_soecp", "[hamiltonian]")
432 {
433 app_log() << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n";
434 app_log() << "!!!! Evaluate SOECPComponent !!!!\n";
435 app_log() << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n";
436 typedef QMCTraits::RealType RealType;
437 typedef QMCTraits::ValueType ValueType;
438 typedef QMCTraits::PosType PosType;
439
440 Communicate* c = OHMMS::Controller;
441
442 //Cell definition:
443
444 CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> Lattice;
445 Lattice.BoxBConds = false; // periodic
446 Lattice.R.diagonal(20);
447 Lattice.LR_dim_cutoff = 15;
448 Lattice.reset();
449
450 auto ions_uptr = std::make_unique<ParticleSet>();
451 auto elec_uptr = std::make_unique<ParticleSet>();
452 ParticleSet& ions(*ions_uptr);
453 ParticleSet& elec(*elec_uptr);
454
455 ions.setName("ion0");
456 ions.create(1);
457 ions.R[0][0] = 0.0;
458 ions.R[0][1] = 0.0;
459 ions.R[0][2] = 0.0;
460
461
462 SpeciesSet& ion_species = ions.getSpeciesSet();
463 int pIdx = ion_species.addSpecies("H");
464 int pChargeIdx = ion_species.addAttribute("charge");
465 int iatnumber = ion_species.addAttribute("atomic_number");
466 ion_species(pChargeIdx, pIdx) = 0;
467 ion_species(iatnumber, pIdx) = 1;
468 ions.Lattice = Lattice;
469 ions.createSK();
470
471
472 elec.Lattice = Lattice;
473 elec.setName("e");
474 elec.create(1);
475 elec.R[0][0] = 0.138;
476 elec.R[0][1] = -0.24;
477 elec.R[0][2] = 0.216;
478 elec.spins[0] = 0.0;
479
480 SpeciesSet& tspecies = elec.getSpeciesSet();
481 int upIdx = tspecies.addSpecies("u");
482 int chargeIdx = tspecies.addAttribute("charge");
483 int massIdx = tspecies.addAttribute("mass");
484 tspecies(chargeIdx, upIdx) = -1;
485 tspecies(massIdx, upIdx) = 1.0;
486
487 elec.createSK();
488
489 ParticleSetPool ptcl = ParticleSetPool(c);
490 ptcl.addParticleSet(std::move(elec_uptr));
491 ptcl.addParticleSet(std::move(ions_uptr));
492
493 ions.resetGroups();
494 elec.resetGroups();
495
496 TrialWaveFunction psi;
497
498 std::vector<PosType> kup, kdn;
499 std::vector<RealType> k2up, k2dn;
500 QMCTraits::IndexType nelec = elec.getTotalNum();
501 REQUIRE(nelec == 1);
502
503 kup.resize(nelec);
504 kup[0] = PosType(1, 1, 1);
505
506 k2up.resize(nelec);
507 //For some goofy reason, EGOSet needs to be initialized with:
508 //1.) A k-vector list (fine).
509 //2.) A list of -|k|^2. To save on expensive - sign multiplication apparently.
510 k2up[0] = -dot(kup[0], kup[0]);
511
512 kdn.resize(nelec);
513 kdn[0] = PosType(2, 2, 2);
514
515 k2dn.resize(nelec);
516 k2dn[0] = -dot(kdn[0], kdn[0]);
517
518 auto spo_up = std::make_unique<EGOSet>(kup, k2up);
519 auto spo_dn = std::make_unique<EGOSet>(kdn, k2dn);
520
521 auto spinor_set = std::make_unique<SpinorSet>();
522 spinor_set->set_spos(std::move(spo_up), std::move(spo_dn));
523 QMCTraits::IndexType norb = spinor_set->getOrbitalSetSize();
524 REQUIRE(norb == 1);
525
526 DiracDeterminant<>* dd = new DiracDeterminant<>(std::move(spinor_set));
527 dd->resize(nelec, norb);
528
529 psi.addComponent(dd);
530
531 //Now we set up the SO ECP component.
532 ECPComponentBuilder ecp("test_read_soecp", c);
533
534 bool okay3 = ecp.read_pp_file("so_ecp_test.xml");
535 REQUIRE(okay3);
536
537 SOECPComponent* sopp = ecp.pp_so.get();
538 REQUIRE(sopp != nullptr);
539 copyGridUnrotatedForTest(*sopp);
540
541 const int myTableIndex = elec.addTable(ions);
542
543 const auto& myTable = elec.getDistTable(myTableIndex);
544
545 // update all distance tables
546 ions.update();
547 elec.update();
548
549 //Need to set up temporary data for this configuration in trial wavefunction. Needed for ratios.
550 double logpsi = psi.evaluateLog(elec);
551
552 RealType Value1(0.0);
553
554 for (int jel = 0; jel < elec.getTotalNum(); jel++)
555 {
556 const auto& dist = myTable.getDistRow(jel);
557 const auto& displ = myTable.getDisplRow(jel);
558 for (int iat = 0; iat < ions.getTotalNum(); iat++)
559 {
560 if (sopp != nullptr && dist[iat] < sopp->getRmax())
561 {
562 Value1 += sopp->evaluateOne(elec, iat, psi, jel, dist[iat], RealType(-1) * displ[iat]);
563 }
564 }
565 }
566 REQUIRE(Value1 == Approx(-0.3214176962));
567 }
568 #endif
569
570
571 } // namespace qmcplusplus
572