QMCPACK
test_vmc_driver.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) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Mark Dewing, markdewing@gmail.com, University of Illinois at Urbana-Champaign
8 //
9 // File created by: Mark Dewing, markdewing@gmail.com, University of Illinois at Urbana-Champaign
10 //////////////////////////////////////////////////////////////////////////////////////
11 
12 
13 #include "catch.hpp"
14 
15 
16 #include "Utilities/ProjectData.h"
17 #include "Utilities/RandomGenerator.h"
18 #include "OhmmsData/Libxml2Doc.h"
19 #include "OhmmsPETE/OhmmsMatrix.h"
20 #include "Particle/ParticleSet.h"
21 #include "Particle/MCWalkerConfiguration.h"
22 #include "Particle/ParticleSetPool.h"
23 #include "QMCHamiltonians/HamiltonianPool.h"
24 #include "QMCWaveFunctions/WaveFunctionPool.h"
25 #include "QMCWaveFunctions/WaveFunctionComponent.h"
26 #include "QMCWaveFunctions/TrialWaveFunction.h"
27 #include "QMCWaveFunctions/ConstantOrbital.h"
28 #include "QMCHamiltonians/BareKineticEnergy.h"
29 #include "Estimators/EstimatorManagerBase.h"
30 #include "Estimators/TraceManager.h"
31 #include "QMCDrivers/VMC/VMC.h"
32 
33 #include <stdio.h>
34 #include <string>
35 
36 
37 using std::string;
38 
39 namespace qmcplusplus
40 {
41 TEST_CASE("VMC", "[drivers][vmc]")
42 {
43  ProjectData project_data;
44  Communicate* c = OHMMS::Controller;
45  c->setName("test");
46  const SimulationCell simulation_cell;
47  ParticleSet ions(simulation_cell);
48  MCWalkerConfiguration elec(simulation_cell);
49 
50  ions.setName("ion");
51  ions.create({1});
52  ions.R[0] = {0.0, 0.0, 0.0};
53 
54  elec.setName("elec");
55  std::vector<int> agroup(1, 2);
56  elec.create(agroup);
57  elec.R[0] = {1.0, 0.0, 0.0};
58  elec.R[1] = {0.0, 0.0, 1.0};
59  elec.createWalkers(1);
60 
61  SpeciesSet& tspecies = elec.getSpeciesSet();
62  int upIdx = tspecies.addSpecies("u");
63  int chargeIdx = tspecies.addAttribute("charge");
64  int massIdx = tspecies.addAttribute("mass");
65  tspecies(chargeIdx, upIdx) = -1;
66  tspecies(massIdx, upIdx) = 1.0;
67 
68  elec.addTable(ions);
69  elec.update();
70 
71  CloneManager::clearClones();
72 
73  TrialWaveFunction psi(project_data.getRuntimeOptions());
74  psi.addComponent(std::make_unique<ConstantOrbital>());
75  psi.registerData(elec, elec[0]->DataSet);
76  elec[0]->DataSet.allocate();
77 
78  using RNG = RandomBase<QMCTraits::FullPrecRealType>;
79  UPtrVector<RNG> rngs(omp_get_max_threads());
80  for (std::unique_ptr<RNG>& rng : rngs)
81  rng = std::make_unique<FakeRandom<QMCTraits::FullPrecRealType>>();
82 
83  QMCHamiltonian h;
84  h.addOperator(std::make_unique<BareKineticEnergy>(elec, psi), "Kinetic");
85  h.addObservables(elec); // get double free error on 'h.Observables' w/o this
86 
87  elec.resetWalkerProperty(); // get memory corruption w/o this
88 
89  VMC vmc_omp(project_data, elec, psi, h, rngs, c, false);
90 
91  const char* vmc_input = R"(<qmc method="vmc" move="pbyp" checkpoint="-1">
92  <parameter name="substeps">1</parameter>
93  <parameter name="steps">1</parameter>
94  <parameter name="blocks">1</parameter>
95  <parameter name="timestep">0.1</parameter>
96  <parameter name="usedrift">no</parameter>
97  </qmc>
98  )";
99  Libxml2Document doc;
100  bool okay = doc.parseFromString(vmc_input);
101  REQUIRE(okay);
102  xmlNodePtr root = doc.getRoot();
103 
104  vmc_omp.process(root); // need to call 'process' for QMCDriver, which in turn calls 'put'
105 
106  vmc_omp.run();
107 
108  // With the constant wavefunction, no moves should be rejected
109  double ar = vmc_omp.acceptRatio();
110  CHECK(ar == Approx(1.0));
111 
112  // Each electron moved sqrt(tau)*gaussian_rng()
113  // See Particle>Base/tests/test_random_seq.cpp for the gaussian random numbers
114  // Values from diffuse.py for moving one step
115 
116  CHECK(elec[0]->R[0][0] == Approx(0.627670258894097));
117  CHECK(elec.R[0][1] == Approx(0.0));
118  CHECK(elec.R[0][2] == Approx(-0.372329741105903));
119 
120  CHECK(elec.R[1][0] == Approx(0.0));
121  CHECK(elec.R[1][1] == Approx(-0.372329741105903));
122  CHECK(elec.R[1][2] == Approx(1.0));
123 }
124 
125 TEST_CASE("SOVMC", "[drivers][vmc]")
126 {
127  ProjectData project_data;
128  Communicate* c = OHMMS::Controller;
129  c->setName("test");
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 
138  elec.setName("elec");
139  std::vector<int> agroup(1, 1);
140  elec.create(agroup);
141  elec.R[0] = {1.0, 0.0, 0.0};
142  elec.spins[0] = 0.0;
143  elec.setSpinor(true);
144  elec.createWalkers(1);
145 
146  SpeciesSet& tspecies = elec.getSpeciesSet();
147  int upIdx = tspecies.addSpecies("u");
148  int chargeIdx = tspecies.addAttribute("charge");
149  int massIdx = tspecies.addAttribute("mass");
150  tspecies(chargeIdx, upIdx) = -1;
151  tspecies(massIdx, upIdx) = 1.0;
152 
153  elec.addTable(ions);
154  elec.update();
155 
156  CloneManager::clearClones();
157 
158  TrialWaveFunction psi(project_data.getRuntimeOptions());
159  psi.addComponent(std::make_unique<ConstantOrbital>());
160  psi.registerData(elec, elec[0]->DataSet);
161  elec[0]->DataSet.allocate();
162 
163  using RNG = RandomBase<QMCTraits::FullPrecRealType>;
164  UPtrVector<RNG> rngs(omp_get_max_threads());
165  for (std::unique_ptr<RNG>& rng : rngs)
166  rng = std::make_unique<FakeRandom<QMCTraits::FullPrecRealType>>();
167 
168  QMCHamiltonian h;
169  h.addOperator(std::make_unique<BareKineticEnergy>(elec, psi), "Kinetic");
170  h.addObservables(elec); // get double free error on 'h.Observables' w/o this
171 
172  elec.resetWalkerProperty(); // get memory corruption w/o this
173 
174  VMC vmc_omp(project_data, elec, psi, h, rngs, c, false);
175 
176  const char* vmc_input = R"(<qmc method="vmc" move="pbyp" checkpoint="-1">
177  <parameter name="substeps">1</parameter>
178  <parameter name="steps">1</parameter>
179  <parameter name="blocks">1</parameter>
180  <parameter name="timestep">0.1</parameter>
181  <parameter name="usedrift">no</parameter>
182  <parameter name="SpinMass">0.25</parameter>
183  </qmc>
184  )";
185  Libxml2Document doc;
186  bool okay = doc.parseFromString(vmc_input);
187  REQUIRE(okay);
188  xmlNodePtr root = doc.getRoot();
189 
190  vmc_omp.process(root); // need to call 'process' for QMCDriver, which in turn calls 'put'
191 
192  vmc_omp.run();
193 
194  // With the constant wavefunction, no moves should be rejected
195  double ar = vmc_omp.acceptRatio();
196  CHECK(ar == Approx(1.0));
197 
198  // Each electron moved sqrt(tau)*gaussian_rng()
199  // See Particle>Base/tests/test_random_seq.cpp for the gaussian random numbers
200  // Values from diffuse.py for moving one step
201 
202  CHECK(elec.R[0][0] == Approx(0.627670258894097));
203  CHECK(elec.R[0][1] == Approx(0.0));
204  CHECK(elec.R[0][2] == Approx(-0.372329741105903));
205 
206  CHECK(elec.spins[0] == Approx(-0.74465948215809097));
207 
208  //Now we're going to test that the step updated the walker variables.
209  CHECK(elec[0]->R[0][0] == Approx(elec.R[0][0]));
210  CHECK(elec[0]->R[0][1] == Approx(elec.R[0][1]));
211  CHECK(elec[0]->R[0][2] == Approx(elec.R[0][2]));
212  CHECK(elec[0]->spins[0] == Approx(elec.spins[0]));
213 }
214 
215 TEST_CASE("SOVMC-alle", "[drivers][vmc]")
216 {
217  ProjectData project_data;
218  Communicate* c = OHMMS::Controller;
219  c->setName("test");
220  const SimulationCell simulation_cell;
221  ParticleSet ions(simulation_cell);
222  MCWalkerConfiguration elec(simulation_cell);
223 
224  ions.setName("ion");
225  ions.create({1});
226  ions.R[0] = {0.0, 0.0, 0.0};
227 
228  elec.setName("elec");
229  std::vector<int> agroup(1, 1);
230  elec.create(agroup);
231  elec.R[0] = {1.0, 0.0, 0.0};
232  elec.spins[0] = 0.0;
233  elec.setSpinor(true);
234  elec.createWalkers(1);
235 
236  SpeciesSet& tspecies = elec.getSpeciesSet();
237  int upIdx = tspecies.addSpecies("u");
238  int chargeIdx = tspecies.addAttribute("charge");
239  int massIdx = tspecies.addAttribute("mass");
240  tspecies(chargeIdx, upIdx) = -1;
241  tspecies(massIdx, upIdx) = 1.0;
242 
243  elec.addTable(ions);
244  elec.update();
245 
246  CloneManager::clearClones();
247 
248  TrialWaveFunction psi(project_data.getRuntimeOptions());
249  psi.addComponent(std::make_unique<ConstantOrbital>());
250  psi.registerData(elec, elec[0]->DataSet);
251  elec[0]->DataSet.allocate();
252 
253  using RNG = RandomBase<QMCTraits::FullPrecRealType>;
254  UPtrVector<RNG> rngs(omp_get_max_threads());
255  for (std::unique_ptr<RNG>& rng : rngs)
256  rng = std::make_unique<FakeRandom<QMCTraits::FullPrecRealType>>();
257 
258  QMCHamiltonian h;
259  h.addOperator(std::make_unique<BareKineticEnergy>(elec, psi), "Kinetic");
260  h.addObservables(elec); // get double free error on 'h.Observables' w/o this
261 
262  elec.resetWalkerProperty(); // get memory corruption w/o this
263 
264  VMC vmc_omp(project_data, elec, psi, h, rngs, c, false);
265 
266  const char* vmc_input = R"(<qmc method="vmc" move="alle" checkpoint="-1">
267  <parameter name="substeps">1</parameter>
268  <parameter name="steps">1</parameter>
269  <parameter name="blocks">1</parameter>
270  <parameter name="timestep">0.1</parameter>
271  <parameter name="usedrift">no</parameter>
272  <parameter name="SpinMass">0.25</parameter>
273  </qmc>
274  )";
275 
276  Libxml2Document doc;
277  bool okay = doc.parseFromString(vmc_input);
278  REQUIRE(okay);
279  xmlNodePtr root = doc.getRoot();
280 
281  vmc_omp.process(root); // need to call 'process' for QMCDriver, which in turn calls 'put'
282 
283  vmc_omp.run();
284 
285  // With the constant wavefunction, no moves should be rejected
286  double ar = vmc_omp.acceptRatio();
287  CHECK(ar == Approx(1.0));
288 
289  // Each electron moved sqrt(tau)*gaussian_rng()
290  // See Particle>Base/tests/test_random_seq.cpp for the gaussian random numbers
291  // Values from diffuse.py for moving one step
292 
293  CHECK(elec.R[0][0] == Approx(0.627670258894097));
294  CHECK(elec.R[0][1] == Approx(0.0));
295  CHECK(elec.R[0][2] == Approx(-0.372329741105903));
296 
297  CHECK(elec.spins[0] == Approx(-0.74465948215809097));
298 
299  //Now we're going to test that the step updated the walker variables.
300  CHECK(elec[0]->R[0][0] == Approx(elec.R[0][0]));
301  CHECK(elec[0]->R[0][1] == Approx(elec.R[0][1]));
302  CHECK(elec[0]->R[0][2] == Approx(elec.R[0][2]));
303  CHECK(elec[0]->spins[0] == Approx(elec.spins[0]));
304 }
305 } // namespace qmcplusplus
qmcplusplus::ParticleSet::setName
void setName(const std::string &aname)
Definition: ParticleSet.h:237
qmcplusplus::VMC::run
bool run() override
Definition: VMC.cpp:60
Libxml2Document
class that handles xmlDoc
Definition: Libxml2Doc.h:76
qmcplusplus::MCWalkerConfiguration
A set of walkers that are to be advanced by Metropolis Monte Carlo.
Definition: MCWalkerConfiguration.h:55
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::ParticleSet::spins
ParticleScalar spins
internal spin variables for dynamical spin calculations
Definition: ParticleSet.h:81
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::ProjectData
class ProjectData
Definition: ProjectData.h:36
qmcplusplus::TEST_CASE
TEST_CASE("complex_helper", "[type_traits]")
Definition: test_complex_helper.cpp:38
qmcplusplus::VMC
Implements a VMC using particle-by-particle move.
Definition: VMC.h:24
Libxml2Document::getRoot
xmlNodePtr getRoot()
Definition: Libxml2Doc.h:88
qmcplusplus::QMCHamiltonian
Collection of Local Energy Operators.
Definition: QMCHamiltonian.h:49
qmcplusplus::UPtrVector
std::vector< std::unique_ptr< T > > UPtrVector
Definition: template_types.hpp:38
qmcplusplus::ParticleSet::update
void update(bool skipSK=false)
update the internal data
Definition: ParticleSet.cpp:349
OHMMS::Controller
Communicate * Controller
Global Communicator for a process.
Definition: Communicate.cpp:35
qmcplusplus::SpeciesSet::addAttribute
int addAttribute(const std::string &aname)
for a new attribute, allocate the data, !More often used to get the index of a species ...
Definition: SpeciesSet.cpp:45
qmcplusplus::ProjectData::getRuntimeOptions
const RuntimeOptions & getRuntimeOptions() const noexcept
Definition: ProjectData.cpp:265
Communicate
Wrapping information on parallelism.
Definition: Communicate.h:68
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
omp_get_max_threads
omp_int_t omp_get_max_threads()
Definition: OpenMP.h:26
qmcplusplus::REQUIRE
REQUIRE(std::filesystem::exists(filename))
WaveFunctionPool.h
Declaration of WaveFunctionPool.
qmcplusplus::ParticleSet::R
ParticlePos R
Position.
Definition: ParticleSet.h:79
EstimatorManagerBase.h
Manager class of scalar estimators.
Communicate::setName
void setName(const std::string &aname)
Definition: Communicate.h:129
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
TrialWaveFunction.h
Declaration of a TrialWaveFunction.
HamiltonianPool.h
Declaration of HamiltonianPool.
qmcplusplus::TrialWaveFunction
Class to represent a many-body trial wave function.
Definition: TrialWaveFunction.h:60
Libxml2Document::parseFromString
bool parseFromString(const std::string_view data)
Definition: Libxml2Doc.cpp:204
qmcplusplus::CHECK
CHECK(log_values[0]==ComplexApprox(std::complex< double >{ 5.603777579195571, -6.1586603331188225 }))
qmcplusplus::CloneManager::acceptRatio
RealType acceptRatio() const
Definition: CloneManager.cpp:247
WaveFunctionComponent.h
Declaration of WaveFunctionComponent.
qmcplusplus::SpeciesSet
Custom container for set of attributes for a set of species.
Definition: SpeciesSet.h:33
qmcplusplus::ParticleSet::setSpinor
void setSpinor(bool is_spinor)
Definition: ParticleSet.h:257
qmcplusplus::MCWalkerConfiguration::createWalkers
void createWalkers(int numWalkers)
create numWalkers Walkers
Definition: MCWalkerConfiguration.cpp:50
RandomGenerator.h
Declare a global Random Number Generator.
MCWalkerConfiguration.h
Declaration of a MCWalkerConfiguration.
qmcplusplus::QMCDriver::process
void process(xmlNodePtr cur) override
initialize with xmlNode
Definition: QMCDriver.cpp:172
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