QMCPACK
test_ewald2d.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) 2022 QMCPACK developers.
6 //
7 // File developed by: Yubo "Paul" Yang, yubo.paul.yang@gmail.com, CCQ @ Flatiron
8 //
9 // File created by: Yubo "Paul" Yang, yubo.paul.yang@gmail.com, CCQ @ Flatiron
10 //////////////////////////////////////////////////////////////////////////////////////
11 #include "catch.hpp"
12 
13 #include "Particle/ParticleSet.h"
14 #include "QMCHamiltonians/CoulombPBCAA.h"
15 #include "LongRange/EwaldHandler2D.h"
16 
17 namespace qmcplusplus
18 {
19 
20 TEST_CASE("Coulomb PBC A-A Ewald2D square", "[hamiltonian]")
21 {
22  LRCoulombSingleton::CoulombHandler = 0;
23  LRCoulombSingleton::this_lr_type = LRCoulombSingleton::STRICT2D;
24  const double vmad_sq = -1.95013246;
25  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
26  lattice.BoxBConds = true; // periodic
27  lattice.ndim = 2;
28  lattice.R.diagonal(1.0);
29  lattice.LR_dim_cutoff = 30.0;
30  lattice.reset();
31 
32  const SimulationCell simulation_cell(lattice);
33  ParticleSet elec(simulation_cell);
34  elec.setName("e");
35  elec.create({1});
36  elec.R[0] = {0.0, 0.0, 0.0};
37 
38  SpeciesSet& tspecies = elec.getSpeciesSet();
39  int upIdx = tspecies.addSpecies("u");
40  int chargeIdx = tspecies.addAttribute("charge");
41  int massIdx = tspecies.addAttribute("mass");
42  tspecies(chargeIdx, upIdx) = -1;
43  tspecies(massIdx, upIdx) = 1.0;
44 
45  elec.createSK();
46  elec.update();
47 
48  CoulombPBCAA caa = CoulombPBCAA(elec, true, false, false);
49 
50  double val = caa.evaluate(elec);
51  CHECK(val == Approx(vmad_sq));
52 }
53 
54 TEST_CASE("Coulomb PBC A-A Ewald2D body center", "[hamiltonian]")
55 {
56  LRCoulombSingleton::CoulombHandler = 0; // !!!! crucial if not first test
57  LRCoulombSingleton::this_lr_type = LRCoulombSingleton::STRICT2D;
58  const double vmad_bc = -2.7579038;
59 
60  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
61  lattice.BoxBConds = true; // periodic
62  lattice.ndim = 2;
63  lattice.R = 0.0;
64  lattice.R.diagonal(1.0);
65  lattice.LR_dim_cutoff = 30.0;
66  lattice.reset();
67 
68  const SimulationCell simulation_cell(lattice);
69  ParticleSet elec(simulation_cell);
70  elec.setName("e");
71  const int npart = 2;
72  elec.create({npart});
73  // initialize fractional coordinates
74  elec.R[0] = {0.0, 0.0, 0.0};
75  elec.R[1] = {0.5, 0.5, 0.0};
76  // convert to Cartesian coordinates
77  for (int i=0;i<npart;i++)
78  elec.R[i] = dot(elec.R[i], lattice.R);
79 
80  SpeciesSet& tspecies = elec.getSpeciesSet();
81  int upIdx = tspecies.addSpecies("u");
82  int chargeIdx = tspecies.addAttribute("charge");
83  int massIdx = tspecies.addAttribute("mass");
84  tspecies(chargeIdx, upIdx) = -1;
85  tspecies(massIdx, upIdx) = 1.0;
86 
87  elec.createSK();
88  elec.addTable(elec); // !!!! crucial with more than 1 particle
89  elec.update();
90 
91  CoulombPBCAA caa = CoulombPBCAA(elec, true, false, false);
92 
93  double val = caa.evaluate(elec);
94  CHECK(val/npart == Approx(vmad_bc));
95 }
96 
97 TEST_CASE("Coulomb PBC A-A Ewald2D triangular", "[hamiltonian]")
98 {
99  LRCoulombSingleton::CoulombHandler = 0; // !!!! crucial if not first test
100  LRCoulombSingleton::this_lr_type = LRCoulombSingleton::STRICT2D;
101  const double vmad_tri = -1.1061025865191676;
102  const double alat = std::sqrt(2.0*M_PI/std::sqrt(3));
103 
104  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
105  lattice.BoxBConds = true; // periodic
106  lattice.ndim = 2;
107  lattice.R = 0.0;
108  lattice.R(0, 0) = alat;
109  lattice.R(1, 0) = -1.0/2*alat;
110  lattice.R(1, 1) = std::sqrt(3)/2*alat;
111  lattice.R(2, 2) = 2*alat;
112  lattice.LR_dim_cutoff = 30.0;
113  lattice.reset();
114 
115  const SimulationCell simulation_cell(lattice);
116  ParticleSet elec(simulation_cell);
117  elec.setName("e");
118  elec.create({1});
119  elec.R[0] = {0.0, 0.0, 0.0};
120 
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.createSK();
129  elec.update();
130 
131  CoulombPBCAA caa = CoulombPBCAA(elec, true, false, false);
132 
133  double val = caa.evaluate(elec);
134  CHECK(val == Approx(vmad_tri));
135 }
136 
137 TEST_CASE("Coulomb PBC A-A Ewald2D honeycomb", "[hamiltonian]")
138 {
139  LRCoulombSingleton::CoulombHandler = 0; // !!!! crucial if not first test
140  LRCoulombSingleton::this_lr_type = LRCoulombSingleton::STRICT2D;
141  const double vmad_hon = -1.510964233;
142  const double alat = std::sqrt(2.0*M_PI/std::sqrt(3));
143 
144  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
145  lattice.BoxBConds = true; // periodic
146  lattice.ndim = 2;
147  lattice.R = 0.0;
148  lattice.R(0, 0) = alat;
149  lattice.R(1, 0) = -1.0/2*alat;
150  lattice.R(1, 1) = std::sqrt(3)/2*alat;
151  lattice.R(2, 2) = 2*alat;
152  lattice.LR_dim_cutoff = 30.0;
153  lattice.reset();
154 
155  const SimulationCell simulation_cell(lattice);
156  ParticleSet elec(simulation_cell);
157  elec.setName("e");
158  const int npart = 2;
159  elec.create({npart});
160  // initialize fractional coordinates
161  elec.R[0] = {0.0, 0.0, 0.0};
162  elec.R[1] = {2./3, 1./3, 0.0};
163  // convert to Cartesian coordinates
164  for (int i=0;i<npart;i++)
165  elec.R[i] = dot(elec.R[i], lattice.R);
166 
167  SpeciesSet& tspecies = elec.getSpeciesSet();
168  int upIdx = tspecies.addSpecies("u");
169  int chargeIdx = tspecies.addAttribute("charge");
170  int massIdx = tspecies.addAttribute("mass");
171  tspecies(chargeIdx, upIdx) = -1;
172  tspecies(massIdx, upIdx) = 1.0;
173 
174  elec.createSK();
175  elec.addTable(elec); // !!!! crucial with more than 1 particle
176  elec.update();
177 
178  CoulombPBCAA caa = CoulombPBCAA(elec, true, false, false);
179 
180  double val = caa.evaluate(elec);
181  CHECK(val/npart == Approx(vmad_hon));
182 }
183 
184 TEST_CASE("Coulomb PBC A-A Ewald2D tri. in rect.", "[hamiltonian]")
185 {
186  using RealType = QMCTraits::RealType;
187  LRCoulombSingleton::CoulombHandler = 0; // !!!! crucial if not first test
188  LRCoulombSingleton::this_lr_type = LRCoulombSingleton::STRICT2D;
189  const double vmad_tri = -1.1061025865191676;
190  const RealType alat = std::sqrt(2.0*M_PI/std::sqrt(3));
191 
192  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
193  lattice.BoxBConds = true; // periodic
194  lattice.ndim = 2;
195  lattice.R = 0.0;
196  lattice.R(0, 0) = alat;
197  lattice.R(1, 1) = std::sqrt(3)*alat;
198  lattice.R(2, 2) = 4.0;
199  lattice.LR_dim_cutoff = 30.0;
200  lattice.reset();
201  lattice.print(app_log());
202 
203  const SimulationCell simulation_cell(lattice);
204  ParticleSet elec(simulation_cell);
205  elec.setName("e");
206  elec.create({2});
207  elec.R[0] = {0.0, 0.0, 0.0};
208  elec.R[1] = {alat/2, static_cast<RealType>(std::sqrt(3))*alat/2, 0.0};
209 
210  SpeciesSet& tspecies = elec.getSpeciesSet();
211  int upIdx = tspecies.addSpecies("u");
212  int chargeIdx = tspecies.addAttribute("charge");
213  int massIdx = tspecies.addAttribute("mass");
214  tspecies(chargeIdx, upIdx) = -1;
215  tspecies(massIdx, upIdx) = 1.0;
216 
217  elec.createSK();
218  elec.addTable(elec); // !!!! crucial with more than 1 particle
219  elec.update();
220  app_log() << elec.R << std::endl;
221 
222  CoulombPBCAA caa = CoulombPBCAA(elec, true, false, false);
223 
224  double val = caa.evaluate(elec)/elec.getTotalNum();
225  CHECK(val == Approx(vmad_tri));
226 }
227 } // 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::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::QMCTraits::RealType
QTBase::RealType RealType
Definition: Configuration.h:58
qmcplusplus::ParticleSet::getTotalNum
size_t getTotalNum() const
Definition: ParticleSet.h:493
qmcplusplus::app_log
std::ostream & app_log()
Definition: OutputManager.h:65
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::sqrt
MakeReturn< UnaryNode< FnSqrt, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t sqrt(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:136
qmcplusplus::CHECK
CHECK(log_values[0]==ComplexApprox(std::complex< double >{ 5.603777579195571, -6.1586603331188225 }))
qmcplusplus::dot
Tensor< typename BinaryReturn< T1, T2, OpMultiply >::Type_t, D > dot(const AntiSymTensor< T1, D > &lhs, const AntiSymTensor< T2, D > &rhs)
Definition: AntiSymTensor.h:520
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