QMCPACK
test_dmc.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) 2017 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory
8 //
9 // File created by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory
10 //////////////////////////////////////////////////////////////////////////////////////
11 
12 
13 #include "catch.hpp"
14 
15 #include "Utilities/RandomGenerator.h"
16 #include "OhmmsData/Libxml2Doc.h"
17 #include "OhmmsPETE/OhmmsMatrix.h"
18 #include "Particle/ParticleSet.h"
19 #include "Particle/MCWalkerConfiguration.h"
20 #include "Particle/ParticleSetPool.h"
21 #include "QMCWaveFunctions/WaveFunctionComponent.h"
22 #include "QMCWaveFunctions/TrialWaveFunction.h"
23 #include "QMCWaveFunctions/ConstantOrbital.h"
24 #include "QMCWaveFunctions/LinearOrbital.h"
25 #include "QMCHamiltonians/BareKineticEnergy.h"
26 #include "Estimators/EstimatorManagerBase.h"
27 #include "Estimators/TraceManager.h"
28 #include "QMCDrivers/DMC/DMCUpdatePbyP.h"
29 #include "QMCDrivers/GreenFunctionModifiers/DriftModifierUNR.h"
30 #include "Utilities/RuntimeOptions.h"
31 
32 
33 #include <stdio.h>
34 #include <string>
35 
36 
37 using std::string;
38 
39 namespace qmcplusplus
40 {
41 TEST_CASE("DMC Particle-by-Particle advanceWalkers ConstantOrbital", "[drivers][dmc]")
42 {
43  const SimulationCell simulation_cell;
44  ParticleSet ions(simulation_cell);
45  MCWalkerConfiguration elec(simulation_cell);
46 
47  ions.setName("ion");
48  ions.create({1});
49  ions.R[0] = {0.0, 0.0, 0.0};
50  elec.setName("elec");
51  std::vector<int> agroup(1);
52  agroup[0] = 2;
53  elec.create(agroup);
54  elec.R[0] = {1.0, 0.0, 0.0};
55  elec.R[1] = {0.0, 0.0, 1.0};
56  elec.createWalkers(1);
57 
58  SpeciesSet& tspecies = elec.getSpeciesSet();
59  int upIdx = tspecies.addSpecies("u");
60  int chargeIdx = tspecies.addAttribute("charge");
61  int massIdx = tspecies.addAttribute("mass");
62  tspecies(chargeIdx, upIdx) = -1;
63  tspecies(massIdx, upIdx) = 1.0;
64 
65  elec.addTable(ions);
66  elec.update();
67 
68  RuntimeOptions runtime_options;
69  TrialWaveFunction psi(runtime_options);
70  auto orb_uptr = std::make_unique<ConstantOrbital>();
71  auto orb = orb_uptr.get();
72  psi.addComponent(std::move(orb_uptr));
73  psi.registerData(elec, elec[0]->DataSet);
74  elec[0]->DataSet.allocate();
75 
76  FakeRandom rg;
77 
78  QMCHamiltonian h;
79  h.addOperator(std::make_unique<BareKineticEnergy>(elec, psi), "Kinetic");
80  h.addObservables(elec); // get double free error on 'h.Observables' w/o this
81 
82  elec.resetWalkerProperty(); // get memory corruption w/o this
83 
84  DMCUpdatePbyPWithRejectionFast dmc(elec, psi, h, rg);
85  EstimatorManagerBase EM;
86  double tau = 0.1;
87  SimpleFixedNodeBranch branch(tau, 1);
88  TraceManager TM;
89  DriftModifierUNR DM;
90  dmc.resetRun(&branch, &EM, &TM, &DM);
91  dmc.startBlock(1);
92 
93  DMCUpdatePbyPWithRejectionFast::WalkerIter_t begin = elec.begin();
94  DMCUpdatePbyPWithRejectionFast::WalkerIter_t end = elec.end();
95  dmc.advanceWalkers(begin, end, true);
96 
97  // With the constant wavefunction, no moves should be rejected
98  REQUIRE(dmc.nReject == 0);
99  REQUIRE(dmc.nAccept == 2);
100 
101  // Each electron moved sqrt(tau)*gaussian_rng()
102  // See ParticleBase/tests/test_random_seq.cpp for the gaussian random numbers
103  CHECK(elec.R[0][0] == Approx(0.6276702589209545));
104  CHECK(elec.R[0][1] == Approx(0.0));
105  CHECK(elec.R[0][2] == Approx(-0.3723297410790455));
106 
107  CHECK(elec.R[1][0] == Approx(0.0));
108  CHECK(elec.R[1][1] == Approx(-0.3723297410790455));
109  CHECK(elec.R[1][2] == Approx(1.0));
110 
111 
112  // Check rejection in case of node-crossing
113 
114  orb->FakeGradRatio = -1.0;
115 
116  dmc.advanceWalkers(begin, end, true);
117 #ifdef QMC_COMPLEX
118  REQUIRE(dmc.nReject == 0);
119  REQUIRE(dmc.nAccept == 4);
120 #else
121  // Should be rejected
122  REQUIRE(dmc.nReject == 2);
123  REQUIRE(dmc.nAccept == 2);
124 #endif
125 }
126 
127 TEST_CASE("DMC Particle-by-Particle advanceWalkers LinearOrbital", "[drivers][dmc]")
128 
129 {
130  const SimulationCell simulation_cell;
131  ParticleSet ions(simulation_cell);
132  MCWalkerConfiguration elec(simulation_cell);
133 
134  ions.setName("ion");
135  ions.create({1});
136  ions.R[0] = {0.0, 0.0, 0.0};
137  elec.setName("elec");
138  std::vector<int> agroup(1);
139  agroup[0] = 2;
140  elec.create(agroup);
141  elec.R[0] = {1.0, 0.0, 0.0};
142  elec.R[1] = {0.0, 0.0, 1.0};
143  elec.createWalkers(1);
144 
145  SpeciesSet& tspecies = elec.getSpeciesSet();
146  int upIdx = tspecies.addSpecies("u");
147  int chargeIdx = tspecies.addAttribute("charge");
148  int massIdx = tspecies.addAttribute("mass");
149  tspecies(chargeIdx, upIdx) = -1;
150  tspecies(massIdx, upIdx) = 1.0;
151 
152  elec.addTable(ions);
153  elec.update();
154 
155  RuntimeOptions runtime_options;
156  TrialWaveFunction psi(runtime_options);
157  psi.addComponent(std::make_unique<LinearOrbital>());
158  psi.registerData(elec, elec[0]->DataSet);
159  elec[0]->DataSet.allocate();
160 
161  FakeRandom rg;
162 
163  QMCHamiltonian h;
164  h.addOperator(std::make_unique<BareKineticEnergy>(elec, psi), "Kinetic");
165  h.addObservables(elec); // get double free error on 'h.Observables' w/o this
166 
167  elec.resetWalkerProperty(); // get memory corruption w/o this
168 
169  DMCUpdatePbyPWithRejectionFast dmc(elec, psi, h, rg);
170  EstimatorManagerBase EM;
171  double tau = 0.1;
172  SimpleFixedNodeBranch branch(tau, 1);
173  TraceManager TM;
174  DriftModifierUNR DM;
175  dmc.resetRun(&branch, &EM, &TM, &DM);
176  dmc.startBlock(1);
177 
178  DMCUpdatePbyPWithRejectionFast::WalkerIter_t begin = elec.begin();
179  DMCUpdatePbyPWithRejectionFast::WalkerIter_t end = elec.end();
180  dmc.advanceWalkers(begin, end, true);
181 
182  REQUIRE(dmc.nReject == 0);
183  REQUIRE(dmc.nAccept == 2);
184 
185  // Each electron moved sqrt(tau)*gaussian_rng()
186  // See ParticleBase/tests/test_random_seq.cpp for the gaussian random numbers
187  // See DMC_propagator notebook for computation of these values
188  CHECK(elec.R[0][0] == Approx(0.695481606677082));
189  CHECK(elec.R[0][1] == Approx(0.135622695565971));
190  CHECK(elec.R[0][2] == Approx(-0.168895697756948));
191 
192  CHECK(elec.R[1][0] == Approx(0.0678113477829853));
193  CHECK(elec.R[1][1] == Approx(-0.236707045539933));
194  CHECK(elec.R[1][2] == Approx(1.20343404334896));
195 }
196 } // namespace qmcplusplus
qmcplusplus::QMCHamiltonian::addObservables
int addObservables(ParticleSet &P)
add each term to the PropertyList for averages
Definition: QMCHamiltonian.cpp:223
qmcplusplus::ParticleSet::setName
void setName(const std::string &aname)
Definition: ParticleSet.h:237
qmcplusplus::MCWalkerConfiguration
A set of walkers that are to be advanced by Metropolis Monte Carlo.
Definition: MCWalkerConfiguration.h:55
qmcplusplus::QMCUpdateBase::WalkerIter_t
MCWalkerConfiguration::iterator WalkerIter_t
Definition: QMCUpdateBase.h:45
qmcplusplus::SpeciesSet::addSpecies
int addSpecies(const std::string &aname)
When a name species does not exist, add a new species.
Definition: SpeciesSet.cpp:33
qmcplusplus
helper functions for EinsplineSetBuilder
Definition: Configuration.h:43
qmcplusplus::QMCHamiltonian::addOperator
void addOperator(std::unique_ptr< OperatorBase > &&h, const std::string &aname, bool physical=true)
add an operator
Definition: QMCHamiltonian.cpp:99
qmcplusplus::TEST_CASE
TEST_CASE("complex_helper", "[type_traits]")
Definition: test_complex_helper.cpp:38
qmcplusplus::QMCHamiltonian
Collection of Local Energy Operators.
Definition: QMCHamiltonian.h:49
qmcplusplus::ParticleSet::update
void update(bool skipSK=false)
update the internal data
Definition: ParticleSet.cpp:349
qmcplusplus::SpeciesSet::addAttribute
int addAttribute(const std::string &aname)
for a new attribute, allocate the data, !More often used to get the index of a species ...
Definition: SpeciesSet.cpp:45
qmcplusplus::QMCUpdateBase::startBlock
void startBlock(int steps)
prepare to start a block
Definition: QMCUpdateBase.cpp:194
qmcplusplus::QMCUpdateBase::resetRun
void resetRun(BranchEngineType *brancher, EstimatorManagerBase *est, TraceManager *traces, const DriftModifierBase *driftmodifer)
reset the QMCUpdateBase parameters
Definition: QMCUpdateBase.cpp:130
qmcplusplus::ParticleSet::addTable
int addTable(const ParticleSet &psrc, DTModes modes=DTModes::ALL_OFF)
add a distance table
Definition: ParticleSet.cpp:307
qmcplusplus::ParticleSet
Specialized paritlce class for atomistic simulations.
Definition: ParticleSet.h:55
qmcplusplus::SimpleFixedNodeBranch
Manages the state of QMC sections and handles population control for DMCs.
Definition: SimpleFixedNodeBranch.h:99
qmcplusplus::REQUIRE
REQUIRE(std::filesystem::exists(filename))
qmcplusplus::TrialWaveFunction::registerData
void registerData(ParticleSet &P, WFBufferType &buf)
register all the wavefunction components in buffer.
Definition: TrialWaveFunction.cpp:958
qmcplusplus::ParticleSet::R
ParticlePos R
Position.
Definition: ParticleSet.h:79
EstimatorManagerBase.h
Manager class of scalar estimators.
qmcplusplus::ParticleSet::getSpeciesSet
SpeciesSet & getSpeciesSet()
retrun the SpeciesSet of this particle set
Definition: ParticleSet.h:231
qmcplusplus::ParticleSet::create
void create(const std::vector< int > &agroup)
create grouped particles
Definition: ParticleSet.cpp:125
qmcplusplus::EstimatorManagerBase
Class to manage a set of ScalarEstimators.
Definition: EstimatorManagerBase.h:47
TrialWaveFunction.h
Declaration of a TrialWaveFunction.
qmcplusplus::TrialWaveFunction
Class to represent a many-body trial wave function.
Definition: TrialWaveFunction.h:60
qmcplusplus::CHECK
CHECK(log_values[0]==ComplexApprox(std::complex< double >{ 5.603777579195571, -6.1586603331188225 }))
WaveFunctionComponent.h
Declaration of WaveFunctionComponent.
qmcplusplus::SpeciesSet
Custom container for set of attributes for a set of species.
Definition: SpeciesSet.h:33
qmcplusplus::QMCUpdateBase::nAccept
IndexType nAccept
counter for number of moves accepted
Definition: QMCUpdateBase.h:63
qmcplusplus::QMCUpdateBase::nReject
IndexType nReject
counter for number of moves rejected
Definition: QMCUpdateBase.h:65
qmcplusplus::MCWalkerConfiguration::createWalkers
void createWalkers(int numWalkers)
create numWalkers Walkers
Definition: MCWalkerConfiguration.cpp:50
qmcplusplus::QMCUpdateBase::advanceWalkers
virtual void advanceWalkers(WalkerIter_t it, WalkerIter_t it_end, bool recompute)
advance walkers executed at each step
Definition: QMCUpdateBase.cpp:359
qmcplusplus::WalkerConfigurations::end
iterator end()
return the last iterator, [begin(), end())
Definition: WalkerConfigurations.h:123
RandomGenerator.h
Declare a global Random Number Generator.
MCWalkerConfiguration.h
Declaration of a MCWalkerConfiguration.
qmcplusplus::TrialWaveFunction::addComponent
void addComponent(std::unique_ptr< WaveFunctionComponent > &&aterm)
add a WaveFunctionComponent
Definition: TrialWaveFunction.cpp:87
ParticleSetPool.h
Declaration of ParticleSetPool.
qmcplusplus::MCWalkerConfiguration::resetWalkerProperty
void resetWalkerProperty(int ncopy=1)
reset the Property container of all the walkers
Definition: MCWalkerConfiguration.cpp:97
qmcplusplus::WalkerConfigurations::begin
iterator begin()
return the first iterator
Definition: WalkerConfigurations.h:121