QMCPACK
test_einset_spinor.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) 2022 QMCPACK developers
6 //
7 // File developed by: Raymond Clay, rclay@sandia.gov, Sandia National Laboratories
8 // Cody A. Melton, cmelton@sandia.gov, Sandia National Laboratories
9 //
10 // File created by: Raymond Clay, rclay@sandia.gov, Sandia National Laboratories
11 //////////////////////////////////////////////////////////////////////////////////////
12 
13 #include "catch.hpp"
14 
15 #include "OhmmsData/Libxml2Doc.h"
16 #include "OhmmsPETE/OhmmsMatrix.h"
17 #include "Particle/ParticleSet.h"
25 
26 #include <stdio.h>
27 #include <string>
28 #include <limits>
29 
30 using std::string;
31 
32 namespace qmcplusplus
33 {
34 //Now we test the spinor set with Einspline orbitals from HDF.
35 #ifdef QMC_COMPLEX
36 TEST_CASE("Einspline SpinorSet from HDF", "[wavefunction]")
37 {
38  app_log() << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n";
39  app_log() << "!!!!! Einspline SpinorSet from HDF !!!!!\n";
40  app_log() << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n";
41 
43  using RealType = SPOSet::RealType;
45 
47  // O2 test example from pwscf non-collinear calculation.
48  lattice.R = {5.10509515, -3.23993545, 0.00000000, 5.10509515, 3.23993545,
49  0.00000000, -6.49690625, 0.00000000, 7.08268015};
50 
51  ParticleSetPool ptcl = ParticleSetPool(c);
52  ptcl.setSimulationCell(lattice);
53  auto ions_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
54  auto elec_uptr = std::make_unique<ParticleSet>(ptcl.getSimulationCell());
55  ParticleSet& ions_(*ions_uptr);
56  ParticleSet& elec_(*elec_uptr);
57 
58  ions_.setName("ion");
59  ptcl.addParticleSet(std::move(ions_uptr));
60  ions_.create({2});
61 
62  ions_.R[0] = {0.00000000, 0.00000000, 1.08659253};
63  ions_.R[1] = {0.00000000, 0.00000000, -1.08659253};
64  elec_.setName("elec");
65  ptcl.addParticleSet(std::move(elec_uptr));
66  elec_.create({3});
67  elec_.R[0] = {0.1, -0.3, 1.0};
68  elec_.R[1] = {-0.1, 0.3, 1.0};
69  elec_.R[2] = {0.1, 0.2, 0.3};
70  elec_.spins[0] = 0.0;
71  elec_.spins[1] = 0.2;
72  elec_.spins[2] = 0.4;
73  elec_.setSpinor(true);
74 
75  SpeciesSet& tspecies = elec_.getSpeciesSet();
76  int upIdx = tspecies.addSpecies("u");
77  int chargeIdx = tspecies.addAttribute("charge");
78  tspecies(chargeIdx, upIdx) = -1;
79 
80  elec_.update();
81  ions_.update();
82 
83 
84  const char* particles = R"(<tmp>
85  <sposet_builder name="A" type="einspline" href="o2_45deg_spins.pwscf.h5" tilematrix="1 0 0 0 1 0 0 0 1" twistnum="0" source="ion" size="3" precision="float" meshfactor="4.0">
86  <sposet name="myspo" size="3">
87  <occupation mode="ground"/>
88  </sposet>
89  </sposet_builder>
90  </tmp>
91 )";
92 
94  bool okay = doc.parseFromString(particles);
95  CHECK(okay);
96 
97  xmlNodePtr root = doc.getRoot();
98 
99  xmlNodePtr ein1 = xmlFirstElementChild(root);
100 
101  SPOSetBuilderFactory fac(c, elec_, ptcl.getPool());
102  const auto spo_builder_ptr = fac.createSPOSetBuilder(ein1);
103  auto& builder = *spo_builder_ptr;
104 
105  auto spo = builder.createSPOSet(ein1);
106  CHECK(spo);
107 
108  SPOSet::ValueMatrix psiM(elec_.R.size(), spo->getOrbitalSetSize());
109  SPOSet::GradMatrix dpsiM(elec_.R.size(), spo->getOrbitalSetSize());
110  SPOSet::ValueMatrix dspsiM(elec_.R.size(), spo->getOrbitalSetSize()); //spin gradient
111  SPOSet::ValueMatrix d2psiM(elec_.R.size(), spo->getOrbitalSetSize());
112 
113  //These are the reference values computed from a spin-polarized calculation,
114  //with the assumption that the coefficients for phi^\uparrow
115  SPOSet::ValueMatrix psiM_up(elec_.R.size(), spo->getOrbitalSetSize());
116  SPOSet::ValueMatrix psiM_down(elec_.R.size(), spo->getOrbitalSetSize());
117  SPOSet::ValueMatrix psiM_ref(elec_.R.size(), spo->getOrbitalSetSize());
118  SPOSet::ValueMatrix dspsiM_ref(elec_.R.size(), spo->getOrbitalSetSize());
119 
120  SPOSet::GradMatrix dpsiM_up(elec_.R.size(), spo->getOrbitalSetSize());
121  SPOSet::GradMatrix dpsiM_down(elec_.R.size(), spo->getOrbitalSetSize());
122  SPOSet::GradMatrix dpsiM_ref(elec_.R.size(), spo->getOrbitalSetSize());
123 
124  SPOSet::ValueMatrix d2psiM_up(elec_.R.size(), spo->getOrbitalSetSize());
125  SPOSet::ValueMatrix d2psiM_down(elec_.R.size(), spo->getOrbitalSetSize());
126  SPOSet::ValueMatrix d2psiM_ref(elec_.R.size(), spo->getOrbitalSetSize());
127 
128 
129  //These reference values were generated as follows:
130  // 1.) Non-Collinear O2 calculation in PBC's performed using Quantum Espresso.
131  // 2.) Spinor wavefunction converted to HDF5 using convertpw4qmc tool. Mainly, this places the up channel in the spin_0 slot, and spin down in spin_1.
132  // 3.) Up and down values for spinor components are directly evaluated using plane waves to make reference values.
133 
134  //reference values, elec 0, orbital 0
135  psiM_up[0][0] = ValueType(3.0546848269585873, 2.6698880339914073);
136  dpsiM_up[0][0][0] = ValueType(-0.1851673255644419, -0.1618419361786101);
137  dpsiM_up[0][0][1] = ValueType(0.60078848166567, 0.5251078699998748);
138  dpsiM_up[0][0][2] = ValueType(-4.882727805715862, -4.267653505749376);
139  d2psiM_up[0][0] = ValueType(-5.949965529898391, -5.200452764159391);
140  psiM_down[0][0] = ValueType(1.26529305723931, 1.105896575232331);
141  dpsiM_down[0][0][0] = ValueType(-0.07669912992155427, -0.0670361472429191);
142  dpsiM_down[0][0][1] = ValueType(0.24885435007410817, 0.21750534167815624);
143  dpsiM_down[0][0][2] = ValueType(-2.0224938024371926, -1.7677088591422807);
144  d2psiM_down[0][0] = ValueType(-2.4645607715652593, -2.1540811306541467);
145 
146  //reference values, elec 0, orbital 1
147  psiM_up[0][1] = ValueType(-0.05877385066521865, -1.652862402896465);
148  dpsiM_up[0][1][0] = ValueType(0.003806703582748908, 0.1070679097294273);
149  dpsiM_up[0][1][1] = ValueType(-0.01231680454381807, -0.34637011425999203);
150  dpsiM_up[0][1][2] = ValueType(0.09976715349787516, 2.8057021736885446);
151  d2psiM_up[0][1] = ValueType(0.11243204605091028, 3.1618388472842645);
152  psiM_down[0][1] = ValueType(0.14190658959478397, 3.9903718588991324);
153  dpsiM_down[0][1][0] = ValueType(-0.009192765914265412, -0.258482966548934);
154  dpsiM_down[0][1][1] = ValueType(0.0297380860827861, 0.8362132640454849);
155  dpsiM_down[0][1][2] = ValueType(-0.24088367433889216, -6.773578493991127);
156  d2psiM_down[0][1] = ValueType(-0.27145824651495754, -7.633377309590959);
157 
158  //reference values, elec 0, orbital 2
159  psiM_up[0][2] = ValueType(1.7341610599591415, 3.26452640029962);
160  dpsiM_up[0][2][0] = ValueType(0.02410454300409052, 0.04537650582197769);
161  dpsiM_up[0][2][1] = ValueType(-0.0812208525894339, -0.15289674365675554);
162  dpsiM_up[0][2][2] = ValueType(2.056759046129918, 3.871811527443685);
163  d2psiM_up[0][2] = ValueType(-0.3715079628152589, -0.6993565364031098);
164  psiM_down[0][2] = ValueType(0.7183159092489255, 1.3522053467852482);
165  dpsiM_down[0][2][0] = ValueType(0.009984891703495969, 0.01879614056538452);
166  dpsiM_down[0][2][1] = ValueType(-0.033643334635896874, -0.06333185576262795);
167  dpsiM_down[0][2][2] = ValueType(0.8519414150210098, 1.603750121892103);
168  d2psiM_down[0][2] = ValueType(-0.15388532808689237, -0.2896807896155573);
169 
170 
171  //reference values, elec 1, orbital 0
172  psiM_up[1][0] = ValueType(3.0526148989188244, 2.6680787492636187);
173  dpsiM_up[1][0][0] = ValueType(0.20449301174627027, 0.17873282170446286);
174  dpsiM_up[1][0][1] = ValueType(-0.6096780888298439, -0.5328779559603193);
175  dpsiM_up[1][0][2] = ValueType(-4.885040183718155, -4.269674660852541);
176  d2psiM_up[1][0] = ValueType(-5.875072106235885, -5.134991421765417);
177  psiM_down[1][0] = ValueType(1.2644358252460446, 1.1051472909800415);
178  dpsiM_down[1][0][0] = ValueType(0.0847027693565092, 0.0740326785381825);
179  dpsiM_down[1][0][1] = ValueType(-0.25253678670740615, -0.22072360765802454);
180  dpsiM_down[1][0][2] = ValueType(-2.023451052801436, -1.768545254958754);
181  d2psiM_down[1][0] = ValueType(-2.433542859102463, -2.126969850545346);
182 
183  //reference values, elec 1, orbital 1
184  psiM_up[1][1] = ValueType(-0.05872929134760467, -1.6516107719315123);
185  dpsiM_up[1][1][0] = ValueType(-0.0042225364734192325, -0.11876835593196035);
186  dpsiM_up[1][1][1] = ValueType(0.012491965861615007, 0.35129150754532346);
187  dpsiM_up[1][1][2] = ValueType(0.09980846579193113, 2.806855260627992);
188  d2psiM_up[1][1] = ValueType(0.11086616211845124, 3.1178291585160025);
189  psiM_down[1][1] = ValueType(0.14179908178203693, 3.9873502499791);
190  dpsiM_down[1][1][0] = ValueType(0.010197668920898767, 0.2867312658960351);
191  dpsiM_down[1][1][1] = ValueType(-0.030160592987572725, -0.8480940968707702);
192  dpsiM_down[1][1][2] = ValueType(-0.24098310461934494, -6.776362513721667);
193  d2psiM_down[1][1] = ValueType(-0.26767894399782044, -7.527130337670782);
194 
195  //reference values, elec 1, orbital 2
196  psiM_up[1][2] = ValueType(1.7338733833257558, 3.263984881354726);
197  dpsiM_up[1][2][0] = ValueType(-0.02165584086872901, -0.0407670903699481);
198  dpsiM_up[1][2][1] = ValueType(0.08288083949305346, 0.15602164581174188);
199  dpsiM_up[1][2][2] = ValueType(2.0621151061966456, 3.881894235760205);
200  d2psiM_up[1][2] = ValueType(-0.3566890854259599, -0.6714605501817572);
201  psiM_down[1][2] = ValueType(0.7181968101586865, 1.3519810682722548);
202  dpsiM_down[1][2][0] = ValueType(-0.00897085696147509, -0.01688677968381685);
203  dpsiM_down[1][2][1] = ValueType(0.03433096009876233, 0.06462627360298884);
204  dpsiM_down[1][2][2] = ValueType(0.8541600134552085, 1.6079267140278);
205  d2psiM_down[1][2] = ValueType(-0.1477482607270697, -0.2781267037329471);
206 
207 
208  //reference values, elec 2, orbital 0
209  psiM_up[2][0] = ValueType(4.774972481925916, 4.173472045741891);
210  dpsiM_up[2][0][0] = ValueType(-0.9457258313862555, -0.8265932416325468);
211  dpsiM_up[2][0][1] = ValueType(-1.899691502097708, -1.6603886891267976);
212  dpsiM_up[2][0][2] = ValueType(0.7771301258291673, 0.6792356362625406);
213  d2psiM_up[2][0] = ValueType(-20.40945585069239, -17.83848971293892);
214  psiM_down[2][0] = ValueType(1.9778599493693798, 1.7286974948556129);
215  dpsiM_down[2][0][0] = ValueType(-0.3917329590645693, -0.3423838452062018);
216  dpsiM_down[2][0][1] = ValueType(-0.786878588427934, -0.6877513351266712);
217  dpsiM_down[2][0][2] = ValueType(0.3218978428488249, 0.281346363058232);
218  d2psiM_down[2][0] = ValueType(-8.45387947982597, -7.388894117402044);
219 
220  //reference values, elec 2, orbital 1
221  psiM_up[2][1] = ValueType(-0.095146182382511, -2.6757440636563343);
222  dpsiM_up[2][1][0] = ValueType(0.01912387482485274, 0.53780199541144);
223  dpsiM_up[2][1][1] = ValueType(0.03838799057297392, 1.0795586887258484);
224  dpsiM_up[2][1][2] = ValueType(-0.013683016882420245, -0.38479709783829663);
225  d2psiM_up[2][1] = ValueType(0.41702609987278866, 11.727776988772089);
226  psiM_down[2][1] = ValueType(0.22972652344682393, 6.459831671158625);
227  dpsiM_down[2][1][0] = ValueType(-0.046172654628017486, -1.2983723819140731);
228  dpsiM_down[2][1][1] = ValueType(-0.09268554947869961, -2.606290867185097);
229  dpsiM_down[2][1][2] = ValueType(0.03303644631176311, 0.9289838072933512);
230  d2psiM_down[2][1] = ValueType(-1.006891760427076, -28.313415815931304);
231 
232  //reference values, elec 2, orbital 2
233  psiM_up[2][2] = ValueType(0.5573944761518197, 1.0492847452220198);
234  dpsiM_up[2][2][0] = ValueType(0.06369314000215545, 0.11990123738728313);
235  dpsiM_up[2][2][1] = ValueType(0.1265010825081423, 0.23813600324436654);
236  dpsiM_up[2][2][2] = ValueType(1.6373025933118952, 3.082192181081695);
237  d2psiM_up[2][2] = ValueType(-0.8650133588132842, -1.6283710474610622);
238  psiM_down[2][2] = ValueType(0.23088102734804172, 0.43462602449526555);
239  dpsiM_down[2][2][0] = ValueType(0.02638262680056288, 0.0496645026299118);
240  dpsiM_down[2][2][1] = ValueType(0.052398844163716374, 0.0986386889965079);
241  dpsiM_down[2][2][2] = ValueType(0.6781955750070273, 1.276680192776833);
242  d2psiM_down[2][2] = ValueType(-0.3583001666840082, -0.6744889374649511);
243 
244 
245  RealType h = 0.001;
246  RealType h2 = 0.1;
247  for (unsigned int iat = 0; iat < 3; iat++)
248  {
249  RealType s = elec_.spins[iat];
250  RealType coss(0.0), sins(0.0);
251 
252  coss = std::cos(s);
253  sins = std::sin(s);
254 
255  ValueType eis(coss, sins);
256  ValueType emis(coss, -sins);
257  ValueType eye(0, 1.0);
258  //Using the reference values for the up and down channels invdividually, we build the total reference spinor value
259  //consistent with the current spin value of particle iat.
260  psiM_ref[iat][0] = eis * psiM_up[iat][0] + emis * psiM_down[iat][0];
261  psiM_ref[iat][1] = eis * psiM_up[iat][1] + emis * psiM_down[iat][1];
262  psiM_ref[iat][2] = eis * psiM_up[iat][2] + emis * psiM_down[iat][2];
263 
264  dspsiM_ref[iat][0] = eye * (eis * psiM_up[iat][0] - emis * psiM_down[iat][0]);
265  dspsiM_ref[iat][1] = eye * (eis * psiM_up[iat][1] - emis * psiM_down[iat][1]);
266  dspsiM_ref[iat][2] = eye * (eis * psiM_up[iat][2] - emis * psiM_down[iat][2]);
267 
268  dpsiM_ref[iat][0] = eis * dpsiM_up[iat][0] + emis * dpsiM_down[iat][0];
269  dpsiM_ref[iat][1] = eis * dpsiM_up[iat][1] + emis * dpsiM_down[iat][1];
270  dpsiM_ref[iat][2] = eis * dpsiM_up[iat][2] + emis * dpsiM_down[iat][2];
271 
272  d2psiM_ref[iat][0] = eis * d2psiM_up[iat][0] + emis * d2psiM_down[iat][0];
273  d2psiM_ref[iat][1] = eis * d2psiM_up[iat][1] + emis * d2psiM_down[iat][1];
274  d2psiM_ref[iat][2] = eis * d2psiM_up[iat][2] + emis * d2psiM_down[iat][2];
275  }
276 
277  //OK. Going to test evaluate_notranspose with spin.
278  spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM, dpsiM, d2psiM);
279 
280  for (unsigned int iat = 0; iat < 3; iat++)
281  {
282  CHECK(psiM[iat][0] == ComplexApprox(psiM_ref[iat][0]).epsilon(h));
283  CHECK(psiM[iat][1] == ComplexApprox(psiM_ref[iat][1]).epsilon(h));
284  CHECK(psiM[iat][2] == ComplexApprox(psiM_ref[iat][2]).epsilon(h));
285 
286  CHECK(dpsiM[iat][0][0] == ComplexApprox(dpsiM_ref[iat][0][0]).epsilon(h));
287  CHECK(dpsiM[iat][0][1] == ComplexApprox(dpsiM_ref[iat][0][1]).epsilon(h));
288  CHECK(dpsiM[iat][0][2] == ComplexApprox(dpsiM_ref[iat][0][2]).epsilon(h));
289 
290  CHECK(dpsiM[iat][1][0] == ComplexApprox(dpsiM_ref[iat][1][0]).epsilon(h));
291  CHECK(dpsiM[iat][1][1] == ComplexApprox(dpsiM_ref[iat][1][1]).epsilon(h));
292  CHECK(dpsiM[iat][1][2] == ComplexApprox(dpsiM_ref[iat][1][2]).epsilon(h));
293 
294  CHECK(dpsiM[iat][2][0] == ComplexApprox(dpsiM_ref[iat][2][0]).epsilon(h));
295  CHECK(dpsiM[iat][2][1] == ComplexApprox(dpsiM_ref[iat][2][1]).epsilon(h));
296  CHECK(dpsiM[iat][2][2] == ComplexApprox(dpsiM_ref[iat][2][2]).epsilon(h));
297 
298 
299  CHECK(d2psiM[iat][0] == ComplexApprox(d2psiM_ref[iat][0]).epsilon(h2));
300  CHECK(d2psiM[iat][1] == ComplexApprox(d2psiM_ref[iat][1]).epsilon(h2));
301  CHECK(d2psiM[iat][2] == ComplexApprox(d2psiM_ref[iat][2]).epsilon(h2));
302  }
303 
304  //Now we're going to test evaluateValue and evaluateVGL:
305 
306  int OrbitalSetSize = spo->getOrbitalSetSize();
307  //temporary arrays for holding the values of the up and down channels respectively.
308  SPOSet::ValueVector psi_work;
309  SPOSet::ValueVector psi_work_up;
310  SPOSet::ValueVector psi_work_down;
311 
312  //temporary arrays for holding the gradients of the up and down channels respectively.
313  SPOSet::GradVector dpsi_work;
314 
315  //temporary arrays for holding the laplacians of the up and down channels respectively.
316  SPOSet::ValueVector d2psi_work;
317 
318  //temporary arrays for holding the spin gradient
319  SPOSet::ValueVector dspsi_work;
320 
321 
322  psi_work.resize(OrbitalSetSize);
323  psi_work_up.resize(OrbitalSetSize);
324  psi_work_down.resize(OrbitalSetSize);
325 
326  dpsi_work.resize(OrbitalSetSize);
327 
328  d2psi_work.resize(OrbitalSetSize);
329 
330  dspsi_work.resize(OrbitalSetSize);
331 
332  //We worked hard to generate nice reference data above. Let's generate a test for evaluateV
333  //and evaluateVGL by perturbing the electronic configuration by dR, and then make
334  //single particle moves that bring it back to our original R reference values.
335 
336  //Our perturbation vector.
338  dR.resize(3);
339  //Values chosen based on divine inspiration. Not important.
340  dR[0][0] = 0.1;
341  dR[0][1] = 0.2;
342  dR[0][2] = 0.1;
343  dR[1][0] = -0.1;
344  dR[1][1] = -0.05;
345  dR[1][2] = 0.05;
346  dR[2][0] = -0.1;
347  dR[2][1] = 0.0;
348  dR[2][2] = 0.0;
349 
350  //The new R of our perturbed particle set.
352  Rnew.resize(3);
353  Rnew = elec_.R + dR;
354  elec_.R = Rnew;
355  elec_.update();
356 
357  //Now we test evaluateValue()
358  VirtualParticleSet vp(elec_, 1);
359  for (unsigned int iat = 0; iat < 3; iat++)
360  {
361  std::vector<ParticleSet::PosType> delta_v(1);
362  delta_v[0] = -dR[iat];
363  vp.makeMoves(elec_, iat, delta_v);
364 
365  psi_work = 0.0;
366  elec_.makeMove(iat, -dR[iat], false);
367  spo->evaluateValue(elec_, iat, psi_work);
368 
369  CHECK(psi_work[0] == ComplexApprox(psiM_ref[iat][0]).epsilon(h));
370  CHECK(psi_work[1] == ComplexApprox(psiM_ref[iat][1]).epsilon(h));
371  CHECK(psi_work[2] == ComplexApprox(psiM_ref[iat][2]).epsilon(h));
372 
373 
374  std::vector<ValueType> ratios(1);
375  SPOSet::ValueVector psi(OrbitalSetSize);
376  SPOSet::ValueVector psiinv(OrbitalSetSize);
377  psiinv = 2.1;
378  std::pair<SPOSet::ValueVector, SPOSet::ValueVector> multiplier;
379  auto& up_factor = multiplier.first;
380  auto& dn_factor = multiplier.second;
381  up_factor.resize(OrbitalSetSize);
382  dn_factor.resize(OrbitalSetSize);
383  up_factor = 5.2;
384  dn_factor = -0.3;
385  spo->evaluateDetSpinorRatios(vp, psi, multiplier, psiinv, ratios);
386  //create reference value
387  RealType s = elec_.spins[iat];
388  RealType coss(0.0), sins(0.0);
389  coss = std::cos(s);
390  sins = std::sin(s);
391  ValueType eis(coss, sins);
392  ValueType emis(coss, -sins);
393 
394  ValueType refVal;
395  for (int iorb = 0; iorb < OrbitalSetSize; iorb++)
396  refVal += static_cast<ValueType>(2.1) *
397  (eis * static_cast<ValueType>(5.2) * psiM_up[iat][iorb] +
398  emis * static_cast<ValueType>(-0.3) * psiM_down[iat][iorb]);
399  CHECK(ratios[0] == ComplexApprox(refVal));
400  elec_.rejectMove(iat);
401  }
402 
403  //Now we test evaluateVGL()
404  for (unsigned int iat = 0; iat < 3; iat++)
405  {
406  psi_work = 0.0;
407  dpsi_work = 0.0;
408  d2psi_work = 0.0;
409  dspsi_work = 0.0;
410 
411  elec_.makeMove(iat, -dR[iat], false);
412  spo->evaluateVGL_spin(elec_, iat, psi_work, dpsi_work, d2psi_work, dspsi_work);
413 
414  CHECK(psi_work[0] == ComplexApprox(psiM_ref[iat][0]).epsilon(h));
415  CHECK(psi_work[1] == ComplexApprox(psiM_ref[iat][1]).epsilon(h));
416  CHECK(psi_work[2] == ComplexApprox(psiM_ref[iat][2]).epsilon(h));
417 
418  CHECK(dpsi_work[0][0] == ComplexApprox(dpsiM_ref[iat][0][0]).epsilon(h));
419  CHECK(dpsi_work[0][1] == ComplexApprox(dpsiM_ref[iat][0][1]).epsilon(h));
420  CHECK(dpsi_work[0][2] == ComplexApprox(dpsiM_ref[iat][0][2]).epsilon(h));
421 
422  CHECK(dpsi_work[1][0] == ComplexApprox(dpsiM_ref[iat][1][0]).epsilon(h));
423  CHECK(dpsi_work[1][1] == ComplexApprox(dpsiM_ref[iat][1][1]).epsilon(h));
424  CHECK(dpsi_work[1][2] == ComplexApprox(dpsiM_ref[iat][1][2]).epsilon(h));
425 
426  CHECK(dpsi_work[2][0] == ComplexApprox(dpsiM_ref[iat][2][0]).epsilon(h));
427  CHECK(dpsi_work[2][1] == ComplexApprox(dpsiM_ref[iat][2][1]).epsilon(h));
428  CHECK(dpsi_work[2][2] == ComplexApprox(dpsiM_ref[iat][2][2]).epsilon(h));
429 
430  CHECK(d2psi_work[0] == ComplexApprox(d2psiM_ref[iat][0]).epsilon(h2));
431  CHECK(d2psi_work[1] == ComplexApprox(d2psiM_ref[iat][1]).epsilon(h2));
432  CHECK(d2psi_work[2] == ComplexApprox(d2psiM_ref[iat][2]).epsilon(h2));
433 
434  CHECK(dspsi_work[0] == ComplexApprox(dspsiM_ref[iat][0]).epsilon(h));
435  CHECK(dspsi_work[1] == ComplexApprox(dspsiM_ref[iat][1]).epsilon(h));
436  CHECK(dspsi_work[2] == ComplexApprox(dspsiM_ref[iat][2]).epsilon(h));
437 
438  elec_.rejectMove(iat);
439  }
440 
441  //Now we test evaluateSpin:
442 
443  for (unsigned int iat = 0; iat < 3; iat++)
444  {
445  psi_work = 0.0;
446  dspsi_work = 0.0;
447 
448  elec_.makeMove(iat, -dR[iat], false);
449  spo->evaluate_spin(elec_, iat, psi_work, dspsi_work);
450 
451  CHECK(psi_work[0] == ComplexApprox(psiM_ref[iat][0]).epsilon(h));
452  CHECK(psi_work[1] == ComplexApprox(psiM_ref[iat][1]).epsilon(h));
453  CHECK(psi_work[2] == ComplexApprox(psiM_ref[iat][2]).epsilon(h));
454 
455  CHECK(dspsi_work[0] == ComplexApprox(dspsiM_ref[iat][0]).epsilon(h));
456  CHECK(dspsi_work[1] == ComplexApprox(dspsiM_ref[iat][1]).epsilon(h));
457  CHECK(dspsi_work[2] == ComplexApprox(dspsiM_ref[iat][2]).epsilon(h));
458 
459  elec_.rejectMove(iat);
460  }
461 
462 
463  // test batched interface
464  // first move elec_ back to original positions for reference
465  Rnew = elec_.R - dR;
466  elec_.R = Rnew;
467  elec_.update();
468 
469  //just going to use zero displacements for particle spins in test
471  dS.resize(3);
472 
473  //now create second walker, with permuted particle positions
474  ParticleSet elec_2(elec_);
475  // permute electrons
476  elec_2.R[0] = elec_.R[1];
477  elec_2.R[1] = elec_.R[2];
478  elec_2.R[2] = elec_.R[0];
479  elec_2.spins[0] = elec_.spins[1];
480  elec_2.spins[1] = elec_.spins[2];
481  elec_2.spins[2] = elec_.spins[0];
482 
483  ResourceCollection pset_res("test_pset_res");
484  elec_.createResource(pset_res);
485 
486  RefVectorWithLeader<ParticleSet> p_list(elec_);
487  p_list.push_back(elec_);
488  p_list.push_back(elec_2);
489 
490  ResourceCollectionTeamLock<ParticleSet> mw_pset_lock(pset_res, p_list);
491 
492  //update all walkers
493  elec_.mw_update(p_list);
494 
495  std::unique_ptr<SPOSet> spo_2(spo->makeClone());
496  RefVectorWithLeader<SPOSet> spo_list(*spo);
497  spo_list.push_back(*spo);
498  spo_list.push_back(*spo_2);
499 
500  //test resource APIs
501  //First resource is created, and then passed to the collection so it should be null
502  ResourceCollection spo_res("test_spo_res");
503  spo->createResource(spo_res);
504  SpinorSet& spinor = spo_list.getCastedLeader<SpinorSet>();
505  REQUIRE(!spinor.isResourceOwned());
506  //team lock calls the acquireResource, so now the leader's resource shouldn't be null
507  ResourceCollectionTeamLock<SPOSet> mw_spo_lock(spo_res, spo_list);
508  REQUIRE(spinor.isResourceOwned());
509 
510  SPOSet::ValueMatrix psiM_2(elec_.R.size(), spo->getOrbitalSetSize());
511  SPOSet::GradMatrix dpsiM_2(elec_.R.size(), spo->getOrbitalSetSize());
512  SPOSet::ValueMatrix d2psiM_2(elec_.R.size(), spo->getOrbitalSetSize());
513 
514  RefVector<SPOSet::ValueMatrix> logdet_list;
515  RefVector<SPOSet::GradMatrix> dlogdet_list;
516  RefVector<SPOSet::ValueMatrix> d2logdet_list;
517 
518  logdet_list.push_back(psiM);
519  logdet_list.push_back(psiM_2);
520  dlogdet_list.push_back(dpsiM);
521  dlogdet_list.push_back(dpsiM_2);
522  d2logdet_list.push_back(d2psiM);
523  d2logdet_list.push_back(d2psiM_2);
524 
525  spo->mw_evaluate_notranspose(spo_list, p_list, 0, 3, logdet_list, dlogdet_list, d2logdet_list);
526  for (unsigned int iat = 0; iat < 3; iat++)
527  {
528  //walker 0
529  CHECK(logdet_list[0].get()[iat][0] == ComplexApprox(psiM_ref[iat][0]).epsilon(h));
530  CHECK(logdet_list[0].get()[iat][1] == ComplexApprox(psiM_ref[iat][1]).epsilon(h));
531  CHECK(logdet_list[0].get()[iat][2] == ComplexApprox(psiM_ref[iat][2]).epsilon(h));
532  CHECK(dlogdet_list[0].get()[iat][0][0] == ComplexApprox(dpsiM_ref[iat][0][0]).epsilon(h));
533  CHECK(dlogdet_list[0].get()[iat][0][1] == ComplexApprox(dpsiM_ref[iat][0][1]).epsilon(h));
534  CHECK(dlogdet_list[0].get()[iat][0][2] == ComplexApprox(dpsiM_ref[iat][0][2]).epsilon(h));
535  CHECK(dlogdet_list[0].get()[iat][1][0] == ComplexApprox(dpsiM_ref[iat][1][0]).epsilon(h));
536  CHECK(dlogdet_list[0].get()[iat][1][1] == ComplexApprox(dpsiM_ref[iat][1][1]).epsilon(h));
537  CHECK(dlogdet_list[0].get()[iat][1][2] == ComplexApprox(dpsiM_ref[iat][1][2]).epsilon(h));
538  CHECK(dlogdet_list[0].get()[iat][2][0] == ComplexApprox(dpsiM_ref[iat][2][0]).epsilon(h));
539  CHECK(dlogdet_list[0].get()[iat][2][1] == ComplexApprox(dpsiM_ref[iat][2][1]).epsilon(h));
540  CHECK(dlogdet_list[0].get()[iat][2][2] == ComplexApprox(dpsiM_ref[iat][2][2]).epsilon(h));
541  CHECK(d2logdet_list[0].get()[iat][0] == ComplexApprox(d2psiM_ref[iat][0]).epsilon(h2));
542  CHECK(d2logdet_list[0].get()[iat][1] == ComplexApprox(d2psiM_ref[iat][1]).epsilon(h2));
543  CHECK(d2logdet_list[0].get()[iat][2] == ComplexApprox(d2psiM_ref[iat][2]).epsilon(h2));
544 
545  //walker 1, permuted from reference
546  CHECK(logdet_list[1].get()[iat][0] == ComplexApprox(psiM_ref[(iat + 1) % 3][0]).epsilon(h));
547  CHECK(logdet_list[1].get()[iat][1] == ComplexApprox(psiM_ref[(iat + 1) % 3][1]).epsilon(h));
548  CHECK(logdet_list[1].get()[iat][2] == ComplexApprox(psiM_ref[(iat + 1) % 3][2]).epsilon(h));
549  CHECK(dlogdet_list[1].get()[iat][0][0] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][0][0]).epsilon(h));
550  CHECK(dlogdet_list[1].get()[iat][0][1] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][0][1]).epsilon(h));
551  CHECK(dlogdet_list[1].get()[iat][0][2] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][0][2]).epsilon(h));
552  CHECK(dlogdet_list[1].get()[iat][1][0] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][1][0]).epsilon(h));
553  CHECK(dlogdet_list[1].get()[iat][1][1] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][1][1]).epsilon(h));
554  CHECK(dlogdet_list[1].get()[iat][1][2] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][1][2]).epsilon(h));
555  CHECK(dlogdet_list[1].get()[iat][2][0] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][2][0]).epsilon(h));
556  CHECK(dlogdet_list[1].get()[iat][2][1] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][2][1]).epsilon(h));
557  CHECK(dlogdet_list[1].get()[iat][2][2] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][2][2]).epsilon(h));
558  CHECK(d2logdet_list[1].get()[iat][0] == ComplexApprox(d2psiM_ref[(iat + 1) % 3][0]).epsilon(h2));
559  CHECK(d2logdet_list[1].get()[iat][1] == ComplexApprox(d2psiM_ref[(iat + 1) % 3][1]).epsilon(h2));
560  CHECK(d2logdet_list[1].get()[iat][2] == ComplexApprox(d2psiM_ref[(iat + 1) % 3][2]).epsilon(h2));
561  }
562 
563  //first, lets displace all the electrons in each walker.
564  for (int iat = 0; iat < 3; iat++)
565  {
566  MCCoords<CoordsType::POS_SPIN> displs(2);
567  displs.positions = {dR[iat], dR[iat]};
568  displs.spins = {dS[iat], dS[iat]};
569 
570  elec_.mw_makeMove(p_list, iat, displs);
571  std::vector<bool> accept = {true, true};
572  elec_.mw_accept_rejectMove<CoordsType::POS_SPIN>(p_list, iat, accept);
573  }
574  elec_.mw_update(p_list);
575 
576  SPOSet::ValueVector psi_work_2(OrbitalSetSize);
577  SPOSet::GradVector dpsi_work_2(OrbitalSetSize);
578  SPOSet::ValueVector d2psi_work_2(OrbitalSetSize);
579 
580  RefVector<SPOSet::ValueVector> psi_v_list = {psi_work, psi_work_2};
581  RefVector<SPOSet::GradVector> dpsi_v_list = {dpsi_work, dpsi_work_2};
582  RefVector<SPOSet::ValueVector> d2psi_v_list = {d2psi_work, d2psi_work_2};
583  SPOSet::OffloadMatrix<SPOSet::ComplexType> mw_dspin;
584  mw_dspin.resize(2, OrbitalSetSize);
585  //check mw_evaluateVGLWithSpin
586  for (int iat = 0; iat < 3; iat++)
587  {
588  //reset values to zero, updates the ref vectors to zero as well
589  psi_work = 0.0;
590  dpsi_work = 0.0;
591  d2psi_work = 0.0;
592  psi_work_2 = 0.0;
593  dpsi_work_2 = 0.0;
594  d2psi_work_2 = 0.0;
595 
596  MCCoords<CoordsType::POS_SPIN> displs(2);
597  displs.positions = {-dR[iat], -dR[iat]};
598  displs.spins = {-dS[iat], -dS[iat]};
599 
600  elec_.mw_makeMove(p_list, iat, displs);
601  spo->mw_evaluateVGLWithSpin(spo_list, p_list, iat, psi_v_list, dpsi_v_list, d2psi_v_list, mw_dspin);
602  //walker 0
603  CHECK(psi_v_list[0].get()[0] == ComplexApprox(psiM_ref[iat][0]).epsilon(h));
604  CHECK(psi_v_list[0].get()[1] == ComplexApprox(psiM_ref[iat][1]).epsilon(h));
605  CHECK(psi_v_list[0].get()[2] == ComplexApprox(psiM_ref[iat][2]).epsilon(h));
606 
607  CHECK(dpsi_v_list[0].get()[0][0] == ComplexApprox(dpsiM_ref[iat][0][0]).epsilon(h));
608  CHECK(dpsi_v_list[0].get()[0][1] == ComplexApprox(dpsiM_ref[iat][0][1]).epsilon(h));
609  CHECK(dpsi_v_list[0].get()[0][2] == ComplexApprox(dpsiM_ref[iat][0][2]).epsilon(h));
610 
611  CHECK(dpsi_v_list[0].get()[1][0] == ComplexApprox(dpsiM_ref[iat][1][0]).epsilon(h));
612  CHECK(dpsi_v_list[0].get()[1][1] == ComplexApprox(dpsiM_ref[iat][1][1]).epsilon(h));
613  CHECK(dpsi_v_list[0].get()[1][2] == ComplexApprox(dpsiM_ref[iat][1][2]).epsilon(h));
614 
615  CHECK(dpsi_v_list[0].get()[2][0] == ComplexApprox(dpsiM_ref[iat][2][0]).epsilon(h));
616  CHECK(dpsi_v_list[0].get()[2][1] == ComplexApprox(dpsiM_ref[iat][2][1]).epsilon(h));
617  CHECK(dpsi_v_list[0].get()[2][2] == ComplexApprox(dpsiM_ref[iat][2][2]).epsilon(h));
618 
619  CHECK(d2psi_v_list[0].get()[0] == ComplexApprox(d2psiM_ref[iat][0]).epsilon(h2));
620  CHECK(d2psi_v_list[0].get()[1] == ComplexApprox(d2psiM_ref[iat][1]).epsilon(h2));
621  CHECK(d2psi_v_list[0].get()[2] == ComplexApprox(d2psiM_ref[iat][2]).epsilon(h2));
622 
623  CHECK(mw_dspin[0][0] == ComplexApprox(dspsiM_ref[iat][0]).epsilon(h));
624  CHECK(mw_dspin[0][1] == ComplexApprox(dspsiM_ref[iat][1]).epsilon(h));
625  CHECK(mw_dspin[0][2] == ComplexApprox(dspsiM_ref[iat][2]).epsilon(h));
626 
627  //walker 1, permuted from reference
628  CHECK(psi_v_list[1].get()[0] == ComplexApprox(psiM_ref[(iat + 1) % 3][0]).epsilon(h));
629  CHECK(psi_v_list[1].get()[1] == ComplexApprox(psiM_ref[(iat + 1) % 3][1]).epsilon(h));
630  CHECK(psi_v_list[1].get()[2] == ComplexApprox(psiM_ref[(iat + 1) % 3][2]).epsilon(h));
631 
632  CHECK(dpsi_v_list[1].get()[0][0] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][0][0]).epsilon(h));
633  CHECK(dpsi_v_list[1].get()[0][1] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][0][1]).epsilon(h));
634  CHECK(dpsi_v_list[1].get()[0][2] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][0][2]).epsilon(h));
635 
636  CHECK(dpsi_v_list[1].get()[1][0] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][1][0]).epsilon(h));
637  CHECK(dpsi_v_list[1].get()[1][1] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][1][1]).epsilon(h));
638  CHECK(dpsi_v_list[1].get()[1][2] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][1][2]).epsilon(h));
639 
640  CHECK(dpsi_v_list[1].get()[2][0] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][2][0]).epsilon(h));
641  CHECK(dpsi_v_list[1].get()[2][1] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][2][1]).epsilon(h));
642  CHECK(dpsi_v_list[1].get()[2][2] == ComplexApprox(dpsiM_ref[(iat + 1) % 3][2][2]).epsilon(h));
643 
644  CHECK(d2psi_v_list[1].get()[0] == ComplexApprox(d2psiM_ref[(iat + 1) % 3][0]).epsilon(h2));
645  CHECK(d2psi_v_list[1].get()[1] == ComplexApprox(d2psiM_ref[(iat + 1) % 3][1]).epsilon(h2));
646  CHECK(d2psi_v_list[1].get()[2] == ComplexApprox(d2psiM_ref[(iat + 1) % 3][2]).epsilon(h2));
647 
648  CHECK(mw_dspin[1][0] == ComplexApprox(dspsiM_ref[(iat + 1) % 3][0]).epsilon(h));
649  CHECK(mw_dspin[1][1] == ComplexApprox(dspsiM_ref[(iat + 1) % 3][1]).epsilon(h));
650  CHECK(mw_dspin[1][2] == ComplexApprox(dspsiM_ref[(iat + 1) % 3][2]).epsilon(h));
651 
652  std::vector<bool> accept = {false, false};
653  elec_.mw_accept_rejectMove<CoordsType::POS_SPIN>(p_list, iat, accept);
654  }
655 
656  //Need to remember that the electrons are distored from initial positions
657  //Move them back for rotation test
658  Rnew = elec_.R - dR;
659  elec_.R = Rnew;
660  elec_.update();
661 
662  //Let's also test orbital rotation
663  //random unitary
664  SPOSet::ValueVector row0 = {ValueType(0.47586834, 0.07134236), ValueType(-0.35882226, -0.6570473),
665  ValueType(-0.30699602, -0.3372662)};
666  SPOSet::ValueVector row1 = {ValueType(-0.09062269, -0.23061815), ValueType(0.52931815, 0.07182243),
667  ValueType(-0.79260564, -0.15825018)};
668  SPOSet::ValueVector row2 = {ValueType(0.3187019, -0.7781333), ValueType(0.315927, -0.23321542),
669  ValueType(0.36577135, 0.07035348)};
670  SPOSet::ValueMatrix rot_mat(3, 3);
671  for (int iorb = 0; iorb < 3; iorb++)
672  {
673  rot_mat(0, iorb) = row0[iorb];
674  rot_mat(1, iorb) = row1[iorb];
675  rot_mat(2, iorb) = row2[iorb];
676  }
677  SPOSet::ValueMatrix psiM_rot_manual(elec_.R.size(), spo->size());
678  for (int i = 0; i < elec_.R.size(); i++)
679  for (int j = 0; j < spo->size(); j++)
680  {
681  psiM_rot_manual[i][j] = 0.;
682  for (int k = 0; k < spo->size(); k++)
683  psiM_rot_manual[i][j] += psiM_ref[i][k] * rot_mat[k][j];
684  }
685 
686  spo->storeParamsBeforeRotation();
687  spo->applyRotation(rot_mat, false);
688  spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM, dpsiM, d2psiM);
689  auto check = checkMatrix(psiM_rot_manual, psiM, true);
690  CHECKED_ELSE(check.result) { FAIL(check.result_message); }
691 }
692 
693 #endif //QMC_COMPLEX
694 
695 
696 } // namespace qmcplusplus
class that handles xmlDoc
Definition: Libxml2Doc.h:76
helper functions for EinsplineSetBuilder
Definition: Configuration.h:43
QTBase::RealType RealType
Definition: Configuration.h:58
std::ostream & app_log()
Definition: OutputManager.h:65
CrystalLattice< OHMMS_PRECISION, OHMMS_DIM > ParticleLayout
Definition: Configuration.h:79
MakeReturn< UnaryNode< FnSin, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t sin(const Vector< T1, C1 > &l)
TEST_CASE("complex_helper", "[type_traits]")
xmlNodePtr getRoot()
Definition: Libxml2Doc.h:88
OrbitalSetTraits< ValueType >::ValueMatrix ValueMatrix
Definition: SPOSet.h:50
CHECKED_ELSE(check_matrix_result.result)
Communicate * Controller
Global Communicator for a process.
Definition: Communicate.cpp:35
A proxy class to the quantum ParticleSet.
OrbitalSetTraits< ValueType >::GradMatrix GradMatrix
Definition: SPOSet.h:52
Wrapping information on parallelism.
Definition: Communicate.h:68
CrystalLattice< OHMMS_PRECISION, OHMMS_DIM > lattice
MakeReturn< UnaryNode< FnCos, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t cos(const Vector< T1, C1 > &l)
QTBase::ValueType ValueType
Definition: Configuration.h:60
REQUIRE(std::filesystem::exists(filename))
OrbitalSetTraits< ValueType >::ValueVector ValueVector
Definition: SPOSet.h:49
QMCTraits::RealType RealType
OrbitalSetTraits< ValueType >::GradVector GradVector
Definition: SPOSet.h:51
bool parseFromString(const std::string_view data)
Definition: Libxml2Doc.cpp:204
CheckMatrixResult checkMatrix(M1 &a_mat, M2 &b_mat, const bool check_all=false, std::optional< const double > eps=std::nullopt)
This function checks equality a_mat and b_mat elements M1, M2 need to have their element type declare...
Definition: checkMatrix.hpp:63
CHECK(log_values[0]==ComplexApprox(std::complex< double >{ 5.603777579195571, -6.1586603331188225 }))
ParticleAttrib< SingleParticlePos > ParticlePos
Definition: Configuration.h:92
LatticeGaussianProduct::ValueType ValueType
Declaration of WaveFunctionComponent.
ParticleAttrib< Scalar_t > ParticleScalar
Definition: Configuration.h:91
Declaration of ParticleSetPool.