1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
8 // Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
9 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
10 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
11 // Mark Dewing, markdewing@gmail.com, University of Illinois at Urbana-Champaign
12 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
13 //
14 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
15 //////////////////////////////////////////////////////////////////////////////////////
16
17
18 #include "Particle/MCWalkerConfiguration.h"
19 #include "ParticleBase/ParticleUtility.h"
20 #include "ParticleBase/RandomSeqGenerator.h"
21 #include "Message/Communicate.h"
22 #include "WaveFunctionTester.h"
23 #include "QMCDrivers/DriftOperators.h"
24 #include "LongRange/StructFact.h"
25 #include "OhmmsData/AttributeSet.h"
26 #include "OhmmsData/ParameterSet.h"
27 #include "QMCWaveFunctions/SPOSet.h"
28 #include "QMCWaveFunctions/Fermion/SlaterDet.h"
29 #include "QMCWaveFunctions/OrbitalSetTraits.h"
30 #include "Numerics/DeterminantOperators.h"
31 #include "Numerics/SymmetryOperations.h"
32 #include <sstream>
33
34 namespace qmcplusplus
35 {
36 using WP = WalkerProperties::Indexes;
37
WaveFunctionTester(MCWalkerConfiguration & w,TrialWaveFunction & psi,QMCHamiltonian & h,ParticleSetPool & ptclPool,Communicate * comm)38 WaveFunctionTester::WaveFunctionTester(MCWalkerConfiguration& w,
39 TrialWaveFunction& psi,
40 QMCHamiltonian& h,
41 ParticleSetPool& ptclPool,
42 Communicate* comm)
43 : QMCDriver(w, psi, h, comm, "WaveFunctionTester"),
44 PtclPool(ptclPool),
45 checkRatio("no"),
46 checkClone("no"),
47 checkHamPbyP("no"),
48 wftricks("no"),
49 checkEloc("no"),
50 checkBasic("yes"),
51 checkRatioV("no"),
52 deltaParam(0.0),
53 toleranceParam(0.0),
54 outputDeltaVsError(false),
55 checkSlaterDet(true)
56 {
57 m_param.add(checkRatio, "ratio");
58 m_param.add(checkClone, "clone");
59 m_param.add(checkHamPbyP, "hamiltonianpbyp");
60 m_param.add(sourceName, "source");
61 m_param.add(wftricks, "orbitalutility");
62 m_param.add(checkEloc, "printEloc");
63 m_param.add(checkBasic, "basic");
64 m_param.add(checkRatioV, "virtual_move");
65 m_param.add(deltaParam, "delta");
66 m_param.add(toleranceParam, "tolerance");
67 m_param.add(checkSlaterDetOption, "sd");
68
69 deltaR.resize(w.getTotalNum());
70 makeGaussRandom(deltaR);
71 }
72
~WaveFunctionTester()73 WaveFunctionTester::~WaveFunctionTester() {}
74
75 /*!
76 * \brief Test the evaluation of the wavefunction, gradient and laplacian
77 by comparing to the numerical evaluation.
78 *
79 Use the finite difference formulas formulas
80 \f[
81 \nabla_i f({\bf R}) = \frac{f({\bf R+\Delta r_i}) - f({\bf R})}{2\Delta r_i}
82 \f]
83 and
84 \f[
85 \nabla_i^2 f({\bf R}) = \sum_{x,y,z} \frac{f({\bf R}+\Delta x_i)
86 - 2 f({\bf R}) + f({\bf R}-\Delta x_i)}{2\Delta x_i^2},
87 \f]
88 where \f$ f = \ln \Psi \f$ and \f$ \Delta r_i \f$ is a
89 small displacement for the ith particle.
90 */
91
run()92 bool WaveFunctionTester::run()
93 {
94 //DistanceTable::create(1);
95 char fname[16];
96 sprintf(fname, "wftest.%03d", OHMMS::Controller->rank());
97 fout.open(fname);
98 fout.precision(15);
99
100 app_log() << "Starting a Wavefunction tester. Additional information in " << fname << std::endl;
101
102 put(qmcNode);
103
104 RandomGenerator_t* Rng1 = new RandomGenerator_t();
105 H.setRandomGenerator(Rng1);
106
107 if (checkSlaterDetOption == "no")
108 checkSlaterDet = false;
109 if (checkRatio == "yes")
110 {
111 //runRatioTest();
112 runRatioTest2();
113 }
114 else if (checkClone == "yes")
115 runCloneTest();
116 else if (checkEloc != "no")
117 printEloc();
118 else if (sourceName.size() != 0)
119 {
120 runGradSourceTest();
121 // runZeroVarianceTest();
122 }
123 else if (checkRatio == "deriv")
124 {
125 makeGaussRandom(deltaR);
126 deltaR *= 0.2;
127 runDerivTest();
128 runDerivNLPPTest();
129 }
130 else if (checkRatio == "derivclone")
131 runDerivCloneTest();
132 else if (wftricks == "rotate")
133 runwftricks();
134 else if (wftricks == "plot")
135 runNodePlot();
136 else if (checkBasic == "yes")
137 runBasicTest();
138 else if (checkRatioV == "yes")
139 runRatioV();
140 else
141 app_log() << "No wavefunction test specified" << std::endl;
142
143 //RealType ene = H.evaluate(W);
144 //app_log() << " Energy " << ene << std::endl;
145 return true;
146 }
147
runCloneTest()148 void WaveFunctionTester::runCloneTest()
149 {
150 app_log() << " ===== runCloneTest =====\n";
151 for (int iter = 0; iter < 4; ++iter)
152 {
153 app_log() << "Clone" << iter << std::endl;
154 ParticleSet* w_clone = new MCWalkerConfiguration(W);
155 TrialWaveFunction* psi_clone = Psi.makeClone(*w_clone);
156 QMCHamiltonian* h_clone = H.makeClone(*w_clone, *psi_clone);
157 h_clone->setPrimary(false);
158 int nat = W.getTotalNum();
159 MCWalkerConfiguration::PropertyContainer_t Properties(0,0,1,WP::MAXPROPERTIES);
160 //pick the first walker
161 MCWalkerConfiguration::Walker_t* awalker = *(W.begin());
162 //copy the properties of the working walker
163 Properties = awalker->Properties;
164 W.R = awalker->R;
165 W.update();
166 ValueType logpsi1 = Psi.evaluateLog(W);
167 RealType eloc1 = H.evaluate(W);
168 w_clone->R = awalker->R;
169 w_clone->update();
170 ValueType logpsi2 = psi_clone->evaluateLog(*w_clone);
171 RealType eloc2 = h_clone->evaluate(*w_clone);
172 app_log() << "Testing walker-by-walker functions " << std::endl;
173 app_log() << "log (original) = " << logpsi1 << " energy = " << eloc1 << std::endl;
174 app_log() << "log (clone) = " << logpsi2 << " energy = " << eloc2 << std::endl;
175 app_log() << "Testing pbyp functions " << std::endl;
176 Walker_t::WFBuffer_t& wbuffer(awalker->DataSet);
177 wbuffer.clear();
178 app_log() << " Walker Buffer State current=" << wbuffer.current() << " size=" << wbuffer.size() << std::endl;
179 Psi.registerData(W, wbuffer);
180 wbuffer.allocate();
181 Psi.copyFromBuffer(W, wbuffer);
182 Psi.evaluateLog(W);
183 logpsi1 = Psi.updateBuffer(W, wbuffer, false);
184 eloc1 = H.evaluate(W);
185 app_log() << " Walker Buffer State current=" << wbuffer.current() << " size=" << wbuffer.size() << std::endl;
186 wbuffer.clear();
187 app_log() << " Walker Buffer State current=" << wbuffer.current() << " size=" << wbuffer.size() << std::endl;
188 psi_clone->registerData(W, wbuffer);
189 wbuffer.allocate();
190 Psi.copyFromBuffer(W, wbuffer);
191 Psi.evaluateLog(W);
192 logpsi2 = Psi.updateBuffer(W, wbuffer, false);
193 eloc2 = H.evaluate(*w_clone);
194 app_log() << " Walker Buffer State current=" << wbuffer.current() << " size=" << wbuffer.size() << std::endl;
195 app_log() << "log (original) = " << logpsi1 << " energy = " << eloc1 << std::endl;
196 app_log() << "log (clone) = " << logpsi2 << " energy = " << eloc2 << std::endl;
197 delete h_clone;
198 delete psi_clone;
199 delete w_clone;
200 }
201 }
202
printEloc()203 void WaveFunctionTester::printEloc()
204 {
205 app_log() << " ===== printEloc =====\n";
206 ParticleSetPool::PoolType::iterator p;
207 for (p = PtclPool.getPool().begin(); p != PtclPool.getPool().end(); p++)
208 app_log() << "ParticelSet = " << p->first << std::endl;
209 // Find source ParticleSet
210 ParticleSetPool::PoolType::iterator pit(PtclPool.getPool().find(sourceName));
211 if (pit == PtclPool.getPool().end())
212 APP_ABORT("Unknown source \"" + sourceName + "\" in printEloc in WaveFunctionTester.");
213 ParticleSet& source = *((*pit).second);
214 app_log() << "Source = " << sourceName << " " << (*pit).first << std::endl;
215 int nel = W.getTotalNum();
216 int ncenter = source.getTotalNum();
217 // int ncenter = 3;
218 // std::cout <<"number of centers: " <<source.getTotalNum() << std::endl;
219 // std::cout <<"0: " <<source.R[0] << std::endl;
220 // std::cout <<"1: " <<source.R[1] << std::endl;
221 // std::cout <<"2: " <<source.R[2] << std::endl;
222 MCWalkerConfiguration::PropertyContainer_t Properties(0,0,1,WP::MAXPROPERTIES);;
223 //pick the first walker
224 MCWalkerConfiguration::Walker_t* awalker = *(W.begin());
225 //copy the properties of the working walker
226 Properties = awalker->Properties;
227 W.R = awalker->R;
228 W.update();
229 //ValueType psi = Psi.evaluate(W);
230 ValueType logpsi = Psi.evaluateLog(W);
231 RealType eloc = H.evaluate(W);
232 app_log() << " Logpsi: " << logpsi << std::endl;
233 app_log() << " HamTest "
234 << " Total " << eloc << std::endl;
235 for (int i = 0; i < H.sizeOfObservables(); i++)
236 app_log() << " HamTest " << H.getObservableName(i) << " " << H.getObservable(i) << std::endl;
237 //RealType psi = Psi.evaluateLog(W);
238 //int iat=0;
239 double maxR = 1000000.0;
240 std::vector<int> closestElectron(ncenter);
241 for (int iat = 0; iat < ncenter; iat++)
242 {
243 maxR = 10000000;
244 for (int k = 0; k < nel; k++)
245 {
246 double dx = std::sqrt((W.R[k][0] - source.R[iat][0]) * (W.R[k][0] - source.R[iat][0]) +
247 (W.R[k][1] - source.R[iat][1]) * (W.R[k][1] - source.R[iat][1]) +
248 (W.R[k][2] - source.R[iat][2]) * (W.R[k][2] - source.R[iat][2]));
249 if (dx < maxR)
250 {
251 maxR = dx;
252 closestElectron[iat] = k;
253 }
254 }
255 }
256 // closestElectron[iat]=1;
257 std::ofstream out("eloc.dat");
258 double x, dx = 1.0 / 499.0;
259 for (int k = 0; k < 500; k++)
260 {
261 x = -0.5 + k * dx;
262 out << x << " ";
263 for (int iat = 0; iat < ncenter; iat++)
264 {
265 PosType tempR = W.R[closestElectron[iat]];
266 W.R[closestElectron[iat]] = source.R[iat];
267 // W.R[closestElectron[iat]]=0.0;
268 W.R[closestElectron[iat]][0] += x;
269 W.update();
270 ValueType logpsi_p = Psi.evaluateLog(W);
271 ValueType ene = H.evaluate(W);
272 out << ene << " ";
273 W.R[closestElectron[iat]] = source.R[iat];
274 // W.R[closestElectron[iat]]=0.0;
275 W.R[closestElectron[iat]][1] += x;
276 W.update();
277 logpsi_p = Psi.evaluateLog(W);
278 ene = H.evaluate(W);
279 out << ene << " ";
280 W.R[closestElectron[iat]] = source.R[iat];
281 // W.R[closestElectron[iat]]=0.0;
282 W.R[closestElectron[iat]][2] += x;
283 W.update();
284 logpsi_p = Psi.evaluateLog(W);
285 ene = H.evaluate(W);
286 out << ene << " ";
287 W.R[closestElectron[iat]] = tempR;
288 }
289 out << std::endl;
290 }
291 out.close();
292 }
293
294
295 class FiniteDifference : public QMCTraits
296 {
297 public:
298 enum FiniteDiffType
299 {
300 FiniteDiff_LowOrder, // use simplest low-order formulas
301 FiniteDiff_Richardson // use Richardson extrapolation
302 };
FiniteDifference(FiniteDiffType fd_type=FiniteDiff_Richardson)303 FiniteDifference(FiniteDiffType fd_type = FiniteDiff_Richardson) : m_RichardsonSize(10), m_fd_type(fd_type) {}
304
305 int m_RichardsonSize;
306
307
308 FiniteDiffType m_fd_type;
309
310 struct PositionChange
311 {
312 int index; // particle index
313 PosType r;
314 };
315 typedef std::vector<PositionChange> PosChangeVector;
316 typedef std::vector<ValueType> ValueVector;
317
318
319 /** Generate points to evaluate */
320 void finiteDifferencePoints(RealType delta, MCWalkerConfiguration& W, PosChangeVector& positions);
321
322 /** Compute finite difference after log psi is computed for each point */
323 void computeFiniteDiff(RealType delta,
324 PosChangeVector& positions,
325 ValueVector& values,
326 ParticleSet::ParticleGradient_t& G_fd,
327 ParticleSet::ParticleLaplacian_t& L_fd);
328
329 void computeFiniteDiffLowOrder(RealType delta,
330 PosChangeVector& positions,
331 ValueVector& values,
332 ParticleSet::ParticleGradient_t& G_fd,
333 ParticleSet::ParticleLaplacian_t& L_fd);
334
335 void computeFiniteDiffRichardson(RealType delta,
336 PosChangeVector& positions,
337 ValueVector& values,
338 ParticleSet::ParticleGradient_t& G_fd,
339 ParticleSet::ParticleLaplacian_t& L_fd);
340 };
341
342
finiteDifferencePoints(RealType delta,MCWalkerConfiguration & W,PosChangeVector & positions)343 void FiniteDifference::finiteDifferencePoints(RealType delta, MCWalkerConfiguration& W, PosChangeVector& positions)
344 {
345 // First position is the central point
346 PositionChange p;
347 p.index = 0;
348 p.r = W.R[0];
349 positions.push_back(p);
350
351 int nat = W.getTotalNum();
352 for (int iat = 0; iat < nat; iat++)
353 {
354 PositionChange p;
355 p.index = iat;
356 PosType r0 = W.R[iat];
357
358 for (int idim = 0; idim < OHMMS_DIM; idim++)
359 {
360 p.r = r0;
361 p.r[idim] = r0[idim] - delta;
362 positions.push_back(p);
363
364 p.r = r0;
365 p.r[idim] = r0[idim] + delta;
366 positions.push_back(p);
367
368 if (m_fd_type == FiniteDiff_Richardson)
369 {
370 RealType dd = delta / 2;
371 for (int nr = 0; nr < m_RichardsonSize; nr++)
372 {
373 p.r = r0;
374 p.r[idim] = r0[idim] - dd;
375 positions.push_back(p);
376
377 p.r = r0;
378 p.r[idim] = r0[idim] + dd;
379 positions.push_back(p);
380
381 dd = dd / 2;
382 }
383 }
384 }
385 }
386 }
387
388
computeFiniteDiff(RealType delta,PosChangeVector & positions,ValueVector & values,ParticleSet::ParticleGradient_t & G_fd,ParticleSet::ParticleLaplacian_t & L_fd)389 void FiniteDifference::computeFiniteDiff(RealType delta,
390 PosChangeVector& positions,
391 ValueVector& values,
392 ParticleSet::ParticleGradient_t& G_fd,
393 ParticleSet::ParticleLaplacian_t& L_fd)
394 {
395 assert(positions.size() == values.size());
396 if (positions.size() == 0)
397 return;
398
399 ValueType logpsi = values[0];
400
401 if (m_fd_type == FiniteDiff_LowOrder)
402 {
403 computeFiniteDiffLowOrder(delta, positions, values, G_fd, L_fd);
404 }
405 else if (m_fd_type == FiniteDiff_Richardson)
406 {
407 computeFiniteDiffRichardson(delta, positions, values, G_fd, L_fd);
408 }
409 }
410
computeFiniteDiffLowOrder(RealType delta,PosChangeVector & positions,ValueVector & values,ParticleSet::ParticleGradient_t & G_fd,ParticleSet::ParticleLaplacian_t & L_fd)411 void FiniteDifference::computeFiniteDiffLowOrder(RealType delta,
412 PosChangeVector& positions,
413 ValueVector& values,
414 ParticleSet::ParticleGradient_t& G_fd,
415 ParticleSet::ParticleLaplacian_t& L_fd)
416 {
417 ValueType logpsi = values[0];
418
419 // lowest order derivative formula
420 ValueType c1 = 1.0 / delta / 2.0;
421 ValueType c2 = 1.0 / delta / delta;
422
423 const RealType twoD(2 * OHMMS_DIM);
424 const int pt_per_deriv = 2; // number of points per derivative
425 for (int pt_i = 1; pt_i < values.size(); pt_i += pt_per_deriv * OHMMS_DIM)
426 {
427 GradType g0;
428 ValueType lap0 = 0.0;
429 for (int idim = 0; idim < OHMMS_DIM; idim++)
430 {
431 int idx = pt_i + idim * pt_per_deriv;
432 ValueType logpsi_m = values[idx];
433 ValueType logpsi_p = values[idx + 1];
434
435 g0[idim] = logpsi_p - logpsi_m;
436 lap0 += logpsi_p + logpsi_m;
437 }
438
439 int iat = positions[pt_i].index;
440 GradType g = c1 * g0;
441 G_fd[iat] = g;
442
443 ValueType lap = c2 * (lap0 - twoD * logpsi);
444 //ValueType lap = c2*(lap0 - 2.0*OHMMS_DIM*logpsi);
445 L_fd[iat] = lap;
446 }
447 }
448
449 // Use Richardson extrapolation to compute the derivatives.
450 // The 'delta' parameter should not be small, as with fixed
451 // order finite difference methods. The algorithm will zoom
452 // in on the right size of parameter.
computeFiniteDiffRichardson(RealType delta,PosChangeVector & positions,ValueVector & values,ParticleSet::ParticleGradient_t & G_fd,ParticleSet::ParticleLaplacian_t & L_fd)453 void FiniteDifference::computeFiniteDiffRichardson(RealType delta,
454 PosChangeVector& positions,
455 ValueVector& values,
456 ParticleSet::ParticleGradient_t& G_fd,
457 ParticleSet::ParticleLaplacian_t& L_fd)
458 {
459 RealType tol = 1e-7;
460 ValueType logpsi = values[0];
461
462 const int pt_per_deriv = 2 * (m_RichardsonSize + 1); // number of points per derivative
463 for (int pt_i = 1; pt_i < values.size(); pt_i += pt_per_deriv * OHMMS_DIM)
464 {
465 GradType g0;
466 GradType gmin;
467 ValueType lmin;
468 std::vector<GradType> g_base(m_RichardsonSize + 1);
469 std::vector<GradType> g_rich(m_RichardsonSize + 1);
470 std::vector<GradType> g_prev(m_RichardsonSize + 1);
471
472 std::vector<ValueType> l_base(m_RichardsonSize + 1);
473 std::vector<ValueType> l_rich(m_RichardsonSize + 1);
474 std::vector<ValueType> l_prev(m_RichardsonSize + 1);
475
476 // Initial gradients and Laplacians at different deltas.
477 RealType dd = delta;
478 const ValueType ctwo(2);
479 for (int inr = 0; inr < m_RichardsonSize + 1; inr++)
480 {
481 RealType twodd = 2 * dd;
482 RealType ddsq = dd * dd;
483 l_base[inr] = 0.0;
484 for (int idim = 0; idim < OHMMS_DIM; idim++)
485 {
486 int idx = pt_i + idim * pt_per_deriv + 2 * inr;
487 ValueType logpsi_m = values[idx];
488 ValueType logpsi_p = values[idx + 1];
489 g_base[inr][idim] = (logpsi_p - logpsi_m) / twodd;
490 l_base[inr] += (logpsi_p + logpsi_m - ctwo * logpsi) / ddsq;
491 //g_base[inr][idim] = (logpsi_p - logpsi_m)/dd/2.0;
492 //l_base[inr] += (logpsi_p + logpsi_m - 2.0*logpsi)/(dd*dd);
493 }
494 dd = dd / 2;
495 }
496
497 // Gradient
498
499 g_prev[0] = g_base[0];
500 RealType fac = 1;
501 bool found = false;
502 for (int inr = 1; inr < m_RichardsonSize + 1; inr++)
503 {
504 g_rich[0] = g_base[inr];
505
506 fac *= 4;
507 for (int j = 1; j < inr + 1; j++)
508 {
509 g_rich[j] = g_rich[j - 1] + (g_rich[j - 1] - g_prev[j - 1]) / (fac - 1);
510 }
511
512 RealType err1 = 0.0;
513 RealType norm = 0.0;
514 for (int idim = 0; idim < OHMMS_DIM; idim++)
515 {
516 err1 += std::abs(g_rich[inr][idim] - g_prev[inr - 1][idim]);
517 norm += std::abs(g_prev[inr - 1][idim]);
518 }
519
520 RealType err_rel = err1 / norm;
521
522 // Not sure about the best stopping criteria
523 if (err_rel < tol)
524 {
525 gmin = g_rich[inr];
526 found = true;
527 break;
528 }
529 g_prev = g_rich;
530 }
531
532 if (!found)
533 {
534 gmin = g_rich[m_RichardsonSize];
535 }
536
537
538 // Laplacian
539 // TODO: eliminate the copied code between the gradient and Laplacian
540 // computations.
541
542 l_prev[0] = l_base[0];
543
544 fac = 1;
545 found = false;
546 for (int inr = 1; inr < m_RichardsonSize + 1; inr++)
547 {
548 l_rich[0] = l_base[inr];
549
550 fac *= 4;
551 for (int j = 1; j < inr + 1; j++)
552 {
553 l_rich[j] = l_rich[j - 1] + (l_rich[j - 1] - l_prev[j - 1]) / (fac - 1);
554 }
555
556 RealType err1 = std::abs(l_rich[inr] - l_prev[inr - 1]);
557 RealType err_rel = std::abs(err1 / l_prev[inr - 1]);
558
559 if (err_rel < tol)
560 {
561 lmin = l_rich[inr];
562 found = true;
563 break;
564 }
565 l_prev = l_rich;
566 }
567
568 if (!found)
569 {
570 lmin = l_rich[m_RichardsonSize];
571 }
572
573 int iat = positions[pt_i].index;
574 G_fd[iat] = gmin;
575 L_fd[iat] = lmin;
576 }
577 }
578
579 // Compute numerical gradient and Laplacian
computeNumericalGrad(RealType delta,ParticleSet::ParticleGradient_t & G_fd,ParticleSet::ParticleLaplacian_t & L_fd)580 void WaveFunctionTester::computeNumericalGrad(RealType delta,
581 ParticleSet::ParticleGradient_t& G_fd, // finite difference
582 ParticleSet::ParticleLaplacian_t& L_fd)
583 {
584 FiniteDifference fd(FiniteDifference::FiniteDiff_LowOrder);
585 //FiniteDifference fd(FiniteDifference::FiniteDiff_Richardson);
586 FiniteDifference::PosChangeVector positions;
587
588 fd.finiteDifferencePoints(delta, W, positions);
589
590 FiniteDifference::ValueVector logpsi_vals;
591 FiniteDifference::PosChangeVector::iterator it;
592
593 for (it = positions.begin(); it != positions.end(); it++)
594 {
595 PosType r0 = W.R[it->index];
596 W.R[it->index] = it->r;
597 W.update();
598 RealType logpsi0 = Psi.evaluateLog(W);
599 #if defined(QMC_COMPLEX)
600 RealType phase0 = Psi.getPhase();
601 ValueType logpsi = std::complex<OHMMS_PRECISION>(logpsi0, phase0);
602 #else
603 ValueType logpsi = logpsi0;
604 #endif
605 logpsi_vals.push_back(logpsi);
606
607 W.R[it->index] = r0;
608 W.update();
609 Psi.evaluateLog(W);
610 }
611
612 fd.computeFiniteDiff(delta, positions, logpsi_vals, G_fd, L_fd);
613 }
614
615 // Usually
616 // lower_iat = 0
617 // upper_iat = nat
checkGradients(int lower_iat,int upper_iat,ParticleSet::ParticleGradient_t & G,ParticleSet::ParticleLaplacian_t & L,ParticleSet::ParticleGradient_t & G_fd,ParticleSet::ParticleLaplacian_t & L_fd,std::stringstream & log,int indent)618 bool WaveFunctionTester::checkGradients(int lower_iat,
619 int upper_iat,
620 ParticleSet::ParticleGradient_t& G,
621 ParticleSet::ParticleLaplacian_t& L,
622 ParticleSet::ParticleGradient_t& G_fd,
623 ParticleSet::ParticleLaplacian_t& L_fd,
624 std::stringstream& log,
625 int indent /* = 0 */)
626 {
627 ParticleSet::Scalar_t rel_tol = 1e-3;
628 ParticleSet::Scalar_t abs_tol = 1e-7;
629 if (toleranceParam > 0.0)
630 {
631 rel_tol = toleranceParam;
632 }
633
634 bool all_okay = true;
635 std::string pad(4 * indent, ' ');
636
637
638 for (int iat = lower_iat; iat < upper_iat; iat++)
639 {
640 ParticleSet::Scalar_t L_err = std::abs(L[iat] - L_fd[iat]);
641 ParticleSet::Scalar_t L_rel_denom = std::max(std::abs(L[iat]), std::abs(L_fd[iat]));
642 ParticleSet::Scalar_t L_err_rel = std::abs(L_err / L_rel_denom);
643
644 if (L_err_rel > rel_tol && L_err > abs_tol)
645 {
646 if (L_err_rel > rel_tol)
647 {
648 log << pad << "Finite difference Laplacian exceeds relative tolerance (" << rel_tol << ") for particle " << iat
649 << std::endl;
650 }
651 else
652 {
653 log << pad << "Finite difference Laplacian exceeds absolute tolerance (" << abs_tol << ") for particle " << iat
654 << std::endl;
655 }
656 log << pad << " Analytic = " << L[iat] << std::endl;
657 log << pad << " Finite diff = " << L_fd[iat] << std::endl;
658 log << pad << " Error = " << L_err << " Relative Error = " << L_err_rel << std::endl;
659 all_okay = false;
660 }
661
662 ParticleSet::Scalar_t G_err_rel[OHMMS_DIM];
663 for (int idim = 0; idim < OHMMS_DIM; idim++)
664 {
665 ParticleSet::Scalar_t G_err = std::abs(G[iat][idim] - G_fd[iat][idim]);
666 ParticleSet::Scalar_t G_rel_denom = std::max(std::abs(G[iat][idim]), std::abs(G_fd[iat][idim]));
667 G_err_rel[idim] = std::abs(G_err / G[iat][idim]);
668
669 if (G_err_rel[idim] > rel_tol && G_err > abs_tol)
670 {
671 if (G_err_rel[idim] > rel_tol)
672 {
673 log << pad << "Finite difference gradient exceeds relative tolerance (" << rel_tol << ") for particle "
674 << iat;
675 }
676 else
677 {
678 log << pad << "Finite difference gradient exceeds absolute tolerance (" << abs_tol << ") for particle "
679 << iat;
680 }
681 log << " component " << idim << std::endl;
682 log << pad << " Analytic = " << G[iat][idim] << std::endl;
683 log << pad << " Finite diff = " << G_fd[iat][idim] << std::endl;
684 log << pad << " Error = " << G_err << " Relative Error = " << G_err_rel[idim] << std::endl;
685 all_okay = false;
686 }
687 }
688
689 fout << pad << "For particle #" << iat << " at " << W.R[iat] << std::endl;
690 fout << pad << "Gradient = " << std::setw(12) << G[iat] << std::endl;
691 fout << pad << " Finite diff = " << std::setw(12) << G_fd[iat] << std::endl;
692 fout << pad << " Error = " << std::setw(12) << G[iat] - G_fd[iat] << std::endl;
693 fout << pad << " Relative Error = ";
694 for (int idim = 0; idim < OHMMS_DIM; idim++)
695 {
696 fout << G_err_rel[idim] << " ";
697 }
698 fout << std::endl << std::endl;
699 fout << pad << "Laplacian = " << std::setw(12) << L[iat] << std::endl;
700 fout << pad << " Finite diff = " << std::setw(12) << L_fd[iat] << std::endl;
701 fout << pad << " Error = " << std::setw(12) << L[iat] - L_fd[iat] << " Relative Error = " << L_err_rel
702 << std::endl
703 << std::endl;
704 }
705 return all_okay;
706 }
707
checkGradientAtConfiguration(MCWalkerConfiguration::Walker_t * W1,std::stringstream & fail_log,bool & ignore)708 bool WaveFunctionTester::checkGradientAtConfiguration(MCWalkerConfiguration::Walker_t* W1,
709 std::stringstream& fail_log,
710 bool& ignore)
711 {
712 int nat = W.getTotalNum();
713 ParticleSet::ParticleGradient_t G(nat), G1(nat);
714 ParticleSet::ParticleLaplacian_t L(nat), L1(nat);
715
716 W.loadWalker(*W1, true);
717
718 // compute analytic values
719 Psi.evaluateLog(W);
720 G = W.G;
721 L = W.L;
722
723 // Use a single delta with a fairly large tolerance
724 //computeNumericalGrad(delta, G1, L1);
725
726 RealType delta = 1.0e-4;
727 if (deltaParam > 0.0)
728 {
729 delta = deltaParam;
730 }
731 FiniteDifference fd(FiniteDifference::FiniteDiff_LowOrder);
732 //RealType delta = 1.0;
733 //FiniteDifference fd(FiniteDifference::FiniteDiff_Richardson);
734
735 FiniteDifference::PosChangeVector positions;
736
737 fd.finiteDifferencePoints(delta, W, positions);
738
739 FiniteDifference::ValueVector logpsi_vals;
740 FiniteDifference::PosChangeVector::iterator it;
741
742 for (it = positions.begin(); it != positions.end(); it++)
743 {
744 PosType r0 = W.R[it->index];
745 W.R[it->index] = it->r;
746 W.update();
747 RealType logpsi0 = Psi.evaluateLog(W);
748 RealType phase0 = Psi.getPhase();
749 #if defined(QMC_COMPLEX)
750 ValueType logpsi = std::complex<OHMMS_PRECISION>(logpsi0, phase0);
751 #else
752 ValueType logpsi = logpsi0;
753 #endif
754 logpsi_vals.push_back(logpsi);
755
756 W.R[it->index] = r0;
757 W.update();
758 Psi.evaluateLog(W);
759 }
760
761 fd.computeFiniteDiff(delta, positions, logpsi_vals, G1, L1);
762
763 fout << "delta = " << delta << std::endl;
764
765 // TODO - better choice of tolerance
766 // TODO - adjust delta and tolerance based on precision of wavefunction
767
768 bool all_okay = checkGradients(0, nat, G, L, G1, L1, fail_log);
769 // RealType tol = 1e-3;
770 // if (toleranceParam> 0.0)
771 // {
772 // tol = toleranceParam;
773 // }
774
775 for (int iorb = 0; iorb < Psi.getOrbitals().size(); iorb++)
776 {
777 WaveFunctionComponent* orb = Psi.getOrbitals()[iorb];
778
779 ParticleSet::ParticleGradient_t G(nat), tmpG(nat), G1(nat);
780 ParticleSet::ParticleLaplacian_t L(nat), tmpL(nat), L1(nat);
781
782
783 LogValueType logpsi1 = orb->evaluateLog(W, G, L);
784
785 fail_log << "WaveFunctionComponent " << iorb << " " << orb->ClassName << " log psi = " << logpsi1 << std::endl;
786
787 FiniteDifference::ValueVector logpsi_vals;
788 FiniteDifference::PosChangeVector::iterator it;
789 for (it = positions.begin(); it != positions.end(); it++)
790 {
791 PosType r0 = W.R[it->index];
792 W.R[it->index] = it->r;
793 W.update();
794 ParticleSet::SingleParticlePos_t zeroR;
795 W.makeMove(it->index, zeroR);
796
797 LogValueType logpsi0 = orb->evaluateLog(W, tmpG, tmpL);
798 #if defined(QMC_COMPLEX)
799 ValueType logpsi(logpsi0.real(), logpsi0.imag());
800 #else
801 ValueType logpsi = std::real(logpsi0);
802 #endif
803 logpsi_vals.push_back(logpsi);
804 W.rejectMove(it->index);
805
806 W.R[it->index] = r0;
807 W.update();
808 }
809 fd.computeFiniteDiff(delta, positions, logpsi_vals, G1, L1);
810
811 fout << " WaveFunctionComponent " << iorb << " " << orb->ClassName << std::endl;
812
813 if (!checkGradients(0, nat, G, L, G1, L1, fail_log, 1))
814 {
815 all_okay = false;
816 }
817
818 if (!checkSlaterDet)
819 continue; // skip SlaterDet check if <backflow> is present
820 // DiracDeterminantWithBackflow::evaluateLog requires a call to BackflowTransformation::evaluate in its owning SlaterDetWithBackflow to work correctly.
821 SlaterDet* sd = dynamic_cast<SlaterDet*>(orb);
822 if (sd)
823 {
824 for (int isd = 0; isd < sd->Dets.size(); isd++)
825 {
826 ParticleSet::ParticleGradient_t G(nat), tmpG(nat), G1(nat);
827 ParticleSet::ParticleLaplacian_t L(nat), tmpL(nat), L1(nat);
828 DiracDeterminantBase& det = *sd->Dets[isd];
829 LogValueType logpsi2 = det.evaluateLog(W, G, L); // this won't work with backflow
830 fail_log << " Slater Determiant " << isd << " (for particles " << det.getFirstIndex() << " to "
831 << det.getLastIndex() << ") log psi = " << logpsi2 << std::endl;
832 // Should really check the condition number on the matrix determinant.
833 // For now, just ignore values that too small.
834 if (std::real(logpsi2) < -40.0)
835 {
836 ignore = true;
837 }
838 FiniteDifference::ValueVector logpsi_vals;
839 FiniteDifference::PosChangeVector::iterator it;
840 for (it = positions.begin(); it != positions.end(); it++)
841 {
842 PosType r0 = W.R[it->index];
843 W.R[it->index] = it->r;
844 W.update();
845
846 LogValueType logpsi0 = det.evaluateLog(W, tmpG, tmpL);
847 #if defined(QMC_COMPLEX)
848 ValueType logpsi(logpsi0.real(), logpsi0.imag());
849 #else
850 ValueType logpsi = std::real(logpsi0);
851 #endif
852 logpsi_vals.push_back(logpsi);
853
854 W.R[it->index] = r0;
855 W.update();
856 }
857 fd.computeFiniteDiff(delta, positions, logpsi_vals, G1, L1);
858
859 if (!checkGradients(det.getFirstIndex(), det.getLastIndex(), G, L, G1, L1, fail_log, 2))
860 {
861 all_okay = false;
862 }
863
864 #if 0
865 // Testing single particle orbitals doesn't work yet - probably something
866 // with setup after setting the position.
867 std::map<std::string, SPOSetPtr>::iterator spo_it = sd->mySPOSet.begin();
868 for (; spo_it != sd->mySPOSet.end(); spo_it++)
869 {
870 SPOSetPtr spo = spo_it->second;
871 fail_log << " SPO set = " << spo_it->first << " name = " << spo->className;
872 fail_log << " orbital set size = " << spo->size();
873 fail_log << " basis set size = " << spo->getBasisSetSize() << std::endl;
874
875 ParticleSet::ParticleGradient_t G(nat), tmpG(nat), G1(nat);
876 ParticleSet::ParticleLaplacian_t L(nat), tmpL(nat), L1(nat);
877 RealType logpsi3 = det.evaluateLog(W, G, L);
878 FiniteDifference::ValueVector logpsi_vals;
879 FiniteDifference::PosChangeVector::iterator it;
880 for (it = positions.begin(); it != positions.end(); it++)
881 {
882 PosType r0 = W.R[it->index];
883 W.R[it->index] = it->r;
884 W.update();
885 ParticleSet::SingleParticlePos_t zeroR;
886 W.makeMove(it->index,zeroR);
887
888 SPOSet::ValueVector_t psi(spo->size());
889
890 spo->evaluate(W, it->index, psi);
891 ValueType logpsi = psi[0];
892 logpsi_vals.push_back(logpsi);
893
894 W.rejectMove(it->index);
895 W.R[it->index] = r0;
896 W.update();
897 }
898 fd.computeFiniteDiff(delta, positions, logpsi_vals, G1, L1);
899
900 if (!checkGradients(det.getFirstIndex(), det.getLastIndex(), G, L, G1, L1, fail_log, 3))
901 {
902 all_okay = false;
903 }
904 }
905 #endif
906 }
907 }
908 }
909 return all_okay;
910 }
911
runBasicTest()912 void WaveFunctionTester::runBasicTest()
913 {
914 int nat = W.getTotalNum();
915 fout << "Numerical gradient and Laplacian test" << std::endl;
916
917 std::stringstream fail_log;
918 bool all_okay = true;
919 int fails = 0;
920 int nconfig = 0;
921 int nignore = 0;
922
923 MCWalkerConfiguration::iterator Wit(W.begin());
924 for (; Wit != W.end(); Wit++)
925 {
926 fout << "Walker # " << nconfig << std::endl;
927 std::stringstream fail_log1;
928 bool ignore = false;
929 bool this_okay = checkGradientAtConfiguration(*Wit, fail_log1, ignore);
930 if (ignore)
931 {
932 nignore++;
933 }
934 if (!this_okay && !ignore)
935 {
936 fail_log << "Walker # " << nconfig << std::endl;
937 fail_log << fail_log1.str();
938 fail_log << std::endl;
939 fails++;
940 all_okay = false;
941 }
942 nconfig++;
943 }
944
945 app_log() << "Number of samples = " << nconfig << std::endl;
946 app_log() << "Number ignored (bad positions) = " << nignore << std::endl << std::endl;
947 app_log() << "Number of fails = " << fails << std::endl << std::endl;
948
949 if (!all_okay)
950 {
951 std::string fail_name("wf_fails.dat");
952 app_log() << "More detail on finite difference failures in " << fail_name << std::endl;
953 std::ofstream eout(fail_name.c_str());
954 eout << fail_log.str();
955 eout.close();
956 }
957
958 app_log() << "Finite difference test: " << (all_okay ? "PASS" : "FAIL") << std::endl;
959
960 // Compute approximation error vs. delta.
961 if (outputDeltaVsError)
962 {
963 double delta = 1.0;
964 bool inputOkay = true;
965
966 std::ofstream dout(DeltaVsError.outputFile.c_str());
967
968 int iat = DeltaVsError.particleIndex;
969 int ig = DeltaVsError.gradientComponentIndex;
970
971 if (iat < 0 || iat >= nat)
972 {
973 dout << "# Particle index (" << iat << ") is out of range (0 - " << nat - 1 << ")" << std::endl;
974 inputOkay = false;
975 }
976
977 if (ig < 0 || ig >= OHMMS_DIM)
978 {
979 dout << "# Gradient component index (" << ig << ") is out of range (0 - " << OHMMS_DIM - 1 << ")" << std::endl;
980 inputOkay = false;
981 }
982
983 if (inputOkay)
984 {
985 dout << "# Particle = " << iat << " Gradient component = " << ig << std::endl;
986 dout << "#" << std::setw(11) << "delta" << std::setw(14) << "G_err_rel" << std::setw(14) << "L_err_rel"
987 << std::endl;
988 ParticleSet::ParticleGradient_t G(nat), G1(nat);
989 ParticleSet::ParticleLaplacian_t L(nat), L1(nat);
990 for (int i = 0; i < 20; i++)
991 {
992 // compute analytic values
993 G = W.G;
994 L = W.L;
995 Psi.evaluateLog(W);
996
997 computeNumericalGrad(delta, G1, L1);
998 ParticleSet::Scalar_t L_err = std::abs(L[iat] - L1[iat]);
999 ParticleSet::Scalar_t L_err_rel = std::abs(L_err / L[iat]);
1000 ParticleSet::Scalar_t G_err = std::abs(G[iat][ig] - G1[iat][ig]);
1001 ParticleSet::Scalar_t G_err_rel = std::abs(G_err / G[iat][ig]);
1002 dout << std::setw(12) << delta;
1003 dout << std::setw(14) << std::abs(G_err_rel) << std::setw(14) << std::abs(L_err_rel);
1004 dout << std::endl;
1005 delta *= std::sqrt(0.1);
1006 }
1007 }
1008 dout.close();
1009 }
1010
1011 fout << "Ratio test" << std::endl;
1012
1013 RealType tol = 1e-3;
1014 RealType ratio_tol = 1e-9;
1015 bool any_ratio_fail = false;
1016 makeGaussRandom(deltaR);
1017 fout << "deltaR:" << std::endl;
1018 fout << deltaR << std::endl;
1019 fout << "Particle Ratio of Ratios Computed Ratio Internal Ratio" << std::endl;
1020 //testing ratio alone
1021 for (int iat = 0; iat < nat; iat++)
1022 {
1023 W.update();
1024 //ValueType psi_p = log(std::abs(Psi.evaluate(W)));
1025 RealType psi_p = Psi.evaluateLog(W);
1026 RealType phase_p = Psi.getPhase();
1027 W.makeMove(iat, deltaR[iat]);
1028 //W.update();
1029 ValueType aratio = Psi.calcRatio(W, iat);
1030 RealType phaseDiff = Psi.getPhaseDiff();
1031 W.rejectMove(iat);
1032 Psi.rejectMove(iat);
1033 W.R[iat] += deltaR[iat];
1034 W.update();
1035 //ValueType psi_m = log(std::abs(Psi.evaluate(W)));
1036 RealType psi_m = Psi.evaluateLog(W);
1037 RealType phase_m = Psi.getPhase();
1038
1039 #if defined(QMC_COMPLEX)
1040 RealType ratioMag = std::exp(psi_m - psi_p);
1041 RealType dphase = phase_m - phase_p;
1042 if (phaseDiff < 0.0)
1043 phaseDiff += 2.0 * M_PI;
1044 if (phaseDiff > 2.0 * M_PI)
1045 phaseDiff -= 2.0 * M_PI;
1046 if (dphase < 0.0)
1047 dphase += 2.0 * M_PI;
1048 if (dphase > 2.0 * M_PI)
1049 dphase -= 2.0 * M_PI;
1050 ValueType ratDiff = std::complex<OHMMS_PRECISION>(ratioMag * std::cos(dphase), ratioMag * std::sin(dphase));
1051 // TODO - test complex ratio against a tolerance
1052 fout << iat << " " << aratio / ratDiff << " " << ratDiff << " " << aratio << std::endl;
1053 fout << " ratioMag " << std::abs(aratio) / ratioMag << " " << ratioMag << std::endl;
1054 fout << " PhaseDiff " << phaseDiff / dphase << " " << phaseDiff << " " << dphase << std::endl;
1055 #else
1056 RealType ratDiff = std::exp(psi_m - psi_p) * std::cos(phase_m - phase_p);
1057 fout << iat << " " << aratio / ratDiff << " " << ratDiff << " " << aratio << std::endl;
1058 if (std::abs(aratio / ratDiff - 1.0) > ratio_tol)
1059 {
1060 app_log() << "Wavefunction ratio exceeds tolerance " << tol << ") for particle " << iat << std::endl;
1061 app_log() << " Internally computed ratio = " << aratio << std::endl;
1062 app_log() << " Separately computed ratio = " << ratDiff << std::endl;
1063 any_ratio_fail = true;
1064 }
1065 #endif
1066 }
1067 app_log() << "Ratio test: " << (any_ratio_fail ? "FAIL" : "PASS") << std::endl;
1068 }
1069
runRatioTest()1070 void WaveFunctionTester::runRatioTest()
1071 {
1072 #if 0
1073 int nat = W.getTotalNum();
1074 ParticleSet::ParticleGradient_t Gp(nat), dGp(nat);
1075 ParticleSet::ParticleLaplacian_t Lp(nat), dLp(nat);
1076 bool checkHam=(checkHamPbyP == "yes");
1077 Tau=0.025;
1078 MCWalkerConfiguration::iterator it(W.begin()), it_end(W.end());
1079 while (it != it_end)
1080 {
1081 makeGaussRandom(deltaR);
1082 Walker_t::WFBuffer_t tbuffer;
1083 W.R = (**it).R+Tau*deltaR;
1084 (**it).R=W.R;
1085 W.update();
1086 RealType logpsi=Psi.registerData(W,tbuffer);
1087 RealType ene;
1088 if (checkHam)
1089 ene = H.registerData(W,tbuffer);
1090 else
1091 ene = H.evaluate(W);
1092 (*it)->DataSet=tbuffer;
1093 //RealType ene = H.evaluate(W);
1094 (*it)->resetProperty(logpsi,Psi.getPhase(),ene,0.0,0.0,1.0);
1095 H.saveProperty((*it)->getPropertyBase());
1096 ++it;
1097 app_log() << " HamTest " << " Total " << ene << std::endl;
1098 for (int i=0; i<H.sizeOfObservables(); i++)
1099 app_log() << " HamTest " << H.getObservableName(i) << " " << H.getObservable(i) << std::endl;
1100 }
1101 fout << " Update using drift " << std::endl;
1102 bool pbyp_mode=true;
1103 for (int iter=0; iter<4; ++iter)
1104 {
1105 int iw=0;
1106 it=W.begin();
1107 while (it != it_end)
1108 {
1109 fout << "\nStart Walker " << iw++ << std::endl;
1110 Walker_t& thisWalker(**it);
1111 W.loadWalker(thisWalker,pbyp_mode);
1112 Walker_t::WFBuffer_t& w_buffer(thisWalker.DataSet);
1113 Psi.copyFromBuffer(W,w_buffer);
1114 H.copyFromBuffer(W,w_buffer);
1115 // Psi.evaluateLog(W);
1116 RealType eold(thisWalker.Properties(WP::LOCALENERGY));
1117 RealType logpsi(thisWalker.Properties(LOGPSI));
1118 RealType emixed(eold), enew(eold);
1119 makeGaussRandom(deltaR);
1120 //mave a move
1121 RealType ratio_accum(1.0);
1122 for (int iat=0; iat<nat; iat++)
1123 {
1124 PosType dr(Tau*deltaR[iat]);
1125 PosType newpos(W.makeMove(iat,dr));
1126 //RealType ratio=Psi.ratio(W,iat,dGp,dLp);
1127 W.dG=0;
1128 W.dL=0;
1129 RealType ratio=Psi.ratio(W,iat,W.dG,W.dL);
1130 Gp = W.G + W.dG;
1131 //RealType enew = H.evaluatePbyP(W,iat);
1132 if (checkHam)
1133 enew = H.evaluatePbyP(W,iat);
1134 if (ratio > Random())
1135 {
1136 fout << " Accepting a move for " << iat << std::endl;
1137 fout << " Energy after a move " << enew << std::endl;
1138 W.G += W.dG;
1139 W.L += W.dL;
1140 W.acceptMove(iat);
1141 Psi.acceptMove(W,iat);
1142 if (checkHam)
1143 H.acceptMove(iat);
1144 ratio_accum *= ratio;
1145 }
1146 else
1147 {
1148 fout << " Rejecting a move for " << iat << std::endl;
1149 W.rejectMove(iat);
1150 Psi.rejectMove(iat);
1151 //H.rejectMove(iat);
1152 }
1153 }
1154 fout << " Energy after pbyp = " << H.getLocalEnergy() << std::endl;
1155 RealType newlogpsi_up = Psi.evaluateLog(W,w_buffer);
1156 W.saveWalker(thisWalker);
1157 RealType ene_up;
1158 if (checkHam)
1159 ene_up= H.evaluate(W,w_buffer);
1160 else
1161 ene_up = H.evaluate(W);
1162 Gp=W.G;
1163 Lp=W.L;
1164 W.R=thisWalker.R;
1165 W.update();
1166 RealType newlogpsi=Psi.updateBuffer(W,w_buffer,false);
1167 RealType ene = H.evaluate(W);
1168 thisWalker.resetProperty(newlogpsi,Psi.getPhase(),ene);
1169 //thisWalker.resetProperty(std::log(psi),Psi.getPhase(),ene);
1170 fout << iter << " Energy by update = "<< ene_up << " " << ene << " " << ene_up-ene << std::endl;
1171 fout << iter << " Ratio " << ratio_accum*ratio_accum
1172 << " | " << std::exp(2.0*(newlogpsi-logpsi)) << " "
1173 << ratio_accum*ratio_accum/std::exp(2.0*(newlogpsi-logpsi)) << std::endl
1174 << " new log(psi) updated " << newlogpsi_up
1175 << " new log(psi) calculated " << newlogpsi
1176 << " old log(psi) " << logpsi << std::endl;
1177 fout << " Gradients " << std::endl;
1178 for (int iat=0; iat<nat; iat++)
1179 fout << W.G[iat]-Gp[iat] << W.G[iat] << std::endl; //W.G[iat] << G[iat] << std::endl;
1180 fout << " Laplacians " << std::endl;
1181 for (int iat=0; iat<nat; iat++)
1182 fout << W.L[iat]-Lp[iat] << " " << W.L[iat] << std::endl;
1183 ++it;
1184 }
1185 }
1186 fout << " Update without drift : for VMC useDrift=\"no\"" << std::endl;
1187 for (int iter=0; iter<4; ++iter)
1188 {
1189 it=W.begin();
1190 int iw=0;
1191 while (it != it_end)
1192 {
1193 fout << "\nStart Walker " << iw++ << std::endl;
1194 Walker_t& thisWalker(**it);
1195 W.loadWalker(thisWalker,pbyp_mode);
1196 Walker_t::WFBuffer_t& w_buffer(thisWalker.DataSet);
1197 //Psi.updateBuffer(W,w_buffer,true);
1198 Psi.copyFromBuffer(W,w_buffer);
1199 RealType eold(thisWalker.Properties(WP::LOCALENERGY));
1200 RealType logpsi(thisWalker.Properties(LOGPSI));
1201 RealType emixed(eold), enew(eold);
1202 //mave a move
1203 RealType ratio_accum(1.0);
1204 for (int substep=0; substep<3; ++substep)
1205 {
1206 makeGaussRandom(deltaR);
1207 for (int iat=0; iat<nat; iat++)
1208 {
1209 PosType dr(Tau*deltaR[iat]);
1210 PosType newpos(W.makeMove(iat,dr));
1211 RealType ratio=Psi.ratio(W,iat);
1212 RealType prob = ratio*ratio;
1213 if (prob > Random())
1214 {
1215 fout << " Accepting a move for " << iat << std::endl;
1216 W.acceptMove(iat);
1217 Psi.acceptMove(W,iat);
1218 ratio_accum *= ratio;
1219 }
1220 else
1221 {
1222 fout << " Rejecting a move for " << iat << std::endl;
1223 W.rejectMove(iat);
1224 Psi.rejectMove(iat);
1225 }
1226 }
1227 RealType logpsi_up = Psi.updateBuffer(W,w_buffer,false);
1228 W.saveWalker(thisWalker);
1229 RealType ene = H.evaluate(W);
1230 thisWalker.resetProperty(logpsi_up,Psi.getPhase(),ene);
1231 }
1232 Gp=W.G;
1233 Lp=W.L;
1234 W.update();
1235 RealType newlogpsi=Psi.evaluateLog(W);
1236 fout << iter << " Ratio " << ratio_accum*ratio_accum
1237 << " | " << std::exp(2.0*(newlogpsi-logpsi)) << " "
1238 << ratio_accum*ratio_accum/std::exp(2.0*(newlogpsi-logpsi)) << std::endl
1239 << " new log(psi) " << newlogpsi
1240 << " old log(psi) " << logpsi << std::endl;
1241 fout << " Gradients " << std::endl;
1242 for (int iat=0; iat<nat; iat++)
1243 {
1244 fout << W.G[iat]-Gp[iat] << W.G[iat] << std::endl; //W.G[iat] << G[iat] << std::endl;
1245 }
1246 fout << " Laplacians " << std::endl;
1247 for (int iat=0; iat<nat; iat++)
1248 {
1249 fout << W.L[iat]-Lp[iat] << " " << W.L[iat] << std::endl;
1250 }
1251 ++it;
1252 }
1253 }
1254 //for(it=W.begin();it != it_end; ++it)
1255 //{
1256 // Walker_t& thisWalker(**it);
1257 // Walker_t::WFBuffer_t& w_buffer((*it)->DataSet);
1258 // w_buffer.rewind();
1259 // W.updateBuffer(**it,w_buffer);
1260 // RealType logpsi=Psi.updateBuffer(W,w_buffer,true);
1261 //}
1262 #endif
1263 }
1264
runRatioTest2()1265 void WaveFunctionTester::runRatioTest2()
1266 {
1267 app_log() << " ===== runRatioTest2 =====\n";
1268 int nat = W.getTotalNum();
1269 ParticleSet::ParticleGradient_t Gp(nat), dGp(nat);
1270 ParticleSet::ParticleLaplacian_t Lp(nat), dLp(nat);
1271 Tau = 0.025;
1272 MCWalkerConfiguration::iterator it(W.begin()), it_end(W.end());
1273 for (; it != it_end; ++it)
1274 {
1275 makeGaussRandom(deltaR);
1276 Walker_t::WFBuffer_t tbuffer;
1277 (**it).R += Tau * deltaR;
1278 W.loadWalker(**it, true);
1279 Psi.registerData(W, tbuffer);
1280 tbuffer.allocate();
1281 Psi.copyFromBuffer(W, tbuffer);
1282 Psi.evaluateLog(W);
1283 RealType logpsi = Psi.updateBuffer(W, tbuffer, false);
1284 RealType ene = H.evaluate(W);
1285 (*it)->DataSet = tbuffer;
1286 //RealType ene = H.evaluate(W);
1287 (*it)->resetProperty(logpsi, Psi.getPhase(), ene, 0.0, 0.0, 1.0);
1288 H.saveProperty((*it)->getPropertyBase());
1289 app_log() << " HamTest "
1290 << " Total " << ene << std::endl;
1291 for (int i = 0; i < H.sizeOfObservables(); i++)
1292 app_log() << " HamTest " << H.getObservableName(i) << " " << H.getObservable(i) << std::endl;
1293 }
1294 for (int iter = 0; iter < 20; ++iter)
1295 {
1296 int iw = 0;
1297 it = W.begin();
1298 //while(it != it_end)
1299 for (; it != it_end; ++it)
1300 {
1301 fout << "\nStart Walker " << iw++ << std::endl;
1302 Walker_t& thisWalker(**it);
1303 W.loadWalker(thisWalker, true);
1304 Walker_t::WFBuffer_t& w_buffer(thisWalker.DataSet);
1305 Psi.copyFromBuffer(W, w_buffer);
1306 RealType eold(thisWalker.Properties(WP::LOCALENERGY));
1307 RealType logpsi(thisWalker.Properties(WP::LOGPSI));
1308 Psi.evaluateLog(W);
1309 ParticleSet::ParticleGradient_t realGrad(W.G);
1310 makeGaussRandom(deltaR);
1311 //mave a move
1312 for (int iat = 0; iat < nat; iat++)
1313 {
1314 TinyVector<ParticleSet::SingleParticleValue_t, OHMMS_DIM> grad_now = Psi.evalGrad(W, iat);
1315 GradType grad_new;
1316 for (int sds = 0; sds < 3; sds++)
1317 fout << realGrad[iat][sds] - grad_now[sds] << " ";
1318 PosType dr(Tau * deltaR[iat]);
1319 W.makeMove(iat, dr);
1320 ValueType ratio2 = Psi.calcRatioGrad(W, iat, grad_new);
1321 W.rejectMove(iat);
1322 Psi.rejectMove(iat);
1323 W.makeMove(iat, dr);
1324 ValueType ratio1 = Psi.calcRatio(W, iat);
1325 //Psi.rejectMove(iat);
1326 W.rejectMove(iat);
1327 fout << " ratio1 = " << ratio1 << " ration2 = " << ratio2 << std::endl;
1328 }
1329 }
1330 }
1331 //for(it=W.begin();it != it_end; ++it)
1332 //{
1333 // Walker_t& thisWalker(**it);
1334 // Walker_t::WFBuffer_t& w_buffer((*it)->DataSet);
1335 // w_buffer.rewind();
1336 // W.updateBuffer(**it,w_buffer);
1337 // RealType logpsi=Psi.updateBuffer(W,w_buffer,true);
1338 //}
1339 }
1340
1341 template<typename T, unsigned D>
randomize(ParticleAttrib<TinyVector<T,D>> & displ,T fac)1342 inline void randomize(ParticleAttrib<TinyVector<T, D>>& displ, T fac)
1343 {
1344 T* rv = &(displ[0][0]);
1345 assignUniformRand(rv, displ.size() * D, Random);
1346 for (int i = 0; i < displ.size() * D; ++i)
1347 rv[i] = fac * (rv[i] - 0.5);
1348 }
1349
runRatioV()1350 void WaveFunctionTester::runRatioV()
1351 {
1352 #if 0
1353 app_log() << "WaveFunctionTester::runRatioV " << std::endl;
1354 int nat = W.getTotalNum();
1355 Tau=0.025;
1356
1357 //create a VP with 8 virtual moves
1358 VirtualParticleSet vp(&W,8);
1359 W.enableVirtualMoves();
1360
1361 //cheating
1362 const ParticleSet& ions=W.DistTables[1]->origin();
1363 DistanceTableData* dt_ie=W.DistTables[1];
1364 double Rmax=2.0;
1365
1366 ParticleSet::ParticlePos_t sphere(8);
1367 std::vector<RealType> ratio_1(8), ratio_v(8);
1368 MCWalkerConfiguration::iterator it(W.begin()), it_end(W.end());
1369 while (it != it_end)
1370 {
1371 makeGaussRandom(deltaR);
1372 Walker_t::WFBuffer_t tbuffer;
1373 W.R = (**it).R+Tau*deltaR;
1374 (**it).R=W.R;
1375 W.update();
1376 RealType logpsi=Psi.registerData(W,tbuffer);
1377
1378 W.initVirtualMoves();
1379
1380 for(int iat=0; iat<ions.getTotalNum(); ++iat)
1381 {
1382 for(int nn=dt_ie->M[iat],iel=0; nn<dt_ie->M[iat+1]; nn++,iel++)
1383 {
1384 RealType r(dt_ie->r(nn));
1385 if(r>Rmax) continue;
1386 randomize(sphere,(RealType)0.5);
1387
1388 for(int k=0; k<sphere.size(); ++k)
1389 {
1390 W.makeMoveOnSphere(iel,sphere[k]);
1391 ratio_1[k]=Psi.ratio(W,iel);
1392 W.rejectMove(iel);
1393 Psi.resetPhaseDiff();
1394 }
1395
1396 vp.makeMoves(iel,sphere);
1397
1398 Psi.evaluateRatios(vp,ratio_v);
1399
1400 app_log() << "IAT = " << iat << " " << iel << std::endl;
1401 for(int k=0; k<sphere.size(); ++k)
1402 {
1403 app_log() << ratio_1[k]/ratio_v[k] << " " << ratio_1[k] << std::endl;
1404 }
1405 app_log() << std::endl;
1406 }
1407 }
1408 ++it;
1409 }
1410 #endif
1411 }
1412
runGradSourceTest()1413 void WaveFunctionTester::runGradSourceTest()
1414 {
1415
1416 app_log() << " ===== runGradSourceTest =====\n";
1417 ParticleSetPool::PoolType::iterator p;
1418 for (p = PtclPool.getPool().begin(); p != PtclPool.getPool().end(); p++)
1419 app_log() << "ParticelSet = " << p->first << std::endl;
1420 // Find source ParticleSet
1421 ParticleSetPool::PoolType::iterator pit(PtclPool.getPool().find(sourceName));
1422 app_log() << pit->first << std::endl;
1423 // if(pit == PtclPool.getPool().end())
1424 // APP_ABORT("Unknown source \"" + sourceName + "\" WaveFunctionTester.");
1425 ParticleSet& source = *((*pit).second);
1426 RealType delta = 0.00001;
1427 ValueType c1 = 1.0 / delta / 2.0;
1428 ValueType c2 = 1.0 / delta / delta;
1429 int nat = W.getTotalNum();
1430 ParticleSet::ParticlePos_t deltaR(nat);
1431 MCWalkerConfiguration::PropertyContainer_t Properties(0,0,1,WP::MAXPROPERTIES);;
1432 //pick the first walker
1433 MCWalkerConfiguration::Walker_t* awalker = *(W.begin());
1434 //copy the properties of the working walker
1435 Properties = awalker->Properties;
1436 //sample a new walker configuration and copy to ParticleSet::R
1437 //makeGaussRandom(deltaR);
1438 W.R = awalker->R;
1439 //W.R += deltaR;
1440 W.update();
1441 //ValueType psi = Psi.evaluate(W);
1442 ValueType logpsi = Psi.evaluateLog(W);
1443 RealType eloc = H.evaluate(W);
1444 H.auxHevaluate(W);
1445 app_log() << " HamTest "
1446 << " Total " << eloc << std::endl;
1447 for (int i = 0; i < H.sizeOfObservables(); i++)
1448 app_log() << " HamTest " << H.getObservableName(i) << " " << H.getObservable(i) << std::endl;
1449 //RealType psi = Psi.evaluateLog(W);
1450 ParticleSet::ParticleGradient_t G(nat), G1(nat);
1451 ParticleSet::ParticleLaplacian_t L(nat), L1(nat);
1452 G = W.G;
1453 L = W.L;
1454
1455 //This code computes d/dR \ln \Psi_T using the evalGradSource method, and
1456 // by finite difference. Results are saved in grad_ion and grad_ion_FD respectively.
1457 // GRAD TEST COMPUTATION
1458 int nions = source.getTotalNum();
1459 ParticleSet::ParticleGradient_t grad_ion(nions), grad_ion_FD(nions);
1460 for (int iat = 0; iat < nions; iat++)
1461 {
1462 grad_ion[iat] = Psi.evalGradSource(W, source, iat);
1463 PosType rI = source.R[iat];
1464 for (int iondim = 0; iondim < 3; iondim++)
1465 {
1466 source.R[iat][iondim] = rI[iondim] + delta;
1467 source.update();
1468 W.update();
1469
1470 ValueType log_p = Psi.evaluateLog(W);
1471
1472 source.R[iat][iondim] = rI[iondim] - delta;
1473 source.update();
1474 W.update();
1475 ValueType log_m = Psi.evaluateLog(W);
1476
1477 //symmetric finite difference formula for gradient.
1478 grad_ion_FD[iat][iondim] = c1 * (log_p - log_m);
1479
1480 //reset everything to how it was.
1481 source.R[iat][iondim] = rI[iondim];
1482 source.update();
1483 W.update();
1484 }
1485 //this lastone makes sure the distance tables correspond to unperturbed source.
1486 }
1487 //END GRAD TEST COMPUTATION
1488
1489 for (int isrc = 0; isrc < 1 /*source.getTotalNum()*/; isrc++)
1490 {
1491 TinyVector<ParticleSet::ParticleGradient_t, OHMMS_DIM> grad_grad;
1492 TinyVector<ParticleSet::ParticleLaplacian_t, OHMMS_DIM> lapl_grad;
1493 TinyVector<ParticleSet::ParticleGradient_t, OHMMS_DIM> grad_grad_FD;
1494 TinyVector<ParticleSet::ParticleLaplacian_t, OHMMS_DIM> lapl_grad_FD;
1495 for (int dim = 0; dim < OHMMS_DIM; dim++)
1496 {
1497 grad_grad[dim].resize(nat);
1498 lapl_grad[dim].resize(nat);
1499 grad_grad_FD[dim].resize(nat);
1500 lapl_grad_FD[dim].resize(nat);
1501 }
1502 Psi.evaluateLog(W);
1503 GradType grad_log = Psi.evalGradSource(W, source, isrc, grad_grad, lapl_grad);
1504 ValueType log = Psi.evaluateLog(W);
1505 //grad_log = Psi.evalGradSource (W, source, isrc);
1506 for (int iat = 0; iat < nat; iat++)
1507 {
1508 PosType r0 = W.R[iat];
1509 GradType gFD[OHMMS_DIM];
1510 GradType lapFD = ValueType();
1511 for (int eldim = 0; eldim < 3; eldim++)
1512 {
1513 W.R[iat][eldim] = r0[eldim] + delta;
1514 W.update();
1515 ValueType log_p = Psi.evaluateLog(W);
1516 GradType gradlogpsi_p = Psi.evalGradSource(W, source, isrc);
1517 W.R[iat][eldim] = r0[eldim] - delta;
1518 W.update();
1519 ValueType log_m = Psi.evaluateLog(W);
1520 GradType gradlogpsi_m = Psi.evalGradSource(W, source, isrc);
1521 lapFD += gradlogpsi_m + gradlogpsi_p;
1522 gFD[eldim] = gradlogpsi_p - gradlogpsi_m;
1523 W.R[iat] = r0;
1524 W.update();
1525 //Psi.evaluateLog(W);
1526 }
1527 const ValueType six(6);
1528 for (int iondim = 0; iondim < OHMMS_DIM; iondim++)
1529 {
1530 for (int eldim = 0; eldim < OHMMS_DIM; eldim++)
1531 grad_grad_FD[iondim][iat][eldim] = c1 * gFD[eldim][iondim];
1532 lapl_grad_FD[iondim][iat] = c2 * (lapFD[iondim] - six * grad_log[iondim]);
1533 }
1534 }
1535 for (int dimsrc = 0; dimsrc < OHMMS_DIM; dimsrc++)
1536 {
1537 for (int iat = 0; iat < nat; iat++)
1538 {
1539 fout << "For particle #" << iat << " at " << W.R[iat] << std::endl;
1540 fout << "Grad Gradient = " << std::setw(12) << grad_grad[dimsrc][iat] << std::endl
1541 << " Finite diff = " << std::setw(12) << grad_grad_FD[dimsrc][iat] << std::endl
1542 << " Error = " << std::setw(12) << grad_grad_FD[dimsrc][iat] - grad_grad[dimsrc][iat] << std::endl
1543 << std::endl;
1544 fout << "Grad Laplacian = " << std::setw(12) << lapl_grad[dimsrc][iat] << std::endl
1545 << " Finite diff = " << std::setw(12) << lapl_grad_FD[dimsrc][iat] << std::endl
1546 << " Error = " << std::setw(12) << lapl_grad_FD[dimsrc][iat] - lapl_grad[dimsrc][iat] << std::endl
1547 << std::endl;
1548 }
1549 }
1550 fout << "==== BEGIN Ion Gradient Check ====\n";
1551 for (int iat = 0; iat < nions; iat++)
1552 {
1553 fout << "For ion #" << iat << " at " << source.R[iat] << std::endl;
1554 fout << "Gradient = " << std::setw(12) << grad_ion[iat] << std::endl
1555 << " Finite diff = " << std::setw(12) << grad_ion_FD[iat] << std::endl
1556 << " Error = " << std::setw(12) << grad_ion_FD[iat] - grad_ion[iat] << std::endl
1557 << std::endl;
1558 }
1559 fout << "==== END Ion Gradient Check ====\n";
1560 fout.flush();
1561 }
1562 }
1563
1564
runZeroVarianceTest()1565 void WaveFunctionTester::runZeroVarianceTest()
1566 {
1567 app_log() << " ===== runZeroVarianceTest =====\n";
1568 ParticleSetPool::PoolType::iterator p;
1569 for (p = PtclPool.getPool().begin(); p != PtclPool.getPool().end(); p++)
1570 app_log() << "ParticelSet = " << p->first << std::endl;
1571 // Find source ParticleSet
1572 ParticleSetPool::PoolType::iterator pit(PtclPool.getPool().find(sourceName));
1573 app_log() << pit->first << std::endl;
1574 // if(pit == PtclPool.getPool().end())
1575 // APP_ABORT("Unknown source \"" + sourceName + "\" WaveFunctionTester.");
1576 ParticleSet& source = *((*pit).second);
1577 int nat = W.getTotalNum();
1578 ParticleSet::ParticlePos_t deltaR(nat);
1579 MCWalkerConfiguration::PropertyContainer_t Properties(0,0,1,WP::MAXPROPERTIES);;
1580 //pick the first walker
1581 MCWalkerConfiguration::Walker_t* awalker = *(W.begin());
1582 //copy the properties of the working walker
1583 Properties = awalker->Properties;
1584 //sample a new walker configuration and copy to ParticleSet::R
1585 //makeGaussRandom(deltaR);
1586 W.R = awalker->R;
1587 //W.R += deltaR;
1588 W.update();
1589 //ValueType psi = Psi.evaluate(W);
1590 ValueType logpsi = Psi.evaluateLog(W);
1591 RealType eloc = H.evaluate(W);
1592 //RealType psi = Psi.evaluateLog(W);
1593 ParticleSet::ParticleGradient_t G(nat), G1(nat);
1594 ParticleSet::ParticleLaplacian_t L(nat), L1(nat);
1595 G = W.G;
1596 L = W.L;
1597 PosType r1(5.0, 2.62, 2.55);
1598 W.R[1] = PosType(4.313, 5.989, 4.699);
1599 W.R[2] = PosType(5.813, 4.321, 4.893);
1600 W.R[3] = PosType(4.002, 5.502, 5.381);
1601 W.R[4] = PosType(5.901, 5.121, 5.311);
1602 W.R[5] = PosType(5.808, 4.083, 5.021);
1603 W.R[6] = PosType(4.750, 5.810, 4.732);
1604 W.R[7] = PosType(4.690, 5.901, 4.989);
1605 for (int i = 1; i < 8; i++)
1606 W.R[i] -= PosType(2.5, 2.5, 2.5);
1607 char fname[32];
1608 sprintf(fname, "ZVtest.%03d.dat", OHMMS::Controller->rank());
1609 FILE* fzout = fopen(fname, "w");
1610 TinyVector<ParticleSet::ParticleGradient_t, OHMMS_DIM> grad_grad;
1611 TinyVector<ParticleSet::ParticleLaplacian_t, OHMMS_DIM> lapl_grad;
1612 TinyVector<ParticleSet::ParticleGradient_t, OHMMS_DIM> grad_grad_FD;
1613 TinyVector<ParticleSet::ParticleLaplacian_t, OHMMS_DIM> lapl_grad_FD;
1614 for (int dim = 0; dim < OHMMS_DIM; dim++)
1615 {
1616 grad_grad[dim].resize(nat);
1617 lapl_grad[dim].resize(nat);
1618 grad_grad_FD[dim].resize(nat);
1619 lapl_grad_FD[dim].resize(nat);
1620 }
1621 for (r1[0] = 0.0; r1[0] < 5.0; r1[0] += 1.0e-4)
1622 {
1623 W.R[0] = r1;
1624 fprintf(fzout, "%1.8e %1.8e %1.8e ", r1[0], r1[1], r1[2]);
1625 RealType log = Psi.evaluateLog(W);
1626 // ValueType psi = std::cos(Psi.getPhase())*std::exp(log);//*W.PropertyList[SIGN];
1627 #if defined(QMC_COMPLEX)
1628 RealType ratioMag = std::exp(log);
1629 ValueType psi =
1630 std::complex<OHMMS_PRECISION>(ratioMag * std::cos(Psi.getPhase()), ratioMag * std::sin(Psi.getPhase()));
1631 #else
1632 ValueType psi = std::cos(Psi.getPhase()) * std::exp(log); //*W.PropertyList[SIGN];
1633 #endif
1634 double E = H.evaluate(W);
1635 //double KE = E - W.PropertyList[LOCALPOTENTIAL];
1636 double KE = -0.5 * (Sum(W.L) + Dot(W.G, W.G));
1637 #if defined(QMC_COMPLEX)
1638 fprintf(fzout, "%16.12e %16.12e %16.12e ", psi.real(), psi.imag(), KE);
1639 #else
1640 fprintf(fzout, "%16.12e %16.12e ", psi, KE);
1641 #endif
1642 for (int isrc = 0; isrc < source.getTotalNum(); isrc++)
1643 {
1644 GradType grad_log = Psi.evalGradSource(W, source, isrc, grad_grad, lapl_grad);
1645 for (int dim = 0; dim < OHMMS_DIM; dim++)
1646 {
1647 double ZV = 0.5 * Sum(lapl_grad[dim]) + Dot(grad_grad[dim], W.G);
1648 #if defined(QMC_COMPLEX)
1649 fprintf(fzout, "%16.12e %16.12e %16.12e ", ZV, grad_log[dim].real(), grad_log[dim].imag());
1650 #else
1651 fprintf(fzout, "%16.12e %16.12e ", ZV, grad_log[dim]);
1652 #endif
1653 }
1654 }
1655 fprintf(fzout, "\n");
1656 }
1657 fclose(fzout);
1658 }
1659
1660
put(xmlNodePtr q)1661 bool WaveFunctionTester::put(xmlNodePtr q)
1662 {
1663 myNode = q;
1664 xmlNodePtr tcur = q->children;
1665 while (tcur != NULL)
1666 {
1667 std::string cname((const char*)(tcur->name));
1668 if (cname == "delta_output")
1669 {
1670 outputDeltaVsError = true;
1671 DeltaVsError.put(tcur);
1672 }
1673 tcur = tcur->next;
1674 }
1675
1676 bool success = putQMCInfo(q);
1677
1678 return success;
1679 }
1680
runDerivTest()1681 void WaveFunctionTester::runDerivTest()
1682 {
1683 app_log() << " ===== runDerivTest =====\n";
1684 app_log() << " Testing derivatives" << std::endl;
1685 int nat = W.getTotalNum();
1686 MCWalkerConfiguration::PropertyContainer_t Properties(0,0,1,WP::MAXPROPERTIES);;
1687 //pick the first walker
1688 MCWalkerConfiguration::Walker_t* awalker = *(W.begin());
1689 //copy the properties of the working walker
1690 Properties = awalker->Properties;
1691 //sample a new walker configuration and copy to ParticleSet::R
1692 W.R = awalker->R + deltaR;
1693
1694 fout << "Position " << std::endl << W.R << std::endl;
1695
1696 //W.R += deltaR;
1697 W.update();
1698 //ValueType psi = Psi.evaluate(W);
1699 ValueType logpsi = Psi.evaluateLog(W);
1700 RealType eloc = H.evaluate(W);
1701 app_log() << " HamTest "
1702 << " Total " << eloc << std::endl;
1703 for (int i = 0; i < H.sizeOfObservables(); i++)
1704 app_log() << " HamTest " << H.getObservableName(i) << " " << H.getObservable(i) << std::endl;
1705 //RealType psi = Psi.evaluateLog(W);
1706 ParticleSet::ParticleGradient_t G(nat), G1(nat);
1707 ParticleSet::ParticleLaplacian_t L(nat), L1(nat);
1708 G = W.G;
1709 L = W.L;
1710 fout << "Gradients" << std::endl;
1711 for (int iat = 0; iat < W.R.size(); iat++)
1712 {
1713 for (int i = 0; i < 3; i++)
1714 fout << W.G[iat][i] << " ";
1715 fout << std::endl;
1716 }
1717 fout << "Laplaians" << std::endl;
1718 for (int iat = 0; iat < W.R.size(); iat++)
1719 {
1720 fout << W.L[iat] << " ";
1721 fout << std::endl;
1722 }
1723 opt_variables_type wfVars, wfvar_prime;
1724 //build optimizables from the wavefunction
1725 wfVars.clear();
1726 Psi.checkInVariables(wfVars);
1727 wfVars.resetIndex();
1728 Psi.checkOutVariables(wfVars);
1729 wfvar_prime = wfVars;
1730 wfVars.print(fout);
1731 int Nvars = wfVars.size();
1732 std::vector<ValueType> Dsaved(Nvars);
1733 std::vector<ValueType> HDsaved(Nvars);
1734 std::vector<RealType> PGradient(Nvars);
1735 std::vector<RealType> HGradient(Nvars);
1736 Psi.resetParameters(wfVars);
1737 logpsi = Psi.evaluateLog(W);
1738
1739 //reuse the sphere
1740 H.setPrimary(false);
1741
1742 eloc = H.evaluate(W);
1743 Psi.evaluateDerivatives(W, wfVars, Dsaved, HDsaved);
1744 RealType FiniteDiff = 1e-6;
1745 QMCTraits::RealType dh = 1.0 / (2.0 * FiniteDiff);
1746 for (int i = 0; i < Nvars; i++)
1747 {
1748 for (int j = 0; j < Nvars; j++)
1749 wfvar_prime[j] = wfVars[j];
1750 wfvar_prime[i] = wfVars[i] + FiniteDiff;
1751 // Psi.checkOutVariables(wfvar_prime);
1752 Psi.resetParameters(wfvar_prime);
1753 Psi.reset();
1754 W.update();
1755 W.G = 0;
1756 W.L = 0;
1757 RealType logpsiPlus = Psi.evaluateLog(W);
1758 H.evaluate(W);
1759 RealType elocPlus = H.getLocalEnergy() - H.getLocalPotential();
1760 wfvar_prime[i] = wfVars[i] - FiniteDiff;
1761 // Psi.checkOutVariables(wfvar_prime);
1762 Psi.resetParameters(wfvar_prime);
1763 Psi.reset();
1764 W.update();
1765 W.G = 0;
1766 W.L = 0;
1767 RealType logpsiMinus = Psi.evaluateLog(W);
1768 H.evaluate(W);
1769 RealType elocMinus = H.getLocalEnergy() - H.getLocalPotential();
1770 PGradient[i] = (logpsiPlus - logpsiMinus) * dh;
1771 HGradient[i] = (elocPlus - elocMinus) * dh;
1772 }
1773 Psi.resetParameters(wfVars);
1774 fout << std::endl << "Deriv Numeric Analytic" << std::endl;
1775 for (int i = 0; i < Nvars; i++)
1776 fout << i << " " << PGradient[i] << " " << std::real(Dsaved[i]) << " " << (PGradient[i] - std::real(Dsaved[i]))
1777 << std::endl;
1778 fout << std::endl << "Hderiv Numeric Analytic" << std::endl;
1779 for (int i = 0; i < Nvars; i++)
1780 fout << i << " " << HGradient[i] << " " << std::real(HDsaved[i]) << " " << (HGradient[i] - std::real(HDsaved[i]))
1781 << std::endl;
1782 }
1783
1784
runDerivNLPPTest()1785 void WaveFunctionTester::runDerivNLPPTest()
1786 {
1787 app_log() << " ===== runDerivNLPPTest =====\n";
1788 char fname[16];
1789 sprintf(fname, "nlpp.%03d", OHMMS::Controller->rank());
1790 std::ofstream nlout(fname);
1791 nlout.precision(15);
1792
1793 app_log() << " Testing derivatives" << std::endl;
1794 int nat = W.getTotalNum();
1795 MCWalkerConfiguration::PropertyContainer_t Properties(0,0,1,WP::MAXPROPERTIES);;
1796 //pick the first walker
1797 MCWalkerConfiguration::Walker_t* awalker = *(W.begin());
1798 //copy the properties of the working walker
1799 Properties = awalker->Properties;
1800 //sample a new walker configuration and copy to ParticleSet::R
1801 W.R = awalker->R + deltaR;
1802
1803 //W.R += deltaR;
1804 W.update();
1805 //ValueType psi = Psi.evaluate(W);
1806 ValueType logpsi = Psi.evaluateLog(W);
1807 RealType eloc = H.evaluate(W);
1808
1809 app_log() << " HamTest "
1810 << " Total " << eloc << std::endl;
1811 for (int i = 0; i < H.sizeOfObservables(); i++)
1812 app_log() << " HamTest " << H.getObservableName(i) << " " << H.getObservable(i) << std::endl;
1813
1814 //RealType psi = Psi.evaluateLog(W);
1815 ParticleSet::ParticleGradient_t G(nat), G1(nat);
1816 ParticleSet::ParticleLaplacian_t L(nat), L1(nat);
1817 G = W.G;
1818 L = W.L;
1819 nlout << "Gradients" << std::endl;
1820 for (int iat = 0; iat < W.R.size(); iat++)
1821 {
1822 for (int i = 0; i < 3; i++)
1823 nlout << W.G[iat][i] << " ";
1824 nlout << std::endl;
1825 }
1826 nlout << "Laplaians" << std::endl;
1827 for (int iat = 0; iat < W.R.size(); iat++)
1828 {
1829 nlout << W.L[iat] << " ";
1830 nlout << std::endl;
1831 }
1832 opt_variables_type wfVars, wfvar_prime;
1833 //build optimizables from the wavefunction
1834 wfVars.clear();
1835 Psi.checkInVariables(wfVars);
1836 wfVars.resetIndex();
1837 Psi.checkOutVariables(wfVars);
1838 wfvar_prime = wfVars;
1839 wfVars.print(nlout);
1840 int Nvars = wfVars.size();
1841 std::vector<ValueType> Dsaved(Nvars);
1842 std::vector<ValueType> HDsaved(Nvars);
1843 std::vector<RealType> PGradient(Nvars);
1844 std::vector<RealType> HGradient(Nvars);
1845 Psi.resetParameters(wfVars);
1846
1847 logpsi = Psi.evaluateLog(W);
1848
1849 //reuse the sphere for non-local pp
1850 H.setPrimary(false);
1851
1852 std::vector<RealType> ene(4), ene_p(4), ene_m(4);
1853 Psi.evaluateDerivatives(W, wfVars, Dsaved, HDsaved);
1854
1855 ene[0] = H.evaluateValueAndDerivatives(W, wfVars, Dsaved, HDsaved, true);
1856 app_log() << "Check the energy " << eloc << " " << H.getLocalEnergy() << " " << ene[0] << std::endl;
1857
1858 RealType FiniteDiff = 1e-6;
1859 QMCTraits::RealType dh = 1.0 / (2.0 * FiniteDiff);
1860 for (int i = 0; i < Nvars; i++)
1861 {
1862 for (int j = 0; j < Nvars; j++)
1863 wfvar_prime[j] = wfVars[j];
1864 wfvar_prime[i] = wfVars[i] + FiniteDiff;
1865 Psi.resetParameters(wfvar_prime);
1866 Psi.reset();
1867 W.update();
1868 W.G = 0;
1869 W.L = 0;
1870 RealType logpsiPlus = Psi.evaluateLog(W);
1871 RealType elocPlus = H.evaluateVariableEnergy(W, true);
1872
1873 //H.evaluate(W);
1874 //RealType elocPlus=H.getLocalEnergy()-H.getLocalPotential();
1875
1876 wfvar_prime[i] = wfVars[i] - FiniteDiff;
1877 Psi.resetParameters(wfvar_prime);
1878 Psi.reset();
1879 W.update();
1880 W.G = 0;
1881 W.L = 0;
1882 RealType logpsiMinus = Psi.evaluateLog(W);
1883 RealType elocMinus = H.evaluateVariableEnergy(W, true);
1884
1885 //H.evaluate(W);
1886 //RealType elocMinus = H.getLocalEnergy()-H.getLocalPotential();
1887
1888 PGradient[i] = (logpsiPlus - logpsiMinus) * dh;
1889 HGradient[i] = (elocPlus - elocMinus) * dh;
1890 }
1891
1892 nlout << std::endl << "Deriv Numeric Analytic" << std::endl;
1893 for (int i = 0; i < Nvars; i++)
1894 nlout << i << " " << PGradient[i] << " " << Dsaved[i] << " " << (PGradient[i] - Dsaved[i]) << std::endl;
1895 nlout << std::endl << "Hderiv Numeric Analytic" << std::endl;
1896 for (int i = 0; i < Nvars; i++)
1897 nlout << i << " " << HGradient[i] << " " << HDsaved[i] << " " << (HGradient[i] - HDsaved[i]) << std::endl;
1898 }
1899
1900
runDerivCloneTest()1901 void WaveFunctionTester::runDerivCloneTest()
1902 {
1903 app_log() << " ===== runDerivCloneTest =====\n";
1904 app_log() << " Testing derivatives clone" << std::endl;
1905 RandomGenerator_t* Rng1 = new RandomGenerator_t();
1906 RandomGenerator_t* Rng2 = new RandomGenerator_t();
1907 (*Rng1) = (*Rng2);
1908 MCWalkerConfiguration* w_clone = new MCWalkerConfiguration(W);
1909 TrialWaveFunction* psi_clone = Psi.makeClone(*w_clone);
1910 QMCHamiltonian* h_clone = H.makeClone(*w_clone, *psi_clone);
1911 h_clone->setRandomGenerator(Rng2);
1912 H.setRandomGenerator(Rng1);
1913 h_clone->setPrimary(true);
1914 int nat = W.getTotalNum();
1915 ParticleSet::ParticlePos_t deltaR(nat);
1916 //pick the first walker
1917 MCWalkerConfiguration::Walker_t* awalker = *(W.begin());
1918 // MCWalkerConfiguration::Walker_t* bwalker = *(w_clone->begin());
1919 // bwalker->R = awalker->R;
1920 W.R = awalker->R;
1921 W.update();
1922 w_clone->R = awalker->R;
1923 w_clone->update();
1924 opt_variables_type wfVars;
1925 //build optimizables from the wavefunction
1926 // wfVars.clear();
1927 Psi.checkInVariables(wfVars);
1928 wfVars.resetIndex();
1929 Psi.checkOutVariables(wfVars);
1930 wfVars.print(fout);
1931 int Nvars = wfVars.size();
1932 opt_variables_type wfvar_prime;
1933 // wfvar_prime.insertFrom(wfVars);
1934 // wfvar_prime.clear();
1935 psi_clone->checkInVariables(wfvar_prime);
1936 wfvar_prime.resetIndex();
1937 for (int j = 0; j < Nvars; j++)
1938 wfvar_prime[j] = wfVars[j];
1939 psi_clone->checkOutVariables(wfvar_prime);
1940 wfvar_prime.print(fout);
1941 psi_clone->resetParameters(wfvar_prime);
1942 Psi.resetParameters(wfVars);
1943 std::vector<ValueType> Dsaved(Nvars, 0), og_Dsaved(Nvars, 0);
1944 std::vector<ValueType> HDsaved(Nvars, 0), og_HDsaved(Nvars, 0);
1945 std::vector<RealType> PGradient(Nvars, 0), og_PGradient(Nvars, 0);
1946 std::vector<RealType> HGradient(Nvars, 0), og_HGradient(Nvars, 0);
1947 ValueType logpsi2 = psi_clone->evaluateLog(*w_clone);
1948 RealType eloc2 = h_clone->evaluate(*w_clone);
1949 psi_clone->evaluateDerivatives(*w_clone, wfvar_prime, Dsaved, HDsaved);
1950 ValueType logpsi1 = Psi.evaluateLog(W);
1951 RealType eloc1 = H.evaluate(W);
1952 Psi.evaluateDerivatives(W, wfVars, og_Dsaved, og_HDsaved);
1953 app_log() << "log (original) = " << logpsi1 << " energy = " << eloc1 << std::endl;
1954 for (int i = 0; i < H.sizeOfObservables(); i++)
1955 app_log() << " HamTest " << H.getObservableName(i) << " " << H.getObservable(i) << std::endl;
1956 app_log() << "log (clone) = " << logpsi2 << " energy = " << eloc2 << std::endl;
1957 for (int i = 0; i < h_clone->sizeOfObservables(); i++)
1958 app_log() << " HamTest " << h_clone->getObservableName(i) << " " << h_clone->getObservable(i) << std::endl;
1959 // app_log()<<" Saved quantities:"<< std::endl;
1960 // for(int i=0;i<Nvars;i++) app_log()<<Dsaved[i]<<" "<<og_Dsaved[i]<<" "<<HDsaved[i]<<" "<<og_HDsaved[i]<< std::endl;
1961 RealType FiniteDiff = 1e-6;
1962 QMCTraits::RealType dh = 1.0 / (2.0 * FiniteDiff);
1963 for (int i = 0; i < Nvars; i++)
1964 {
1965 for (int j = 0; j < Nvars; j++)
1966 wfvar_prime[j] = wfVars[j];
1967 wfvar_prime[i] = wfVars[i] + FiniteDiff;
1968 psi_clone->resetParameters(wfvar_prime);
1969 psi_clone->reset();
1970 w_clone->update();
1971 w_clone->G = 0;
1972 w_clone->L = 0;
1973 RealType logpsiPlus = psi_clone->evaluateLog(*w_clone);
1974 h_clone->evaluate(*w_clone);
1975 RealType elocPlus = h_clone->getLocalEnergy() - h_clone->getLocalPotential();
1976 wfvar_prime[i] = wfVars[i] - FiniteDiff;
1977 psi_clone->resetParameters(wfvar_prime);
1978 psi_clone->reset();
1979 w_clone->update();
1980 w_clone->G = 0;
1981 w_clone->L = 0;
1982 RealType logpsiMinus = psi_clone->evaluateLog(*w_clone);
1983 h_clone->evaluate(*w_clone);
1984 RealType elocMinus = h_clone->getLocalEnergy() - h_clone->getLocalPotential();
1985 PGradient[i] = (logpsiPlus - logpsiMinus) * dh;
1986 HGradient[i] = (elocPlus - elocMinus) * dh;
1987 }
1988 fout << "CLONE" << std::endl;
1989 fout << std::endl << " Deriv Numeric Analytic" << std::endl;
1990 for (int i = 0; i < Nvars; i++)
1991 fout << i << " " << PGradient[i] << " " << std::real(Dsaved[i]) << " "
1992 << (PGradient[i] - std::real(Dsaved[i])) / PGradient[i] << std::endl;
1993 fout << std::endl << " Hderiv Numeric Analytic" << std::endl;
1994 for (int i = 0; i < Nvars; i++)
1995 fout << i << " " << HGradient[i] << " " << std::real(HDsaved[i]) << " "
1996 << (HGradient[i] - std::real(HDsaved[i])) / HGradient[i] << std::endl;
1997 for (int i = 0; i < Nvars; i++)
1998 {
1999 for (int j = 0; j < Nvars; j++)
2000 wfvar_prime[j] = wfVars[j];
2001 wfvar_prime[i] = wfVars[i] + FiniteDiff;
2002 Psi.resetParameters(wfvar_prime);
2003 Psi.reset();
2004 W.update();
2005 W.G = 0;
2006 W.L = 0;
2007 RealType logpsiPlus = Psi.evaluateLog(W);
2008 H.evaluate(W);
2009 RealType elocPlus = H.getLocalEnergy() - H.getLocalPotential();
2010 wfvar_prime[i] = wfVars[i] - FiniteDiff;
2011 Psi.resetParameters(wfvar_prime);
2012 Psi.reset();
2013 W.update();
2014 W.G = 0;
2015 W.L = 0;
2016 RealType logpsiMinus = Psi.evaluateLog(W);
2017 H.evaluate(W);
2018 RealType elocMinus = H.getLocalEnergy() - H.getLocalPotential();
2019 PGradient[i] = (logpsiPlus - logpsiMinus) * dh;
2020 HGradient[i] = (elocPlus - elocMinus) * dh;
2021 }
2022 fout << "ORIGINAL" << std::endl;
2023 fout << std::endl << " Deriv Numeric Analytic" << std::endl;
2024 for (int i = 0; i < Nvars; i++)
2025 fout << i << " " << PGradient[i] << " " << std::real(Dsaved[i]) << " "
2026 << (PGradient[i] - std::real(Dsaved[i])) / PGradient[i] << std::endl;
2027 fout << std::endl << " Hderiv Numeric Analytic" << std::endl;
2028 for (int i = 0; i < Nvars; i++)
2029 fout << i << " " << HGradient[i] << " " << std::real(HDsaved[i]) << " "
2030 << (HGradient[i] - std::real(HDsaved[i])) / HGradient[i] << std::endl;
2031 }
runwftricks()2032 void WaveFunctionTester::runwftricks()
2033 {
2034 std::vector<WaveFunctionComponent*>& Orbitals = Psi.getOrbitals();
2035 app_log() << " Total of " << Orbitals.size() << " orbitals." << std::endl;
2036 int SDindex(0);
2037 for (int i = 0; i < Orbitals.size(); i++)
2038 if ("SlaterDet" == Orbitals[i]->ClassName)
2039 SDindex = i;
2040 SPOSetPtr Phi = dynamic_cast<SlaterDet*>(Orbitals[SDindex])->getPhi();
2041 int NumOrbitals = Phi->getBasisSetSize();
2042 app_log() << "Basis set size: " << NumOrbitals << std::endl;
2043 std::vector<int> SPONumbers(0, 0);
2044 std::vector<int> irrepRotations(0, 0);
2045 std::vector<int> Grid(0, 0);
2046 xmlNodePtr kids = myNode->children;
2047 std::string doProj("yes");
2048 std::string doRotate("yes");
2049 std::string sClass("C2V");
2050 ParameterSet aAttrib;
2051 aAttrib.add(doProj, "projection");
2052 aAttrib.add(doRotate, "rotate");
2053 aAttrib.put(myNode);
2054 while (kids != NULL)
2055 {
2056 std::string cname((const char*)(kids->name));
2057 if (cname == "orbitals")
2058 {
2059 putContent(SPONumbers, kids);
2060 }
2061 else if (cname == "representations")
2062 {
2063 putContent(irrepRotations, kids);
2064 }
2065 else if (cname == "grid")
2066 putContent(Grid, kids);
2067 kids = kids->next;
2068 }
2069 ParticleSet::ParticlePos_t R_cart(1);
2070 R_cart.setUnit(PosUnit::Cartesian);
2071 ParticleSet::ParticlePos_t R_unit(1);
2072 R_unit.setUnit(PosUnit::Lattice);
2073 // app_log()<<" My crystals basis set is:"<< std::endl;
2074 // std::vector<std::vector<RealType> > BasisMatrix(3, std::vector<RealType>(3,0.0));
2075 //
2076 // for (int i=0;i<3;i++)
2077 // {
2078 // R_unit[0][0]=0;
2079 // R_unit[0][1]=0;
2080 // R_unit[0][2]=0;
2081 // R_unit[0][i]=1;
2082 // W.convert2Cart(R_unit,R_cart);
2083 // app_log()<<"basis_"<<i<<": ("<<R_cart[0][0]<<", "<<R_cart[0][1]<<", "<<R_cart[0][2]<<")"<< std::endl;
2084 // for (int j=0;j<3;j++) BasisMatrix[j][i]=R_cart[0][j];
2085 // }
2086 int Nrotated(SPONumbers.size());
2087 app_log() << " Projected orbitals: ";
2088 for (int i = 0; i < Nrotated; i++)
2089 app_log() << SPONumbers[i] << " ";
2090 app_log() << std::endl;
2091 //indexing trick
2092 // for(int i=0;i<Nrotated;i++) SPONumbers[i]-=1;
2093 SymmetryBuilder SO;
2094 SO.put(myNode);
2095 SymmetryGroup symOp(*SO.getSymmetryGroup());
2096 // SO.changeBasis(InverseBasisMatrix);
2097 OrbitalSetTraits<ValueType>::ValueVector_t values;
2098 values.resize(NumOrbitals);
2099 RealType overG0(1.0 / Grid[0]);
2100 RealType overG1(1.0 / Grid[1]);
2101 RealType overG2(1.0 / Grid[2]);
2102 std::vector<RealType> NormPhi(Nrotated, 0.0);
2103 int totsymops = symOp.getSymmetriesSize();
2104 Matrix<RealType> SymmetryOrbitalValues;
2105 SymmetryOrbitalValues.resize(Nrotated, totsymops);
2106 int ctabledim = symOp.getClassesSize();
2107 Matrix<double> projs(Nrotated, ctabledim);
2108 Matrix<double> orthoProjs(Nrotated, Nrotated);
2109 std::vector<RealType> brokenSymmetryCharacter(totsymops);
2110 for (int k = 0; k < Nrotated; k++)
2111 for (int l = 0; l < totsymops; l++)
2112 brokenSymmetryCharacter[l] += irrepRotations[k] * symOp.getsymmetryCharacter(l, irrepRotations[k] - 1);
2113 // app_log()<<"bsc: ";
2114 // for(int l=0;l<totsymops;l++) app_log()<<brokenSymmetryCharacter[l]<<" ";
2115 // app_log()<< std::endl;
2116 // for(int l=0;l<totsymops;l++) brokenSymmetryCharacter[l]+=0.5;
2117 if ((doProj == "yes") || (doRotate == "yes"))
2118 {
2119 OrbitalSetTraits<ValueType>::ValueVector_t identityValues(values.size());
2120 //Loop over grid
2121 for (int i = 0; i < Grid[0]; i++)
2122 for (int j = 0; j < Grid[1]; j++)
2123 for (int k = 0; k < Grid[2]; k++)
2124 {
2125 //Loop over symmetry classes and small group operators
2126 for (int l = 0; l < totsymops; l++)
2127 {
2128 R_unit[0][0] = overG0 * RealType(i); // R_cart[0][0]=0;
2129 R_unit[0][1] = overG1 * RealType(j); // R_cart[0][1]=0;
2130 R_unit[0][2] = overG2 * RealType(k); // R_cart[0][2]=0;
2131 // for(int a=0; a<3; a++) for(int b=0;b<3;b++) R_cart[0][a]+=BasisMatrix[a][b]*R_unit[0][b];
2132 W.convert2Cart(R_unit, R_cart);
2133 symOp.TransformSinglePosition(R_cart, l);
2134 W.R[0] = R_cart[0];
2135 values = 0.0;
2136 //evaluate orbitals
2137 // Phi->evaluate(W,0,values);
2138 Psi.evaluateLog(W);
2139 // YYYY: is the following two lines still maintained?
2140 //for(int n=0; n<NumOrbitals; n++)
2141 // values[n] = Phi->t_logpsi(0,n);
2142 if (l == 0)
2143 {
2144 identityValues = values;
2145 #if defined(QMC_COMPLEX)
2146 for (int n = 0; n < Nrotated; n++)
2147 NormPhi[n] += totsymops * real(values[SPONumbers[n]] * values[SPONumbers[n]]);
2148 #else
2149 for (int n = 0; n < Nrotated; n++)
2150 NormPhi[n] += totsymops * (values[SPONumbers[n]] * values[SPONumbers[n]]);
2151 #endif
2152 }
2153 //now we have phi evaluated at the rotated/inverted/whichever coordinates
2154 for (int n = 0; n < Nrotated; n++)
2155 {
2156 int N = SPONumbers[n];
2157 #if defined(QMC_COMPLEX)
2158 RealType phi2 = real(values[N] * identityValues[N]);
2159 #else
2160 RealType phi2 = (values[N] * identityValues[N]);
2161 #endif
2162 SymmetryOrbitalValues(n, l) += phi2;
2163 }
2164 for (int n = 0; n < Nrotated; n++)
2165 for (int p = 0; p < Nrotated; p++)
2166 {
2167 int N = SPONumbers[n];
2168 int P = SPONumbers[p];
2169 #if defined(QMC_COMPLEX)
2170 orthoProjs(n, p) += 0.5 * real(identityValues[N] * values[P] + identityValues[P] * values[N]) *
2171 brokenSymmetryCharacter[l];
2172 #else
2173 orthoProjs(n, p) +=
2174 0.5 * (identityValues[N] * values[P] + identityValues[P] * values[N]) * brokenSymmetryCharacter[l];
2175 #endif
2176 }
2177 }
2178 }
2179 for (int n = 0; n < Nrotated; n++)
2180 for (int l = 0; l < totsymops; l++)
2181 SymmetryOrbitalValues(n, l) /= NormPhi[n];
2182 for (int n = 0; n < Nrotated; n++)
2183 for (int l = 0; l < Nrotated; l++)
2184 orthoProjs(n, l) /= std::sqrt(NormPhi[n] * NormPhi[l]);
2185 // if (true){
2186 // app_log()<< std::endl;
2187 // for(int n=0;n<Nrotated;n++) {
2188 // for(int l=0;l<totsymops;l++) app_log()<<SymmetryOrbitalValues(n,l)<<" ";
2189 // app_log()<< std::endl;
2190 // }
2191 // app_log()<< std::endl;
2192 // }
2193 for (int n = 0; n < Nrotated; n++)
2194 {
2195 if (false)
2196 app_log() << " orbital #" << SPONumbers[n] << std::endl;
2197 for (int i = 0; i < ctabledim; i++)
2198 {
2199 double proj(0);
2200 for (int j = 0; j < totsymops; j++)
2201 proj += symOp.getsymmetryCharacter(j, i) * SymmetryOrbitalValues(n, j);
2202 if (false)
2203 app_log() << " Rep " << i << ": " << proj;
2204 projs(n, i) = proj < 1e-4 ? 0 : proj;
2205 }
2206 if (false)
2207 app_log() << std::endl;
2208 }
2209 if (true)
2210 {
2211 app_log() << "Printing Projection Matrix" << std::endl;
2212 for (int n = 0; n < Nrotated; n++)
2213 {
2214 for (int l = 0; l < ctabledim; l++)
2215 app_log() << projs(n, l) << " ";
2216 app_log() << std::endl;
2217 }
2218 app_log() << std::endl;
2219 }
2220 if (true)
2221 {
2222 app_log() << "Printing Coefficient Matrix" << std::endl;
2223 for (int n = 0; n < Nrotated; n++)
2224 {
2225 for (int l = 0; l < ctabledim; l++)
2226 app_log() << std::sqrt(projs(n, l)) << " ";
2227 app_log() << std::endl;
2228 }
2229 app_log() << std::endl;
2230 }
2231 if (doRotate == "yes")
2232 {
2233 // app_log()<<"Printing Broken Symmetry Projection Matrix"<< std::endl;
2234 // for(int n=0;n<Nrotated;n++) {
2235 // for(int l=0;l<Nrotated;l++) app_log()<<orthoProjs(n,l)<<" ";
2236 // app_log()<< std::endl;
2237 // }
2238 char JOBU('A');
2239 char JOBVT('A');
2240 int vdim = Nrotated;
2241 Vector<double> Sigma(vdim);
2242 Matrix<double> U(vdim, vdim);
2243 Matrix<double> VT(vdim, vdim);
2244 int lwork = 8 * Nrotated;
2245 std::vector<double> work(lwork, 0);
2246 int info(0);
2247 dgesvd(&JOBU, &JOBVT, &vdim, &vdim, orthoProjs.data(), &vdim, Sigma.data(), U.data(), &vdim, VT.data(), &vdim,
2248 &(work[0]), &lwork, &info);
2249 app_log() << "Printing Rotation Matrix" << std::endl;
2250 for (int n = 0; n < vdim; n++)
2251 {
2252 for (int l = 0; l < vdim; l++)
2253 app_log() << VT(l, n) << " ";
2254 app_log() << std::endl;
2255 }
2256 app_log() << std::endl << "Printing Eigenvalues" << std::endl;
2257 for (int n = 0; n < vdim; n++)
2258 app_log() << Sigma[n] << " ";
2259 app_log() << std::endl;
2260 }
2261 }
2262 }
2263
runNodePlot()2264 void WaveFunctionTester::runNodePlot()
2265 {
2266 app_log() << " ===== runNodePlot =====\n";
2267 xmlNodePtr kids = myNode->children;
2268 std::string doEnergy("no");
2269 ParameterSet aAttrib;
2270 aAttrib.add(doEnergy, "energy");
2271 aAttrib.put(myNode);
2272 std::vector<int> Grid;
2273 while (kids != NULL)
2274 {
2275 std::string cname((const char*)(kids->name));
2276 if (cname == "grid")
2277 putContent(Grid, kids);
2278 kids = kids->next;
2279 }
2280 ParticleSet::ParticlePos_t R_cart(1);
2281 R_cart.setUnit(PosUnit::Cartesian);
2282 ParticleSet::ParticlePos_t R_unit(1);
2283 R_unit.setUnit(PosUnit::Lattice);
2284 Walker_t& thisWalker(**(W.begin()));
2285 W.loadWalker(thisWalker, true);
2286 Walker_t::WFBuffer_t& w_buffer(thisWalker.DataSet);
2287 Psi.copyFromBuffer(W, w_buffer);
2288 #if OHMMS_DIM == 2
2289 assert(Grid.size() == 2);
2290 char fname[16];
2291 // sprintf(fname,"loc.xy");
2292 // std::ofstream e_out(fname);
2293 // e_out.precision(6);
2294 // e_out<<"#e x y"<< std::endl;
2295 int nat = W.getTotalNum();
2296 int nup = W.getTotalNum() / 2; //std::max(W.getSpeciesSet().findSpecies("u"),W.getSpeciesSet().findSpecies("d"));
2297 // for(int iat(0);iat<nat;iat++)
2298 // e_out<<iat<<" "<<W[0]->R[iat][0]<<" "<<W[0]->R[iat][1]<< std::endl;
2299 RealType overG0(1.0 / Grid[0]);
2300 RealType overG1(1.0 / Grid[1]);
2301 for (int iat(0); iat < nat; iat++)
2302 {
2303 W.update();
2304 std::stringstream fn;
2305 fn << RootName.c_str() << ".ratios." << iat << ".py";
2306 std::ofstream plot_out(fn.str().c_str());
2307 plot_out.precision(6);
2308 // plot_out<<"#e x y ratio"<< std::endl;
2309 R_unit[0][0] = 1.0;
2310 R_unit[0][1] = 1.0;
2311 W.convert2Cart(R_unit, R_cart);
2312 RealType xmax = R_cart[0][0];
2313 RealType ymax = R_cart[0][1];
2314 plot_out << "import matplotlib\n";
2315 plot_out << "import numpy as np\n";
2316 plot_out << "import matplotlib.cm as cm\n";
2317 plot_out << "import matplotlib.mlab as mlab\n";
2318 plot_out << "import matplotlib.pyplot as plt\n";
2319 plot_out << std::endl;
2320 plot_out << "matplotlib.rcParams['xtick.direction'] = 'out'\n";
2321 plot_out << "matplotlib.rcParams['ytick.direction'] = 'out'\n";
2322 plot_out << std::endl;
2323 plot_out << "x = np.arange(0, " << xmax << ", " << xmax * overG0 << ")\n";
2324 plot_out << "y = np.arange(0, " << ymax << ", " << ymax * overG1 << ")\n";
2325 plot_out << "X, Y = np.meshgrid(x, y)\n";
2326 plot_out << "Z = [";
2327 for (int i = 0; i < Grid[0]; i++)
2328 {
2329 plot_out << "[ ";
2330 for (int j = 0; j < Grid[1]; j++)
2331 {
2332 R_unit[0][0] = overG0 * RealType(i);
2333 R_unit[0][1] = overG1 * RealType(j);
2334 W.convert2Cart(R_unit, R_cart);
2335 PosType dr(R_cart[0] - W.R[iat]);
2336 W.makeMove(iat, dr);
2337 ValueType aratio = Psi.calcRatio(W, iat);
2338 W.rejectMove(iat);
2339 Psi.rejectMove(iat);
2340 plot_out << aratio << ", ";
2341 }
2342 plot_out << "], ";
2343 }
2344 plot_out << "]\n";
2345 plot_out << "up_y=[";
2346 for (int ix(0); ix < nup; ix++)
2347 {
2348 RealType yy(W[0]->R[ix][0]);
2349 while (yy > xmax)
2350 yy -= xmax;
2351 while (yy < 0)
2352 yy += xmax;
2353 plot_out << yy << ", ";
2354 }
2355 plot_out << "]\n";
2356 plot_out << "up_x=[";
2357 for (int ix(0); ix < nup; ix++)
2358 {
2359 RealType yy(W[0]->R[ix][1]);
2360 while (yy > ymax)
2361 yy -= ymax;
2362 while (yy < 0)
2363 yy += ymax;
2364 plot_out << yy << ", ";
2365 }
2366 plot_out << "]\n";
2367 plot_out << "dn_y=[";
2368 for (int ix(nup); ix < nat; ix++)
2369 {
2370 RealType yy(W[0]->R[ix][0]);
2371 while (yy > xmax)
2372 yy -= xmax;
2373 while (yy < 0)
2374 yy += xmax;
2375 plot_out << yy << ", ";
2376 }
2377 plot_out << "]\n";
2378 plot_out << "dn_x=[";
2379 for (int ix(nup); ix < nat; ix++)
2380 {
2381 RealType yy(W[0]->R[ix][1]);
2382 while (yy > ymax)
2383 yy -= ymax;
2384 while (yy < 0)
2385 yy += ymax;
2386 plot_out << yy << ", ";
2387 }
2388 plot_out << "]\n";
2389 plot_out << "matplotlib.rcParams['contour.negative_linestyle'] = 'solid'\n";
2390 plot_out << "plt.figure()\n";
2391 plot_out << "CS = plt.contourf(X, Y, Z, 5, cmap=cm.gray)\n";
2392 plot_out << "CS2 = plt.contour(X, Y, Z, colors='k',levels=[0])\n";
2393 plot_out << "PTu = plt.scatter(up_x,up_y, c='r', marker='o')\n";
2394 plot_out << "PTd = plt.scatter(dn_x,dn_y, c='b', marker='d')\n";
2395 plot_out << "plt.clabel(CS2, fontsize=9, inline=1)\n";
2396 plot_out << "plt.title('2D Nodal Structure')\n";
2397 plot_out << "plt.xlim(0," << ymax * (1.0 - overG1) << ")\n";
2398 plot_out << "plt.ylim(0," << xmax * (1.0 - overG0) << ")\n";
2399 fn.str("");
2400 fn << RootName.c_str() << ".ratios." << iat << ".png";
2401 plot_out << "plt.savefig('" << fn.str().c_str() << "', bbox_inches='tight', pad_inches=0.01 )\n";
2402 }
2403 #elif OHMMS_DIM == 3
2404 // assert(Grid.size()==3);
2405 //
2406 // RealType overG0(1.0/Grid[0]);
2407 // RealType overG1(1.0/Grid[1]);
2408 // RealType overG2(1.0/Grid[2]);
2409 // int iat(0);
2410 // W.update();
2411 // plot_out<<"#e x y z ratio"<< std::endl;
2412 //
2413 // for(int i=0;i<Grid[0];i++)
2414 // for(int j=0;j<Grid[1];j++)
2415 // for(int k=0;k<Grid[2];k++)
2416 // {
2417 // R_unit[iat][0]=overG0*RealType(i);
2418 // R_unit[iat][1]=overG1*RealType(j);
2419 // R_unit[iat][2]=overG2*RealType(k);
2420 // W.convert2Cart(R_unit,R_cart);
2421 // PosType dr(R_cart[iat]-W.R[iat]);
2422 //
2423 // W.makeMove(iat,dr);
2424 // RealType aratio = Psi.ratio(W,iat);
2425 // W.rejectMove(iat);
2426 // Psi.rejectMove(iat);
2427 // plot_out<<iat<<" "<<R_cart[iat][0]<<" "<<R_cart[iat][1]<<" "<<R_cart[iat][2]<<" "<<aratio<<" "<< std::endl;
2428 // }
2429 #endif
2430 }
2431
FiniteDiffErrData()2432 FiniteDiffErrData::FiniteDiffErrData() : particleIndex(0), gradientComponentIndex(0), outputFile("delta.dat") {}
2433
put(xmlNodePtr q)2434 bool FiniteDiffErrData::put(xmlNodePtr q)
2435 {
2436 ParameterSet param;
2437 param.add(outputFile, "file");
2438 param.add(particleIndex, "particle_index");
2439 param.add(gradientComponentIndex, "gradient_index");
2440 bool s = param.put(q);
2441 return s;
2442 }
2443
2444
2445 } // namespace qmcplusplus
2446