QMCPACK
test_coulomb_pbcAB_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 "QMCHamiltonians/CoulombPBCAB.h"
19 #include "QMCHamiltonians/CoulombPBCAA.h"
20 #include "LongRange/EwaldHandler3D.h"
21 
22 #include <stdio.h>
23 #include <string>
24 
25 using std::string;
26 
27 namespace qmcplusplus
28 {
29 TEST_CASE("Coulomb PBC A-B Ewald3D", "[hamiltonian]")
30 {
31  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
32  lattice.BoxBConds = true; // periodic
33  lattice.R.diagonal(1.0);
34  lattice.reset();
35 
36  const SimulationCell simulation_cell(lattice);
37  ParticleSet ions(simulation_cell);
38  ParticleSet elec(simulation_cell);
39 
40  ions.setName("ion");
41  ions.create({1});
42  ions.R[0] = {0.0, 0.0, 0.0};
43  SpeciesSet& ion_species = ions.getSpeciesSet();
44  int pIdx = ion_species.addSpecies("H");
45  int pChargeIdx = ion_species.addAttribute("charge");
46  ion_species(pChargeIdx, pIdx) = 1;
47  ions.createSK();
48  ions.update();
49 
50 
51  elec.setName("elec");
52  elec.create({1});
53  elec.R[0] = {0.5, 0.0, 0.0};
54  SpeciesSet& tspecies = elec.getSpeciesSet();
55  int upIdx = tspecies.addSpecies("u");
56  int chargeIdx = tspecies.addAttribute("charge");
57  int massIdx = tspecies.addAttribute("mass");
58  tspecies(chargeIdx, upIdx) = -1;
59  tspecies(massIdx, upIdx) = 1.0;
60 
61  elec.resetGroups();
62  elec.createSK();
63  elec.addTable(ions);
64  elec.update();
65 
66 
67  LRCoulombSingleton::CoulombHandler = std::make_unique<EwaldHandler3D>(ions);
68  LRCoulombSingleton::CoulombHandler->initBreakup(ions);
69 
70 
71  CoulombPBCAB cab(ions, elec);
72 
73  // Self energy plus Background charge term
74  double consts = cab.evalConsts(elec);
75  CHECK(consts == Approx(0.0523598776 * 2)); //not validated
76 
77  double val_ei = cab.evaluate(elec);
78  CHECK(val_ei == Approx(-0.008302 + 0.0523598776 * 2)); //Not validated
79 
80  CoulombPBCAA caa_elec(elec, true, false, false);
81  CoulombPBCAA caa_ion(ions, false, false, false);
82  double val_ee = caa_elec.evaluate(elec);
83  double val_ii = caa_ion.evaluate(ions);
84  double sum = val_ee + val_ii + val_ei;
85 
86  CHECK(val_ee == Approx(-1.418927));
87  CHECK(val_ii == Approx(-1.418927));
88  CHECK(sum == Approx(-2.741436)); // Can be validated via Ewald summation elsewhere
89  // -2.74136517454081
90 
91  LRCoulombSingleton::CoulombHandler.reset(nullptr);
92 }
93 
94 TEST_CASE("Coulomb PBC A-B BCC H Ewald3D", "[hamiltonian]")
95 {
96  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
97  lattice.BoxBConds = true; // periodic
98  lattice.R.diagonal(3.77945227);
99  lattice.reset();
100 
101  const SimulationCell simulation_cell(lattice);
102  ParticleSet ions(simulation_cell);
103  ParticleSet elec(simulation_cell);
104 
105  ions.setName("ion");
106  ions.create({2});
107  ions.R[0] = {0.0, 0.0, 0.0};
108  ions.R[1] = {1.88972614, 1.88972614, 1.88972614};
109  SpeciesSet& ion_species = ions.getSpeciesSet();
110  int pIdx = ion_species.addSpecies("H");
111  int pChargeIdx = ion_species.addAttribute("charge");
112  ion_species(pChargeIdx, pIdx) = 1;
113  ions.createSK();
114  ions.update();
115 
116 
117  elec.setName("elec");
118  elec.create({2});
119  elec.R[0] = {0.5, 0.0, 0.0};
120  elec.R[1] = {0.0, 0.5, 0.0};
121  SpeciesSet& tspecies = elec.getSpeciesSet();
122  int upIdx = tspecies.addSpecies("u");
123  int chargeIdx = tspecies.addAttribute("charge");
124  int massIdx = tspecies.addAttribute("mass");
125  tspecies(chargeIdx, upIdx) = -1;
126  tspecies(massIdx, upIdx) = 1.0;
127 
128  elec.resetGroups();
129  elec.createSK();
130  elec.addTable(ions);
131  elec.update();
132 
133 
134  LRCoulombSingleton::CoulombHandler = std::make_unique<EwaldHandler3D>(ions);
135  LRCoulombSingleton::CoulombHandler->initBreakup(ions);
136 
137  CoulombPBCAB cab(ions, elec);
138 
139  // Background charge term
140  double consts = cab.evalConsts(elec);
141  CHECK(consts == Approx(0.0277076538 * 4)); //not validated
142 
143 
144  double val_ei = cab.evaluate(elec);
145  CHECK(val_ei == Approx(-2.223413 + 0.0277076538 * 4)); //Not validated
146 
147 
148  CoulombPBCAA caa_elec(elec, false, false, false);
149  CoulombPBCAA caa_ion(ions, false, false, false);
150  double val_ee = caa_elec.evaluate(elec);
151  double val_ii = caa_ion.evaluate(ions);
152  double sum = val_ee + val_ii + val_ei;
153 
154  CHECK(val_ee == Approx(-0.012808 - 0.0277076538 * 2));
155  CHECK(val_ii == Approx(-0.907659 - 0.0277076538 * 2));
156  CHECK(sum == Approx(-3.143880)); // Can be validated via Ewald summation elsewhere
157  // -3.14349127313640
158 
159  LRCoulombSingleton::CoulombHandler.reset(nullptr);
160 }
161 
162 } // 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::CoulombPBCAB::evaluate
Return_t evaluate(ParticleSet &P) override
Evaluate the local energy contribution of this component.
Definition: CoulombPBCAB.cpp:122
qmcplusplus::CoulombPBCAA
Calculates the AA Coulomb potential using PBCs.
Definition: CoulombPBCAA.h:38
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::CoulombPBCAB
Calculates the AA Coulomb potential using PBCs.
Definition: CoulombPBCAB.h:38
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
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
qmcplusplus::lattice
CrystalLattice< OHMMS_PRECISION, OHMMS_DIM > lattice
Definition: test_SpinDensityNew.cpp:142
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::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
qmcplusplus::CHECK
CHECK(log_values[0]==ComplexApprox(std::complex< double >{ 5.603777579195571, -6.1586603331188225 }))
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::CoulombPBCAA::evaluate
Return_t evaluate(ParticleSet &P) override
Evaluate the local energy contribution of this component.
Definition: CoulombPBCAA.cpp:201
qmcplusplus::CoulombPBCAB::evalConsts
Return_t evalConsts(const ParticleSet &P, bool report=true)
Evaluates madelung and background contributions to total energy.
Definition: CoulombPBCAB.cpp:390