QMCPACK
test_ewald_quasi2d.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 
16 namespace qmcplusplus
17 {
18 
19 TEST_CASE("Coulomb PBC A-A Ewald Quasi2D exception", "[hamiltonian]")
20 {
21  LRCoulombSingleton::CoulombHandler = 0;
22  LRCoulombSingleton::this_lr_type = LRCoulombSingleton::QUASI2D;
23  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
24  lattice.BoxBConds = true;
25  lattice.BoxBConds[2] = false; // ppn
26  lattice.ndim = 2;
27  lattice.R.diagonal(1.0);
28  lattice.LR_dim_cutoff = 1.0;
29  lattice.reset();
30 
31  const SimulationCell simulation_cell(lattice);
32  ParticleSet elec(simulation_cell);
33  elec.setName("e");
34  elec.create({1});
35  elec.R[0] = {0.0, 0.0, 0.0};
36 
37  SpeciesSet& tspecies = elec.getSpeciesSet();
38  int upIdx = tspecies.addSpecies("u");
39  int chargeIdx = tspecies.addAttribute("charge");
40  int massIdx = tspecies.addAttribute("mass");
41  tspecies(chargeIdx, upIdx) = -1;
42  tspecies(massIdx, upIdx) = 1.0;
43 
44  elec.createSK();
45  elec.update();
46 
47  CHECK_THROWS(CoulombPBCAA(elec, true, false, false));
48 }
49 
50 TEST_CASE("Coulomb PBC A-A Ewald Quasi2D square", "[hamiltonian]")
51 {
52  LRCoulombSingleton::CoulombHandler = 0;
53  LRCoulombSingleton::this_lr_type = LRCoulombSingleton::QUASI2D;
54  const double vmad_sq = -1.95013246;
55  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
56  lattice.BoxBConds = true;
57  lattice.BoxBConds[2] = false; // ppn
58  lattice.ndim = 2;
59  lattice.R.diagonal(1.0);
60  lattice.LR_dim_cutoff = 30.0;
61  lattice.reset();
62 
63  const SimulationCell simulation_cell(lattice);
64  ParticleSet elec(simulation_cell);
65  elec.setName("e");
66  elec.create({1});
67  elec.R[0] = {0.0, 0.0, 0.0};
68 
69  SpeciesSet& tspecies = elec.getSpeciesSet();
70  int upIdx = tspecies.addSpecies("u");
71  int chargeIdx = tspecies.addAttribute("charge");
72  int massIdx = tspecies.addAttribute("mass");
73  tspecies(chargeIdx, upIdx) = -1;
74  tspecies(massIdx, upIdx) = 1.0;
75 
76  elec.createSK();
77  elec.update();
78 
79  CoulombPBCAA caa(elec, true, false, false);
80 
81  double val = caa.evaluate(elec);
82  CHECK(val == Approx(vmad_sq));
83 }
84 
85 TEST_CASE("Coulomb PBC A-A Ewald Quasi2D triangular", "[hamiltonian]")
86 {
87  LRCoulombSingleton::CoulombHandler = 0; // !!!! crucial if not first test
88  LRCoulombSingleton::this_lr_type = LRCoulombSingleton::QUASI2D;
89  const double vmad_tri = -1.106102587;
90  const double alat = std::sqrt(2.0*M_PI/std::sqrt(3));
91 
92  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
93  lattice.BoxBConds = true;
94  lattice.BoxBConds[2] = false; // ppn
95  lattice.ndim = 2;
96  lattice.R = 0.0;
97  lattice.R(0, 0) = alat;
98  lattice.R(1, 0) = -1.0/2*alat;
99  lattice.R(1, 1) = std::sqrt(3)/2*alat;
100  lattice.R(2, 2) = 2*alat;
101  lattice.LR_dim_cutoff = 30.0;
102  lattice.reset();
103 
104  const SimulationCell simulation_cell(lattice);
105  ParticleSet elec(simulation_cell);
106  elec.setName("e");
107  elec.create({1});
108  elec.R[0] = {0.0, 0.0, 0.0};
109 
110  SpeciesSet& tspecies = elec.getSpeciesSet();
111  int upIdx = tspecies.addSpecies("u");
112  int chargeIdx = tspecies.addAttribute("charge");
113  int massIdx = tspecies.addAttribute("mass");
114  tspecies(chargeIdx, upIdx) = -1;
115  tspecies(massIdx, upIdx) = 1.0;
116 
117  elec.createSK();
118  elec.update();
119 
120  CoulombPBCAA caa(elec, true, false, false);
121 
122  double val = caa.evaluate(elec);
123  CHECK(val == Approx(vmad_tri));
124 }
125 
126 TEST_CASE("Coulomb PBC A-A Ewald Quasi2D staggered square", "[hamiltonian]")
127 {
128  LRCoulombSingleton::CoulombHandler = 0; // !!!! crucial if not first test
129  LRCoulombSingleton::this_lr_type = LRCoulombSingleton::QUASI2D;
130  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
131  lattice.BoxBConds = true;
132  lattice.BoxBConds[2] = false; // ppn
133  lattice.ndim = 2;
134  lattice.R = 0.0;
135  lattice.R.diagonal(1.0);
136  lattice.R(2,2) = 10;
137  lattice.LR_dim_cutoff = 30.0;
138  lattice.reset();
139 
140  const SimulationCell simulation_cell(lattice);
141  ParticleSet elec(simulation_cell);
142  int npart = 2;
143  elec.setName("e");
144  elec.create({npart});
145  elec.R[0] = {0.0, 0.0, 0.0};
146  elec.R[1] = {0.5, 0.5, 0.0};
147 
148  SpeciesSet& tspecies = elec.getSpeciesSet();
149  int upIdx = tspecies.addSpecies("u");
150  int chargeIdx = tspecies.addAttribute("charge");
151  int massIdx = tspecies.addAttribute("mass");
152  tspecies(chargeIdx, upIdx) = -1;
153  tspecies(massIdx, upIdx) = 1.0;
154 
155  elec.createSK();
156  elec.addTable(elec); // !!!! crucial to compute distance table
157 
158  CoulombPBCAA caa(elec, true, false, false);
159 
160  const int ntest = 4;
161  const double zheight[ntest] = {0, 0.1, 0.5, 3.0};
162  const double vmad_sq = -1.95013246;
163  const double vmad_bc = -2.7579038;
164  const double vmad_ref[ntest] = {vmad_bc, -2.4846003745840357, -2.019557388069959, vmad_sq};
165  double val;
166  for (int itest=0; itest<ntest; itest++)
167  {
168  elec.R[1][2] = zheight[itest];
169  elec.update();
170  val = caa.evaluate(elec);
171  CHECK(val/npart == Approx(vmad_ref[itest]));
172  }
173 }
174 
175 TEST_CASE("Coulomb PBC A-A Ewald Quasi2D staggered square 2x2", "[hamiltonian]")
176 {
177  LRCoulombSingleton::CoulombHandler = 0; // !!!! crucial if not first test
178  LRCoulombSingleton::this_lr_type = LRCoulombSingleton::QUASI2D;
179  //const double vmad_sq = -1.95013246;
180  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
181  lattice.BoxBConds = true;
182  lattice.BoxBConds[2] = false; // ppn
183  lattice.ndim = 2;
184  lattice.R.diagonal(2.0);
185  lattice.R(2,2) = 10;
186  lattice.LR_dim_cutoff = 30.0;
187  lattice.reset();
188 
189  const SimulationCell simulation_cell(lattice);
190  ParticleSet elec(simulation_cell);
191  int npart = 8;
192  elec.setName("e");
193  elec.create({npart});
194  elec.R[0] = {0.0, 0.0, 0.0};
195  elec.R[1] = {1.0, 0.0, 0.0};
196  elec.R[2] = {0.0, 1.0, 0.0};
197  elec.R[3] = {1.0, 1.0, 0.0};
198  elec.R[4] = {0.5, 0.5, 0.0};
199  elec.R[5] = {1.5, 0.5, 0.0};
200  elec.R[6] = {0.5, 1.5, 0.0};
201  elec.R[7] = {1.5, 1.5, 0.0};
202 
203  SpeciesSet& tspecies = elec.getSpeciesSet();
204  int upIdx = tspecies.addSpecies("u");
205  int chargeIdx = tspecies.addAttribute("charge");
206  int massIdx = tspecies.addAttribute("mass");
207  tspecies(chargeIdx, upIdx) = -1;
208  tspecies(massIdx, upIdx) = 1.0;
209 
210  elec.createSK();
211  elec.addTable(elec); // !!!! crucial to compute distance table
212 
213  CoulombPBCAA caa(elec, true, false, false);
214 
215  const int ntest = 4;
216  const double zheight[ntest] = {0, 0.1, 0.5, 3.0};
217  const double vmad_sq = -1.95013246;
218  const double vmad_bc = -2.7579038;
219  const double vmad_ref[ntest] = {vmad_bc, -2.4846003745840357, -2.019557388069959, vmad_sq};
220  double val;
221  for (int itest=0; itest<ntest; itest++)
222  {
223  for (int i=npart/2;i<npart;i++)
224  elec.R[i][2] = zheight[itest];
225  elec.update();
226  val = caa.evaluate(elec);
227  CHECK(val/npart == Approx(vmad_ref[itest]));
228  }
229 }
230 
231 TEST_CASE("Coulomb PBC A-A Ewald Quasi2D staggered triangle", "[hamiltonian]")
232 {
233  LRCoulombSingleton::CoulombHandler = 0; // !!!! crucial if not first test
234  LRCoulombSingleton::this_lr_type = LRCoulombSingleton::QUASI2D;
235  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
236  lattice.BoxBConds = true;
237  lattice.BoxBConds[2] = false; // ppn
238  lattice.ndim = 2;
239  const double alat = std::sqrt(2.0*M_PI/std::sqrt(3));
240  lattice.R = 0.0;
241  lattice.R(0, 0) = alat;
242  lattice.R(1, 0) = -1.0/2*alat;
243  lattice.R(1, 1) = std::sqrt(3)/2*alat;
244  lattice.R(2, 2) = 2*alat;
245  lattice.LR_dim_cutoff = 30.0;
246  lattice.reset();
247 
248  const SimulationCell simulation_cell(lattice);
249  ParticleSet elec(simulation_cell);
250  int npart = 2;
251  elec.setName("e");
252  elec.create({npart});
253  // initialize fractional coordinates
254  elec.R[0] = {0.0, 0.0, 0.0};
255  elec.R[1] = {2./3, 1./3, 0.0};
256  // convert to Cartesian coordinates
257  for (int i=0;i<npart;i++)
258  elec.R[i] = dot(elec.R[i], lattice.R);
259 
260  SpeciesSet& tspecies = elec.getSpeciesSet();
261  int upIdx = tspecies.addSpecies("u");
262  int chargeIdx = tspecies.addAttribute("charge");
263  int massIdx = tspecies.addAttribute("mass");
264  tspecies(chargeIdx, upIdx) = -1;
265  tspecies(massIdx, upIdx) = 1.0;
266 
267  elec.createSK();
268  elec.addTable(elec); // !!!! crucial to compute distance table
269 
270  CoulombPBCAA caa(elec, true, false, false);
271 
272  const int ntest = 4;
273  const double zheight[ntest] = {0, 0.1, 0.5, 3.0};
274  const double vmad_hon = -1.510964233;
275  const double vmad_tri = -1.106102587;
276  const double vmad_ref[ntest] = {vmad_hon, -1.4193042644, -1.2005504968, vmad_tri};
277  double val;
278  for (int itest=0; itest<ntest; itest++)
279  {
280  elec.R[1][2] = zheight[itest];
281  elec.update();
282  val = caa.evaluate(elec);
283  CHECK(val/npart == Approx(vmad_ref[itest]));
284  }
285 }
286 
287 TEST_CASE("Coulomb PBC A-A Ewald Quasi2D staggered triangle 2x2", "[hamiltonian]")
288 {
289  LRCoulombSingleton::CoulombHandler = 0; // !!!! crucial if not first test
290  LRCoulombSingleton::this_lr_type = LRCoulombSingleton::QUASI2D;
291  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
292  lattice.BoxBConds = true;
293  lattice.BoxBConds[2] = false; // ppn
294  lattice.ndim = 2;
295  const double rs = 10;
296  const double alat = rs*2*std::sqrt(2.0*M_PI/std::sqrt(3));
297  lattice.R = 0.0;
298  lattice.R(0, 0) = alat;
299  lattice.R(1, 0) = -1.0/2*alat;
300  lattice.R(1, 1) = std::sqrt(3)/2*alat;
301  lattice.R(2, 2) = 2*alat;
302  lattice.LR_dim_cutoff = 30.0;
303  lattice.reset();
304 
305  const SimulationCell simulation_cell(lattice);
306  ParticleSet elec(simulation_cell);
307  int npart = 8;
308  elec.setName("e");
309  elec.create({npart});
310  // initialize fractional coordinates
311  TinyVector<double, 3> r0 = {0.0, 0.0, 0.0};
312  TinyVector<double, 3> r1 = {2./3, 1./3, 0.0};
313  elec.R[0] = {0.0, 0.0, 0.0};
314  elec.R[1] = {0.5, 0.0, 0.0};
315  elec.R[2] = {0.0, 0.5, 0.0};
316  elec.R[3] = {0.5, 0.5, 0.0};
317  elec.R[4] = {0.0, 0.0, 0.0};
318  elec.R[5] = {0.5, 0.0, 0.0};
319  elec.R[6] = {0.0, 0.5, 0.0};
320  elec.R[7] = {0.5, 0.5, 0.0};
321  for (int i=0;i<npart/2;i++)
322  elec.R[i] += r0/2;
323  for (int i=npart/2;i<npart;i++)
324  elec.R[i] += r1/2;
325  // convert to Cartesian coordinates
326  for (int i=0;i<npart;i++)
327  {
328  elec.R[i] = dot(elec.R[i], lattice.R);
329  }
330 
331  SpeciesSet& tspecies = elec.getSpeciesSet();
332  int upIdx = tspecies.addSpecies("u");
333  int chargeIdx = tspecies.addAttribute("charge");
334  int massIdx = tspecies.addAttribute("mass");
335  tspecies(chargeIdx, upIdx) = -1;
336  tspecies(massIdx, upIdx) = 1.0;
337 
338  elec.createSK();
339  elec.addTable(elec); // !!!! crucial to compute distance table
340 
341  CoulombPBCAA caa(elec, true, false, false);
342 
343  const int ntest = 4;
344  TinyVector<double, ntest> zheight = {0, 0.1, 0.5, 3.0};
345  zheight *= rs;
346  const double vmad_hon = -1.510964233;
347  const double vmad_tri = -1.106102587;
348  TinyVector<double, ntest> vmad_ref = {vmad_hon, -1.4193042644, -1.2005504968, vmad_tri};
349  vmad_ref /= rs;
350  double val;
351  for (int itest=0; itest<ntest; itest++)
352  {
353  for (int i=npart/2;i<npart;i++)
354  elec.R[i][2] = zheight[itest];
355  elec.update();
356  val = caa.evaluate(elec);
357  CHECK(val/npart == Approx(vmad_ref[itest]));
358  }
359 }
360 
361 } // 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::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::for
for(int i=0;i< size_test;++i) CHECK(Approx(gauss_random_vals[offset_for_rs+i])
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