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) 2018 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 16 #include "Utilities/RandomGenerator.h" 17 #include "OhmmsData/Libxml2Doc.h" 18 #include "OhmmsPETE/OhmmsMatrix.h" 19 #include "Particle/ParticleSet.h" 20 #include "Particle/MCWalkerConfiguration.h" 21 #include "Particle/ParticleSetPool.h" 22 #include "QMCHamiltonians/HamiltonianPool.h" 23 #include "QMCWaveFunctions/WaveFunctionPool.h" 24 #include "QMCWaveFunctions/WaveFunctionComponent.h" 25 #include "QMCWaveFunctions/TrialWaveFunction.h" 26 #include "QMCWaveFunctions/ConstantOrbital.h" 27 #include "QMCHamiltonians/BareKineticEnergy.h" 28 #include "Estimators/EstimatorManagerBase.h" 29 #include "Estimators/TraceManager.h" 30 #include "QMCDrivers/DMC/DMC.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", "[drivers][dmc]") 42 { 43 Communicate* c; 44 c = OHMMS::Controller; 45 46 ParticleSet ions; 47 MCWalkerConfiguration elec; 48 49 ions.setName("ion"); 50 ions.create(1); 51 ions.R[0][0] = 0.0; 52 ions.R[0][1] = 0.0; 53 ions.R[0][2] = 0.0; 54 55 elec.setName("elec"); 56 std::vector<int> agroup(1); 57 agroup[0] = 2; 58 elec.create(agroup); 59 elec.R[0][0] = 1.0; 60 elec.R[0][1] = 0.0; 61 elec.R[0][2] = 0.0; 62 elec.R[1][0] = 0.0; 63 elec.R[1][1] = 0.0; 64 elec.R[1][2] = 1.0; 65 elec.createWalkers(1); 66 67 SpeciesSet& tspecies = elec.getSpeciesSet(); 68 int upIdx = tspecies.addSpecies("u"); 69 int chargeIdx = tspecies.addAttribute("charge"); 70 int massIdx = tspecies.addAttribute("mass"); 71 tspecies(chargeIdx, upIdx) = -1; 72 tspecies(massIdx, upIdx) = 1.0; 73 74 elec.addTable(ions); 75 elec.update(); 76 77 CloneManager::clear_for_unit_tests(); 78 79 TrialWaveFunction psi; 80 ConstantOrbital* orb = new ConstantOrbital; 81 psi.addComponent(orb); 82 psi.registerData(elec, elec.WalkerList[0]->DataSet); 83 elec.WalkerList[0]->DataSet.allocate(); 84 85 FakeRandom rg; 86 87 QMCHamiltonian h; 88 BareKineticEnergy<double>* p_bke = new BareKineticEnergy<double>(elec); 89 h.addOperator(p_bke, "Kinetic"); 90 h.addObservables(elec); // get double free error on 'h.Observables' w/o this 91 92 elec.resetWalkerProperty(); // get memory corruption w/o this 93 94 DMC dmc_omp(elec, psi, h, c, false); 95 96 const char* dmc_input = "<qmc method=\"dmc\"> \ 97 <parameter name=\"steps\">1</parameter> \ 98 <parameter name=\"blocks\">1</parameter> \ 99 <parameter name=\"timestep\">0.1</parameter> \ 100 </qmc> \ 101 "; 102 Libxml2Document* doc = new Libxml2Document; 103 bool okay = doc->parseFromString(dmc_input); 104 REQUIRE(okay); 105 xmlNodePtr root = doc->getRoot(); 106 107 dmc_omp.process(root); // need to call 'process' for QMCDriver, which in turn calls 'put' 108 109 dmc_omp.run(); 110 111 // With the constant wavefunction, no moves should be rejected 112 double ar = dmc_omp.acceptRatio(); 113 REQUIRE(ar == Approx(1.0)); 114 115 // Each electron moved sqrt(tau)*gaussian_rng() 116 // See Particle>Base/tests/test_random_seq.cpp for the gaussian random numbers 117 // Values from diffuse.py for moving one step 118 119 REQUIRE(elec[0]->R[0][0] == Approx(0.627670258894097)); 120 REQUIRE(elec.R[0][1] == Approx(0.0)); 121 REQUIRE(elec.R[0][2] == Approx(-0.372329741105903)); 122 123 REQUIRE(elec.R[1][0] == Approx(0.0)); 124 REQUIRE(elec.R[1][1] == Approx(-0.372329741105903)); 125 REQUIRE(elec.R[1][2] == Approx(1.0)); 126 127 delete doc; 128 delete p_bke; 129 } 130 131 TEST_CASE("SODMC", "[drivers][dmc]") 132 { 133 Communicate* c; 134 c = OHMMS::Controller; 135 136 ParticleSet ions; 137 MCWalkerConfiguration elec; 138 139 ions.setName("ion"); 140 ions.create(1); 141 ions.R[0][0] = 0.0; 142 ions.R[0][1] = 0.0; 143 ions.R[0][2] = 0.0; 144 145 elec.setName("elec"); 146 std::vector<int> agroup(1); 147 agroup[0] = 1; 148 elec.create(agroup); 149 elec.R[0][0] = 1.0; 150 elec.R[0][1] = 0.0; 151 elec.R[0][2] = 0.0; 152 elec.spins[0] = 0.0; 153 elec.createWalkers(1); 154 155 SpeciesSet& tspecies = elec.getSpeciesSet(); 156 int upIdx = tspecies.addSpecies("u"); 157 int chargeIdx = tspecies.addAttribute("charge"); 158 int massIdx = tspecies.addAttribute("mass"); 159 tspecies(chargeIdx, upIdx) = -1; 160 tspecies(massIdx, upIdx) = 1.0; 161 162 elec.addTable(ions); 163 elec.update(); 164 165 CloneManager::clear_for_unit_tests(); 166 167 TrialWaveFunction psi; 168 ConstantOrbital* orb = new ConstantOrbital; 169 psi.addComponent(orb); 170 psi.registerData(elec, elec.WalkerList[0]->DataSet); 171 elec.WalkerList[0]->DataSet.allocate(); 172 173 FakeRandom rg; 174 175 QMCHamiltonian h; 176 BareKineticEnergy<double>* p_bke = new BareKineticEnergy<double>(elec); 177 h.addOperator(p_bke, "Kinetic"); 178 h.addObservables(elec); // get double free error on 'h.Observables' w/o this 179 180 elec.resetWalkerProperty(); // get memory corruption w/o this 181 182 DMC dmc_omp(elec, psi, h, c, false); 183 184 const char* dmc_input = "<qmc method=\"dmc\"> \ 185 <parameter name=\"steps\">1</parameter> \ 186 <parameter name=\"blocks\">1</parameter> \ 187 <parameter name=\"timestep\">0.1</parameter> \ 188 <parameter name=\"SpinMoves\">yes</parameter> \ 189 <parameter name=\"SpinMass\">0.25</parameter> \ 190 </qmc> \ 191 "; 192 Libxml2Document* doc = new Libxml2Document; 193 bool okay = doc->parseFromString(dmc_input); 194 REQUIRE(okay); 195 xmlNodePtr root = doc->getRoot(); 196 197 dmc_omp.process(root); // need to call 'process' for QMCDriver, which in turn calls 'put' 198 199 dmc_omp.run(); 200 201 // With the constant wavefunction, no moves should be rejected 202 double ar = dmc_omp.acceptRatio(); 203 REQUIRE(ar == Approx(1.0)); 204 205 // Each electron moved sqrt(tau)*gaussian_rng() 206 // See Particle>Base/tests/test_random_seq.cpp for the gaussian random numbers 207 // Values from diffuse.py for moving one step 208 209 REQUIRE(elec[0]->R[0][0] == Approx(0.627670258894097)); 210 REQUIRE(elec.R[0][1] == Approx(0.0)); 211 REQUIRE(elec.R[0][2] == Approx(-0.372329741105903)); 212 213 REQUIRE(elec.spins[0] == Approx(-0.74465948215809097)); 214 215 delete doc; 216 delete p_bke; 217 } 218 } // namespace qmcplusplus 219