QMCPACK
test_force_ewald.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 #include "OhmmsData/Libxml2Doc.h"
16 #include "OhmmsPETE/OhmmsMatrix.h"
17 #include "Particle/ParticleSet.h"
18 #include "Particle/ParticleSetPool.h"
19 #include "QMCHamiltonians/ForceChiesaPBCAA.h"
20 #include "QMCHamiltonians/ForceCeperley.h"
21 #include "QMCWaveFunctions/TrialWaveFunction.h"
22 #include "LongRange/EwaldHandler3D.h"
23 #include "Utilities/RuntimeOptions.h"
24 
25 #include <stdio.h>
26 #include <string>
27 
28 using std::string;
29 
30 namespace qmcplusplus
31 {
32 // PBC case
33 TEST_CASE("Chiesa Force BCC H Ewald3D", "[hamiltonian]")
34 {
35  Communicate* c = OHMMS::Controller;
36 
37  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
38  lattice.BoxBConds = true; // periodic
39  lattice.R.diagonal(3.77945227);
40  lattice.LR_dim_cutoff = 40;
41  lattice.reset();
42 
43  const SimulationCell simulation_cell(lattice);
44  ParticleSet ions(simulation_cell);
45  ParticleSet elec(simulation_cell);
46 
47  ions.setName("ion");
48  ions.create({2});
49  ions.R[0] = {0.0, 0.0, 0.0};
50  ions.R[1] = {1.6, 1.6, 1.88972614};
51  SpeciesSet& ion_species = ions.getSpeciesSet();
52  int pIdx = ion_species.addSpecies("H");
53  int pChargeIdx = ion_species.addAttribute("charge");
54  ion_species(pChargeIdx, pIdx) = 1;
55  ions.createSK();
56 
57 
58  elec.setName("elec");
59  elec.create({2});
60  elec.R[0] = {0.5, 0.0, 0.0};
61  elec.R[1] = {0.0, 0.5, 0.0};
62  SpeciesSet& tspecies = elec.getSpeciesSet();
63  int upIdx = tspecies.addSpecies("u");
64  int chargeIdx = tspecies.addAttribute("charge");
65  int massIdx = tspecies.addAttribute("mass");
66  tspecies(chargeIdx, upIdx) = -1;
67  tspecies(massIdx, upIdx) = 1.0;
68 
69  elec.createSK();
70 
71  ParticleSetPool ptcl = ParticleSetPool(c);
72 
73  ions.resetGroups();
74 
75  // The call to resetGroups is needed transfer the SpeciesSet
76  // settings to the ParticleSet
77  elec.resetGroups();
78 
79  LRCoulombSingleton::CoulombHandler = std::make_unique<EwaldHandler3D>(ions);
80  LRCoulombSingleton::CoulombHandler->initBreakup(ions);
81  LRCoulombSingleton::CoulombDerivHandler = std::make_unique<EwaldHandler3D>(ions);
82  LRCoulombSingleton::CoulombDerivHandler->initBreakup(ions);
83 
84  ForceChiesaPBCAA force(ions, elec);
85  force.setAddIonIon(false);
86  force.InitMatrix();
87 
88  elec.update();
89  force.evaluate(elec);
90 
91  //Ion-Ion forces are validated against Quantum Espresso's ewald method:
92  CHECK(force.getForcesIonIon()[0][0] == Approx(-0.0228366));
93  CHECK(force.getForcesIonIon()[0][1] == Approx(-0.0228366));
94  CHECK(force.getForcesIonIon()[0][2] == Approx(0.0000000));
95  CHECK(force.getForcesIonIon()[1][0] == Approx(0.0228366));
96  CHECK(force.getForcesIonIon()[1][1] == Approx(0.0228366));
97  CHECK(force.getForcesIonIon()[1][2] == Approx(0.0000000));
98 
99  //Electron-Ion forces are unvalidated externally:
100  CHECK(force.getForces()[0][0] == Approx(3.959178977));
101  CHECK(force.getForces()[0][1] == Approx(3.959178977));
102  CHECK(force.getForces()[0][2] == Approx(0.000000000));
103  CHECK(force.getForces()[1][0] == Approx(-0.078308730));
104  CHECK(force.getForces()[1][1] == Approx(-0.078308730));
105  CHECK(force.getForces()[1][2] == Approx(0.000000000));
106 
107  LRCoulombSingleton::CoulombHandler.reset(nullptr);
108  LRCoulombSingleton::CoulombDerivHandler.reset(nullptr);
109 }
110 
111 // test SR and LR pieces separately
112 TEST_CASE("fccz sr lr clone", "[hamiltonian]")
113 {
114  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
115  lattice.BoxBConds = true; // periodic
116  lattice.R.diagonal(3.77945227);
117  lattice.LR_dim_cutoff = 40;
118  lattice.reset();
119 
120  const SimulationCell simulation_cell(lattice);
121  ParticleSet ions(simulation_cell);
122  ParticleSet elec(simulation_cell);
123 
124  ions.setName("ion");
125  ions.create({2});
126  ions.R[0] = {0.0, 0.0, 0.0};
127  ions.R[1] = {1.6, 1.6, 1.88972614};
128  SpeciesSet& ion_species = ions.getSpeciesSet();
129  int pIdx = ion_species.addSpecies("H");
130  int pChargeIdx = ion_species.addAttribute("charge");
131  ion_species(pChargeIdx, pIdx) = 1;
132  ions.createSK();
133 
134 
135  elec.setName("elec");
136  elec.create({2});
137  elec.R[0] = {0.5, 0.0, 0.0};
138  elec.R[1] = {0.0, 0.5, 0.0};
139  SpeciesSet& tspecies = elec.getSpeciesSet();
140  int upIdx = tspecies.addSpecies("u");
141  int chargeIdx = tspecies.addAttribute("charge");
142  int massIdx = tspecies.addAttribute("mass");
143  tspecies(chargeIdx, upIdx) = -1;
144  tspecies(massIdx, upIdx) = 1.0;
145  elec.createSK();
146 
147  // The call to resetGroups is needed transfer the SpeciesSet
148  // settings to the ParticleSet
149  ions.resetGroups();
150  elec.resetGroups();
151 
152  LRCoulombSingleton::CoulombHandler = std::make_unique<EwaldHandler3D>(ions);
153  LRCoulombSingleton::CoulombHandler->initBreakup(ions);
154  LRCoulombSingleton::CoulombDerivHandler = std::make_unique<EwaldHandler3D>(ions);
155  LRCoulombSingleton::CoulombDerivHandler->initBreakup(ions);
156 
157  ForceChiesaPBCAA force(ions, elec);
158  force.setAddIonIon(false);
159  force.InitMatrix();
160 
161  elec.update();
162  // test SR part
163  force.evaluateSR(elec);
164  CHECK(force.getForces()[0][0] == Approx(1.6938118975));
165  CHECK(force.getForces()[0][1] == Approx(1.6938118975));
166  CHECK(force.getForces()[0][2] == Approx(0.000000000));
167  CHECK(force.getForces()[1][0] == Approx(0.000000000));
168  CHECK(force.getForces()[1][1] == Approx(0.000000000));
169  CHECK(force.getForces()[1][2] == Approx(0.000000000));
170  // test LR part
171  force.setForces(0);
172  force.evaluateLR(elec);
173  CHECK(force.getForces()[0][0] == Approx(2.2653670795));
174  CHECK(force.getForces()[0][1] == Approx(2.2653670795));
175  CHECK(force.getForces()[0][2] == Approx(0.000000000));
176  CHECK(force.getForces()[1][0] == Approx(-0.078308730));
177  CHECK(force.getForces()[1][1] == Approx(-0.078308730));
178  CHECK(force.getForces()[1][2] == Approx(0.000000000));
179 
180  // test cloning !!!! makeClone is not testable
181  // example call path:
182  // QMCDrivers/CloneManager::makeClones
183  // QMCHamiltonian::makeClone
184  // OperatorBase::add2Hamiltonian -> ForceChiesaPBCAA::makeClone
185  RuntimeOptions runtime_options;
186  TrialWaveFunction psi(runtime_options);
187  std::unique_ptr<ForceChiesaPBCAA> clone(dynamic_cast<ForceChiesaPBCAA*>(force.makeClone(elec, psi).release()));
188  clone->evaluate(elec);
189  REQUIRE(clone->getAddIonIon() == force.getAddIonIon());
190  CHECK(clone->getForcesIonIon()[0][0] == Approx(-0.0228366));
191  CHECK(clone->getForcesIonIon()[0][1] == Approx(-0.0228366));
192  CHECK(clone->getForcesIonIon()[0][2] == Approx(0.0000000));
193  CHECK(clone->getForcesIonIon()[1][0] == Approx(0.0228366));
194  CHECK(clone->getForcesIonIon()[1][1] == Approx(0.0228366));
195  CHECK(clone->getForcesIonIon()[1][2] == Approx(0.0000000));
196 
197  LRCoulombSingleton::CoulombHandler.reset(nullptr);
198  LRCoulombSingleton::CoulombDerivHandler.reset(nullptr);
199 }
200 
201 // 3 H atoms randomly distributed in a box
202 TEST_CASE("fccz h3", "[hamiltonian]")
203 {
204  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
205  lattice.BoxBConds = true; // periodic
206  lattice.R.diagonal(3.77945227);
207  lattice.LR_dim_cutoff = 40;
208  lattice.reset();
209 
210  const SimulationCell simulation_cell(lattice);
211  ParticleSet ions(simulation_cell);
212  ParticleSet elec(simulation_cell);
213 
214  ions.setName("ion");
215  ions.create({3});
216  ions.R[0] = {0.0, 0.0, 0.0};
217  ions.R[1] = {1.6, 1.6, 1.88972614};
218  ions.R[2] = {1.4, 0.0, 0.0};
219  SpeciesSet& ion_species = ions.getSpeciesSet();
220  int pIdx = ion_species.addSpecies("H");
221  int pChargeIdx = ion_species.addAttribute("charge");
222  ion_species(pChargeIdx, pIdx) = 1;
223  ions.createSK();
224 
225  elec.setName("elec");
226  elec.create({2});
227  elec.R[0] = {0.5, 0.0, 0.0};
228  elec.R[1] = {0.0, 0.5, 0.0};
229  SpeciesSet& tspecies = elec.getSpeciesSet();
230  int upIdx = tspecies.addSpecies("u");
231  int chargeIdx = tspecies.addAttribute("charge");
232  int massIdx = tspecies.addAttribute("mass");
233  tspecies(chargeIdx, upIdx) = -1;
234  tspecies(massIdx, upIdx) = 1.0;
235  elec.createSK();
236 
237  // The call to resetGroups is needed transfer the SpeciesSet
238  // settings to the ParticleSet
239  ions.resetGroups();
240  elec.resetGroups();
241 
242  LRCoulombSingleton::CoulombHandler = std::make_unique<EwaldHandler3D>(ions);
243  LRCoulombSingleton::CoulombHandler->initBreakup(ions);
244  LRCoulombSingleton::CoulombDerivHandler = std::make_unique<EwaldHandler3D>(ions);
245  LRCoulombSingleton::CoulombDerivHandler->initBreakup(ions);
246 
247  ForceChiesaPBCAA force(ions, elec);
248  force.setAddIonIon(false);
249 
250  //Ion-Ion forces are validated against Quantum Espresso's ewald method:
251  CHECK(force.getForcesIonIon()[0][0] == Approx(-0.37660901));
252  CHECK(force.getForcesIonIon()[0][1] == Approx(-0.02283659));
253  CHECK(force.getForcesIonIon()[0][2] == Approx(0.0000000));
254  CHECK(force.getForcesIonIon()[1][0] == Approx(0.04012282));
255  CHECK(force.getForcesIonIon()[1][1] == Approx(0.066670175));
256  CHECK(force.getForcesIonIon()[1][2] == Approx(0.0000000));
257  CHECK(force.getForcesIonIon()[2][0] == Approx(0.336486185));
258  CHECK(force.getForcesIonIon()[2][1] == Approx(-0.04383358));
259  CHECK(force.getForcesIonIon()[2][2] == Approx(0.0000000));
260 
261  elec.update();
262  force.InitMatrix();
263  force.evaluate(elec);
264  //Electron-Ion forces are unvalidated externally:
265  CHECK(force.getForces()[0][0] == Approx(3.959178977));
266  CHECK(force.getForces()[0][1] == Approx(3.959178977));
267  CHECK(force.getForces()[0][2] == Approx(0.000000000));
268  CHECK(force.getForces()[1][0] == Approx(-0.078308730));
269  CHECK(force.getForces()[1][1] == Approx(-0.078308730));
270  CHECK(force.getForces()[1][2] == Approx(0.000000000));
271  CHECK(force.getForces()[2][0] == Approx(-1.4341388802));
272  CHECK(force.getForces()[2][1] == Approx(0.1379375923));
273  CHECK(force.getForces()[2][2] == Approx(0.000000000));
274 
275  LRCoulombSingleton::CoulombHandler.reset(nullptr);
276  LRCoulombSingleton::CoulombDerivHandler.reset(nullptr);
277 }
278 
279 } // namespace qmcplusplus
qmcplusplus::CrystalLattice
a class that defines a supercell in D-dimensional Euclean space.
Definition: CrystalLattice.h:55
qmcplusplus::ParticleSet::setName
void setName(const std::string &aname)
Definition: ParticleSet.h:237
qmcplusplus::LRCoulombSingleton::CoulombDerivHandler
static std::unique_ptr< LRHandlerType > CoulombDerivHandler
Stores the force/stress optimized LR handler.
Definition: LRCoulombSingleton.h:52
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::ForceBase::setForces
void setForces(const ParticleSet::ParticlePos &forces)
Definition: ForceBase.cpp:139
qmcplusplus::TEST_CASE
TEST_CASE("complex_helper", "[type_traits]")
Definition: test_complex_helper.cpp: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::LRCoulombSingleton::CoulombHandler
static std::unique_ptr< LRHandlerType > CoulombHandler
Stores the energ optimized LR handler.
Definition: LRCoulombSingleton.h:50
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
Communicate
Wrapping information on parallelism.
Definition: Communicate.h:68
qmcplusplus::lattice
CrystalLattice< OHMMS_PRECISION, OHMMS_DIM > lattice
Definition: test_SpinDensityNew.cpp:142
qmcplusplus::ForceChiesaPBCAA::evaluate
Return_t evaluate(ParticleSet &P) override
Evaluate the local energy contribution of this component.
Definition: ForceChiesaPBCAA.cpp:173
qmcplusplus::ParticleSet
Specialized paritlce class for atomistic simulations.
Definition: ParticleSet.h:55
qmcplusplus::REQUIRE
REQUIRE(std::filesystem::exists(filename))
qmcplusplus::ParticleSetPool
Manage a collection of ParticleSet objects.
Definition: ParticleSetPool.h:34
qmcplusplus::ParticleSet::R
ParticlePos R
Position.
Definition: ParticleSet.h:79
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.
qmcplusplus::ForceBase::getForcesIonIon
const ParticleSet::ParticlePos & getForcesIonIon() const noexcept
Definition: ForceBase.h:68
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 }))
qmcplusplus::ForceChiesaPBCAA::makeClone
std::unique_ptr< OperatorBase > makeClone(ParticleSet &qp, TrialWaveFunction &psi) final
Definition: ForceChiesaPBCAA.cpp:232
qmcplusplus::ParticleSet::createSK
void createSK()
create Structure Factor with PBCs
Definition: ParticleSet.BC.cpp:29
qmcplusplus::SpeciesSet
Custom container for set of attributes for a set of species.
Definition: SpeciesSet.h:33
qmcplusplus::ForceBase::getAddIonIon
bool getAddIonIon() const noexcept
Definition: ForceBase.h:61
qmcplusplus::ForceBase::setAddIonIon
void setAddIonIon(bool val) noexcept
Definition: ForceBase.h:62
ParticleSetPool.h
Declaration of ParticleSetPool.
qmcplusplus::ForceBase::getForces
const ParticleSet::ParticlePos & getForces() const noexcept
Definition: ForceBase.h:64