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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
8 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
9 // Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory
10 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
11 //
12 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
13 //////////////////////////////////////////////////////////////////////////////////////
14
15
16 #ifndef QMCPLUSPLUS_HEGGRID_H
17 #define QMCPLUSPLUS_HEGGRID_H
18
19 #include "Lattice/CrystalLattice.h"
20 #include <map>
21
22 namespace qmcplusplus
23 {
24 template<typename T, unsigned D>
25 struct kpdata
26 {
27 TinyVector<T, D> k;
28 T k2;
29 int g;
30 };
31
32
33 template<typename T, unsigned D>
kpdata_comp(const kpdata<T,D> & left,const kpdata<T,D> & right)34 bool kpdata_comp(const kpdata<T, D>& left, const kpdata<T, D>& right)
35 {
36 return left.k2 < right.k2;
37 }
38
39
40 template<class T, unsigned D>
41 struct HEGGrid
42 {};
43
44
45 //three-d specialization
46 template<class T>
47 struct HEGGrid<T, 3>
48 {
49 typedef CrystalLattice<T, 3> PL_t;
50 typedef typename PL_t::SingleParticlePos_t PosType;
51 typedef typename PL_t::Scalar_t RealType;
52
53 ///number of kpoints of a half sphere excluding gamma
54 int NumKptsHalf;
55 ///maxmim ksq
56 T MaxKsq;
57 PL_t& Lattice;
58 std::map<int, std::vector<PosType>*> rs;
59
60
61 std::vector<PosType> kpt;
62 std::vector<T> mk2;
63 std::vector<int> deg;
64 std::vector<int> n_within_shell;
65 PosType twist;
66
67
68 typedef kpdata<T, 3> kpdata_t;
69 typedef std::vector<kpdata_t> kpoints_t;
70
71 kpoints_t* kpoints_grid;
72 int nctmp;
73
74
75 HEGGrid(PL_t& lat) : Lattice(lat), twist(0.0), kpoints_grid(0), nctmp(-1)
76 {
77 n_within_shell.resize(31);
78 n_within_shell[0] = 1;
79 n_within_shell[1] = 7;
80 n_within_shell[2] = 19;
81 n_within_shell[3] = 27;
82 n_within_shell[4] = 33;
83 n_within_shell[5] = 57;
84 n_within_shell[6] = 81;
85 n_within_shell[7] = 93;
86 n_within_shell[8] = 123;
87 n_within_shell[9] = 147;
88 n_within_shell[10] = 171;
89 n_within_shell[11] = 179;
90 n_within_shell[12] = 203;
91 n_within_shell[13] = 251;
92 n_within_shell[14] = 257;
93 n_within_shell[15] = 305;
94 n_within_shell[16] = 341;
95 n_within_shell[17] = 365;
96 n_within_shell[18] = 389;
97 n_within_shell[19] = 437;
98 n_within_shell[20] = 461;
99 n_within_shell[21] = 485;
100 n_within_shell[22] = 515;
101 n_within_shell[23] = 587;
102 n_within_shell[24] = 619;
103 n_within_shell[25] = 691;
104 n_within_shell[26] = 739;
105 n_within_shell[27] = 751;
106 n_within_shell[28] = 799;
107 n_within_shell[29] = 847;
108 n_within_shell[30] = 895;
109 }
110
111 ~HEGGrid()
112 {
113 typename std::map<int, std::vector<PosType>*>::iterator it(rs.begin()), it_end(rs.end());
114 while (it != it_end)
115 {
116 delete (*it).second;
117 ++it;
118 }
119 clear_kpoints();
120 }
121
122 /** return the estimated number of grid in each direction */
123 inline int getNC(int nup) const { return static_cast<int>(std::pow(static_cast<T>(nup), 1.0 / 3.0)) / 2 + 1; }
124
125 /** return the estimated number of grid in each direction (upper bound) */
126 inline int get_nc(int nstates) const
127 {
128 return static_cast<int>(std::pow(static_cast<T>(nstates), 1.0 / 3.0) * .7) + 1;
129 }
130
131 //return the number of k-points upto nsh-shell
132 inline int getNumberOfKpoints(int nsh) const
133 {
134 if (nsh < n_within_shell.size())
135 return n_within_shell[nsh];
136 else
137 return -1;
138 }
139
140
141 inline int getShellFromStates(int nst)
142 {
143 for (int i = 0; i < n_within_shell.size(); i++)
144 if (n_within_shell[i] == nst)
145 return i;
146 return -1;
147 }
148
149 //return the shell index for nkpt k-points
150 inline int getShellIndex(int nkpt) const
151 {
152 std::vector<int>::const_iterator loc = std::upper_bound(n_within_shell.begin(), n_within_shell.end(), nkpt);
153 if (loc < n_within_shell.end())
154 return loc - n_within_shell.begin() - 1;
155 else
156 return getNC(nkpt);
157 }
158
159 /** return the cell size for the number of particles and rs
160 * @param nptcl number of particles
161 * @param rs_in rs
162 */
163 inline T getCellLength(int nptcl, T rs_in) const { return std::pow(4.0 * M_PI * nptcl / 3.0, 1.0 / 3.0) * rs_in; }
164
165 void sortGrid(int nc)
166 {
167 int first_ix2, first_ix3;
168 for (int ix1 = 0; ix1 <= nc; ix1++)
169 {
170 if (ix1 == 0)
171 first_ix2 = 0;
172 else
173 first_ix2 = -nc;
174 for (int ix2 = first_ix2; ix2 <= nc; ix2++)
175 {
176 if (ix1 == 0 && ix2 == 0)
177 first_ix3 = 1;
178 else
179 first_ix3 = -nc;
180 for (int ix3 = first_ix3; ix3 <= nc; ix3++)
181 {
182 int ih = ix1 * ix1 + ix2 * ix2 + ix3 * ix3;
183 typename std::map<int, std::vector<PosType>*>::iterator it = rs.find(ih);
184 if (it == rs.end())
185 {
186 std::vector<PosType>* ns = new std::vector<PosType>;
187 ns->push_back(PosType(ix1, ix2, ix3));
188 rs[ih] = ns;
189 }
190 else
191 {
192 (*it).second->push_back(PosType(ix1, ix2, ix3));
193 }
194 }
195 }
196 }
197 }
198
199 void createGrid(int nc, int nkpts)
200 {
201 if (rs.empty())
202 sortGrid(nc);
203 NumKptsHalf = nkpts;
204 kpt.resize(nkpts);
205 mk2.resize(nkpts);
206 int ikpt = 0;
207 //int checkNum=0;
208 //int ke=0;
209 MaxKsq = 0.0;
210 int rsbin = 0;
211 typename std::map<int, std::vector<PosType>*>::iterator rs_it(rs.begin()), rs_end(rs.end());
212 while (ikpt < nkpts && rs_it != rs_end)
213 {
214 typename std::vector<PosType>::iterator ns_it((*rs_it).second->begin()), ns_end((*rs_it).second->end());
215 T minus_ksq = -Lattice.ksq(*ns_it);
216 while (ikpt < nkpts && ns_it != ns_end)
217 {
218 kpt[ikpt] = Lattice.k_cart(*ns_it);
219 mk2[ikpt] = minus_ksq;
220 ++ikpt;
221 ++ns_it;
222 }
223 ++rsbin;
224 ++rs_it;
225 }
226 MaxKsq = Lattice.ksq(*((*rs_it).second->begin()));
227 app_log() << "List of kpoints (half-sphere) " << std::endl;
228 for (int ik = 0; ik < kpt.size(); ik++)
229 {
230 app_log() << ik << " " << kpt[ik] << " " << mk2[ik] << std::endl;
231 }
232 }
233
234
235 void clear_kpoints()
236 {
237 if (kpoints_grid != 0)
238 delete kpoints_grid;
239 }
240
241
242 void create_kpoints(int nc, const PosType& tw, T tol = 1e-6)
243 {
244 if (kpoints_grid == 0)
245 kpoints_grid = new kpoints_t;
246 else if (nc <= nctmp)
247 return;
248 nctmp = nc;
249 kpoints_t& kpoints = *kpoints_grid;
250 app_log() << " resizing kpoint grid" << std::endl;
251 app_log() << " current size = " << kpoints.size() << std::endl;
252 // make space for the kpoint grid
253 int nkpoints = pow(2 * (nc + 1) + 1, 3);
254 kpoints.resize(nkpoints);
255 app_log() << " cubic size = " << kpoints.size() << std::endl;
256 typename kpoints_t::iterator kptmp, kp = kpoints.begin(), kp_end = kpoints.end();
257 // make the kpoint grid
258 T k2max = std::numeric_limits<RealType>::max();
259 for (int i0 = -nc - 1; i0 <= nc + 1; ++i0)
260 for (int i1 = -nc - 1; i1 <= nc + 1; ++i1)
261 for (int i2 = -nc - 1; i2 <= nc + 1; ++i2)
262 {
263 PosType k(i0 + tw[0], i1 + tw[1], i2 + tw[2]);
264 kp->k = Lattice.k_cart(k);
265 kp->k2 = Lattice.ksq(k);
266 if (std::abs(i0) == (nc + 1) || std::abs(i1) == (nc + 1) || std::abs(i2) == (nc + 1))
267 k2max = std::min(k2max, kp->k2);
268 ++kp;
269 }
270 // sort kpoints by magnitude
271 sort(kpoints.begin(), kpoints.end(), kpdata_comp<T, 3>);
272 // eliminate kpoints outside of inscribing sphere
273 int nkp = 0;
274 kp = kpoints.begin();
275 while (kp != kp_end && kp->k2 < k2max + 1e-3)
276 {
277 nkp++;
278 ++kp;
279 }
280 kpoints.resize(nkp);
281 app_log() << " new spherical size = " << kpoints.size() << std::endl;
282 kp_end = kpoints.end();
283 // count degeneracies
284 kp = kpoints.begin();
285 while (kp != kp_end)
286 {
287 T k2 = kp->k2;
288 kptmp = kp;
289 int g = 1;
290 ++kptmp;
291 // look ahead to count
292 while (kptmp != kp_end && std::abs(kptmp->k2 - k2) < tol)
293 {
294 g++;
295 ++kptmp;
296 }
297 kp->g = g;
298 // run over degenerate states to assign
299 for (int n = 0; n < g - 1; ++n)
300 (++kp)->g = g;
301 ++kp;
302 }
303 //app_log()<<"create_kpoints"<< std::endl;
304 //app_log()<<" nkpoints = "<<nkpoints<< std::endl;
305 //app_log()<<" kpoints"<< std::endl;
306 //for(kp=kpoints.begin();kp!=kp_end;++kp)
307 // app_log()<<" "<<kp->k2<<" "<<kp->g<<" "<<kp->k<< std::endl;
308 //APP_ABORT("end create_kpoints");
309 }
310
311
312 void createGrid(int nc, int nkpts, const PosType& twistAngle)
313 {
314 twist = twistAngle;
315 create_kpoints(nc, twistAngle);
316 kpoints_t& kpoints = *kpoints_grid;
317 if (nkpts > kpoints.size())
318 APP_ABORT("HEGGrid::createGrid requested more kpoints than created");
319 kpt.resize(nkpts);
320 mk2.resize(nkpts);
321 deg.resize(nkpts);
322 for (int i = 0; i < nkpts; ++i)
323 {
324 const kpdata_t& kp = kpoints[i];
325 kpt[i] = kp.k;
326 mk2[i] = -kp.k2;
327 deg[i] = kp.g;
328 }
329 app_log() << "List of kpoints with twist = " << twistAngle << std::endl;
330 for (int ik = 0; ik < kpt.size(); ik++)
331 app_log() << ik << " " << kpt[ik] << " " << -mk2[ik] << std::endl;
332 }
333
334
335 void createGrid(const std::vector<int>& states, T tol = 1e-6) { createGrid(states, twist, tol); }
336
337
338 void createGrid(const std::vector<int>& states, const PosType& twistAngle, T tol = 1e-6)
339 {
340 int smax = 0;
341 for (int i = 0; i < states.size(); ++i)
342 smax = std::max(smax, states[i]);
343 smax++;
344 create_kpoints(get_nc(smax), twistAngle, tol);
345 kpoints_t& kpoints = *kpoints_grid;
346 if (smax > kpoints.size())
347 APP_ABORT("HEGGrid::createGrid(states) requested more kpoints than created");
348 int nkpts = states.size();
349 kpt.resize(nkpts);
350 mk2.resize(nkpts);
351 deg.resize(nkpts);
352 for (int i = 0; i < states.size(); ++i)
353 {
354 const kpdata_t& kp = kpoints[states[i]];
355 kpt[i] = kp.k;
356 mk2[i] = -kp.k2;
357 deg[i] = kp.g;
358 }
359 }
360 };
361
362
363 //two-d specialization
364 template<class T>
365 struct HEGGrid<T, 2>
366 {
367 typedef CrystalLattice<T, 2> PL_t;
368 typedef typename PL_t::SingleParticlePos_t PosType;
369
370 ///number of kpoints of a half sphere excluding gamma
371 int NumKptsHalf;
372 ///maxmim ksq
373 T MaxKsq;
374 PL_t& Lattice;
375 std::map<int, std::vector<PosType>*> rs;
376 std::vector<PosType> kpt;
377 std::vector<T> mk2;
378 std::vector<int> deg;
379 std::vector<int> n_within_shell;
380 PosType twist;
381
382 typedef kpdata<T, 2> kpdata_t;
383 typedef std::vector<kpdata_t> kpoints_t;
384
385 kpoints_t* kpoints_grid;
386
387 HEGGrid(PL_t& lat) : Lattice(lat), kpoints_grid(0)
388 {
389 n_within_shell.resize(220);
390 //fill this in
391 n_within_shell[0] = 1;
392 n_within_shell[1] = 5;
393 n_within_shell[2] = 9;
394 n_within_shell[3] = 13;
395 n_within_shell[4] = 21;
396 n_within_shell[5] = 25;
397 n_within_shell[6] = 29;
398 n_within_shell[7] = 37;
399 n_within_shell[8] = 45;
400 n_within_shell[9] = 49;
401 n_within_shell[10] = 57;
402 n_within_shell[11] = 61;
403 n_within_shell[12] = 69;
404 n_within_shell[13] = 81;
405 n_within_shell[14] = 89;
406 n_within_shell[15] = 97;
407 n_within_shell[16] = 101;
408 n_within_shell[17] = 109;
409 n_within_shell[18] = 113;
410 n_within_shell[19] = 121;
411 n_within_shell[20] = 129;
412 n_within_shell[21] = 137;
413 n_within_shell[22] = 145;
414 n_within_shell[23] = 149;
415 n_within_shell[24] = 161;
416 n_within_shell[25] = 169;
417 n_within_shell[26] = 177;
418 n_within_shell[27] = 185;
419 n_within_shell[28] = 193;
420 n_within_shell[29] = 197;
421 n_within_shell[30] = 213;
422 n_within_shell[31] = 221;
423 n_within_shell[32] = 225;
424 n_within_shell[33] = 233;
425 n_within_shell[34] = 241;
426 n_within_shell[35] = 249;
427 n_within_shell[36] = 253;
428 n_within_shell[37] = 261;
429 n_within_shell[38] = 277;
430 n_within_shell[39] = 285;
431 n_within_shell[40] = 293;
432 n_within_shell[41] = 301;
433 n_within_shell[42] = 305;
434 n_within_shell[43] = 317;
435 n_within_shell[44] = 325;
436 n_within_shell[45] = 333;
437 n_within_shell[46] = 341;
438 n_within_shell[47] = 349;
439 n_within_shell[48] = 357;
440 n_within_shell[49] = 365;
441 n_within_shell[50] = 373;
442 n_within_shell[51] = 377;
443 n_within_shell[52] = 385;
444 n_within_shell[53] = 401;
445 n_within_shell[54] = 405;
446 n_within_shell[55] = 421;
447 n_within_shell[56] = 429;
448 n_within_shell[57] = 437;
449 n_within_shell[58] = 441;
450 n_within_shell[59] = 457;
451 n_within_shell[60] = 465;
452 n_within_shell[61] = 473;
453 n_within_shell[62] = 481;
454 n_within_shell[63] = 489;
455 n_within_shell[64] = 497;
456 n_within_shell[65] = 505;
457 n_within_shell[66] = 509;
458 n_within_shell[67] = 517;
459 n_within_shell[68] = 529;
460 n_within_shell[69] = 545;
461 n_within_shell[70] = 553;
462 n_within_shell[71] = 561;
463 n_within_shell[72] = 569;
464 n_within_shell[73] = 577;
465 n_within_shell[74] = 593;
466 n_within_shell[75] = 601;
467 n_within_shell[76] = 609;
468 n_within_shell[77] = 613;
469 n_within_shell[78] = 621;
470 n_within_shell[79] = 633;
471 n_within_shell[80] = 641;
472 n_within_shell[81] = 657;
473 n_within_shell[82] = 665;
474 n_within_shell[83] = 673;
475 n_within_shell[84] = 681;
476 n_within_shell[85] = 697;
477 n_within_shell[86] = 709;
478 n_within_shell[87] = 717;
479 n_within_shell[88] = 725;
480 n_within_shell[89] = 733;
481 n_within_shell[90] = 741;
482 n_within_shell[91] = 749;
483 n_within_shell[92] = 757;
484 n_within_shell[93] = 761;
485 n_within_shell[94] = 769;
486 n_within_shell[95] = 777;
487 n_within_shell[96] = 793;
488 n_within_shell[97] = 797;
489 n_within_shell[98] = 805;
490 n_within_shell[99] = 821;
491 n_within_shell[100] = 829;
492 n_within_shell[101] = 845;
493 n_within_shell[102] = 853;
494 n_within_shell[103] = 861;
495 n_within_shell[104] = 869;
496 n_within_shell[105] = 877;
497 n_within_shell[106] = 885;
498 n_within_shell[107] = 889;
499 n_within_shell[108] = 901;
500 n_within_shell[109] = 917;
501 n_within_shell[110] = 925;
502 n_within_shell[111] = 933;
503 n_within_shell[112] = 941;
504 n_within_shell[113] = 949;
505 n_within_shell[114] = 965;
506 n_within_shell[115] = 973;
507 n_within_shell[116] = 981;
508 n_within_shell[117] = 989;
509 n_within_shell[118] = 997;
510 n_within_shell[119] = 1005;
511 n_within_shell[120] = 1009;
512 n_within_shell[121] = 1033;
513 n_within_shell[122] = 1041;
514 n_within_shell[123] = 1049;
515 n_within_shell[124] = 1057;
516 n_within_shell[125] = 1069;
517 n_within_shell[126] = 1085;
518 n_within_shell[127] = 1093;
519 n_within_shell[128] = 1101;
520 n_within_shell[129] = 1109;
521 n_within_shell[130] = 1117;
522 n_within_shell[131] = 1125;
523 n_within_shell[132] = 1129;
524 n_within_shell[133] = 1137;
525 n_within_shell[134] = 1153;
526 n_within_shell[135] = 1161;
527 n_within_shell[136] = 1177;
528 n_within_shell[137] = 1185;
529 n_within_shell[138] = 1201;
530 n_within_shell[139] = 1209;
531 n_within_shell[140] = 1217;
532 n_within_shell[141] = 1225;
533 n_within_shell[142] = 1229;
534 n_within_shell[143] = 1237;
535 n_within_shell[144] = 1245;
536 n_within_shell[145] = 1257;
537 n_within_shell[146] = 1265;
538 n_within_shell[147] = 1273;
539 n_within_shell[148] = 1281;
540 n_within_shell[149] = 1289;
541 n_within_shell[150] = 1305;
542 n_within_shell[151] = 1313;
543 n_within_shell[152] = 1321;
544 n_within_shell[153] = 1329;
545 n_within_shell[154] = 1353;
546 n_within_shell[155] = 1361;
547 n_within_shell[156] = 1369;
548 n_within_shell[157] = 1373;
549 n_within_shell[158] = 1389;
550 n_within_shell[159] = 1405;
551 n_within_shell[160] = 1413;
552 n_within_shell[161] = 1425;
553 n_within_shell[162] = 1433;
554 n_within_shell[163] = 1441;
555 n_within_shell[164] = 1449;
556 n_within_shell[165] = 1457;
557 n_within_shell[166] = 1465;
558 n_within_shell[167] = 1473;
559 n_within_shell[168] = 1481;
560 n_within_shell[169] = 1489;
561 n_within_shell[170] = 1505;
562 n_within_shell[171] = 1513;
563 n_within_shell[172] = 1517;
564 n_within_shell[173] = 1533;
565 n_within_shell[174] = 1541;
566 n_within_shell[175] = 1549;
567 n_within_shell[176] = 1565;
568 n_within_shell[177] = 1581;
569 n_within_shell[178] = 1597;
570 n_within_shell[179] = 1605;
571 n_within_shell[180] = 1609;
572 n_within_shell[181] = 1617;
573 n_within_shell[182] = 1633;
574 n_within_shell[183] = 1641;
575 n_within_shell[184] = 1649;
576 n_within_shell[185] = 1653;
577 n_within_shell[186] = 1669;
578 n_within_shell[187] = 1685;
579 n_within_shell[188] = 1693;
580 n_within_shell[189] = 1701;
581 n_within_shell[190] = 1709;
582 n_within_shell[191] = 1725;
583 n_within_shell[192] = 1733;
584 n_within_shell[193] = 1741;
585 n_within_shell[194] = 1749;
586 n_within_shell[195] = 1757;
587 n_within_shell[196] = 1765;
588 n_within_shell[197] = 1781;
589 n_within_shell[198] = 1789;
590 n_within_shell[199] = 1793;
591 n_within_shell[200] = 1801;
592 n_within_shell[201] = 1813;
593 n_within_shell[202] = 1829;
594 n_within_shell[203] = 1837;
595 n_within_shell[204] = 1853;
596 n_within_shell[205] = 1861;
597 n_within_shell[206] = 1869;
598 n_within_shell[207] = 1877;
599 n_within_shell[208] = 1885;
600 n_within_shell[209] = 1893;
601 n_within_shell[210] = 1901;
602 n_within_shell[211] = 1917;
603 n_within_shell[212] = 1925;
604 n_within_shell[213] = 1933;
605 n_within_shell[214] = 1941;
606 n_within_shell[215] = 1961;
607 n_within_shell[216] = 1969;
608 n_within_shell[217] = 1977;
609 n_within_shell[218] = 1993;
610 n_within_shell[219] = 2001;
611 }
612
613 ~HEGGrid()
614 {
615 typename std::map<int, std::vector<PosType>*>::iterator it(rs.begin()), it_end(rs.end());
616 while (it != it_end)
617 {
618 delete (*it).second;
619 ++it;
620 }
621 }
622
623 /** return the estimated number of grid in each direction */
624 inline int getNC(int nup) const { return static_cast<int>(std::pow(static_cast<T>(nup), 1.0 / 2)) / 2 + 1; }
625
626 //return the number of k-points upto nsh-shell
627 inline int getNumberOfKpoints(int nsh) const
628 {
629 if (nsh < n_within_shell.size())
630 return n_within_shell[nsh];
631 else
632 return -1;
633 }
634
635 //return the shell index for nkpt k-points
636 inline int getShellIndex(int nkpt) const
637 {
638 std::vector<int>::const_iterator loc = std::upper_bound(n_within_shell.begin(), n_within_shell.end(), nkpt);
639 if (loc < n_within_shell.end())
640 return loc - n_within_shell.begin() - 1;
641 else
642 return getNC(nkpt);
643 }
644
645 inline int getShellFromStates(int nst)
646 {
647 for (int i = 0; i < n_within_shell.size(); i++)
648 if (n_within_shell[i] == nst)
649 return i;
650 return -1;
651 }
652
653 inline T getCellLength(int nptcl, T rs_in) const { return std::sqrt(M_PI * nptcl) * rs_in; }
654
655 void sortGrid(int nc)
656 {
657 int first_ix2, first_ix3;
658 for (int ix1 = 0; ix1 <= nc; ix1++)
659 {
660 if (ix1 == 0)
661 first_ix2 = 1;
662 else
663 first_ix2 = -nc;
664 for (int ix2 = first_ix2; ix2 <= nc; ix2++)
665 {
666 int ih = ix1 * ix1 + ix2 * ix2;
667 typename std::map<int, std::vector<PosType>*>::iterator it = rs.find(ih);
668 if (it == rs.end())
669 {
670 std::vector<PosType>* ns = new std::vector<PosType>;
671 ns->push_back(PosType(ix1, ix2));
672 rs[ih] = ns;
673 }
674 else
675 {
676 (*it).second->push_back(PosType(ix1, ix2));
677 }
678 }
679 }
680 }
681
682 void createGrid(int nc, int nkpts)
683 {
684 if (rs.empty())
685 sortGrid(nc);
686 NumKptsHalf = nkpts;
687 kpt.resize(nkpts);
688 mk2.resize(nkpts);
689 int ikpt = 0;
690 //int checkNum=0;
691 //int ke=0;
692 MaxKsq = 0.0;
693 int rsbin = 0;
694 typename std::map<int, std::vector<PosType>*>::iterator rs_it(rs.begin()), rs_end(rs.end());
695 while (ikpt < nkpts && rs_it != rs_end)
696 {
697 typename std::vector<PosType>::iterator ns_it((*rs_it).second->begin()), ns_end((*rs_it).second->end());
698 T minus_ksq = -Lattice.ksq(*ns_it);
699 while (ikpt < nkpts && ns_it != ns_end)
700 {
701 kpt[ikpt] = Lattice.k_cart(*ns_it);
702 mk2[ikpt] = minus_ksq;
703 ++ikpt;
704 ++ns_it;
705 }
706 ++rsbin;
707 ++rs_it;
708 }
709 MaxKsq = Lattice.ksq(*((*rs_it).second->begin()));
710 app_log() << "List of kpoints (half-sphere) " << std::endl;
711 for (int ik = 0; ik < kpt.size(); ik++)
712 {
713 app_log() << ik << " " << kpt[ik] << " " << mk2[ik] << std::endl;
714 }
715 }
716
717
718 void createGrid(const PosType& twistAngle)
719 {
720 twist = twistAngle;
721 //unfold and add gamma
722 int nkpts = 2 * NumKptsHalf + 1;
723 kpt.resize(nkpts);
724 mk2.resize(nkpts);
725 app_log() << "Check this " << NumKptsHalf << " " << nkpts << std::endl;
726 abort();
727 //add gamma
728 int ikpt = 0;
729 kpt[ikpt] = Lattice.k_cart(twistAngle);
730 mk2[ikpt] = -Lattice.ksq(twistAngle);
731 ++ikpt;
732 typename std::map<int, std::vector<PosType>*>::iterator rs_it(rs.begin()), rs_end(rs.end());
733 while (ikpt < nkpts && rs_it != rs_end)
734 {
735 typename std::vector<PosType>::iterator ns_it((*rs_it).second->begin()), ns_end((*rs_it).second->end());
736 while (ikpt < nkpts && ns_it != ns_end)
737 {
738 //add twist+k
739 PosType k0(twistAngle + (*ns_it));
740 T ksq = Lattice.ksq(k0);
741 kpt[ikpt] = Lattice.k_cart(k0);
742 mk2[ikpt] = -ksq;
743 ++ikpt;
744 //add twist-k
745 k0 = twistAngle - (*ns_it);
746 ksq = Lattice.ksq(k0);
747 kpt[ikpt] = Lattice.k_cart(k0);
748 mk2[ikpt] = -ksq;
749 ++ikpt;
750 ++ns_it;
751 }
752 ++rs_it;
753 }
754 app_log() << "List of kpoints with twist = " << twistAngle << std::endl;
755 for (int ik = 0; ik < kpt.size(); ik++)
756 {
757 app_log() << ik << " " << kpt[ik] << " " << -mk2[ik] << std::endl;
758 }
759 }
760
761
762 void createGrid(const std::vector<int>& states, T tol = 1e-6)
763 {
764 APP_ABORT("HEGGrid::createGrid(states) has not been implemented");
765 }
766 };
767 } // namespace qmcplusplus
768
769 #endif
770