VMPS++
Loading...
Searching...
No Matches
ParamCollection.h
Go to the documentation of this file.
1#ifndef PARAMCOLLECTION
2#define PARAMCOLLECTION
3
4#include <set>
5
6#include <Eigen/Dense>
7using namespace Eigen;
8
9#include "ParamHandler.h"
10#include "CuthillMcKeeCompressor.h" // from ALGS
11#include <boost/rational.hpp>
12#include "termcolor.hpp"
13#include "DmrgTypedefs.h" // for SUB_LATTICE
14
15ArrayXXd create_1D_OBC (size_t L, double lambda1=1., double lambda2=0.)
16{
17 ArrayXXd res(L,L); res.setZero();
18
19 res.matrix().diagonal<1>().setConstant(lambda1);
20 res.matrix().diagonal<-1>().setConstant(lambda1);
21
22 res.matrix().diagonal<2>().setConstant(lambda2);
23 res.matrix().diagonal<-2>().setConstant(lambda2);
24
25 return res;
26}
27
28// Simple generation of periodic boundary conditions with NN and NNN coupling
29// If COMPRESSED=true, the ring is flattened so that adjacent sites are as close together as possible
30ArrayXXd create_1D_PBC (size_t L, double lambda1=1., double lambda2=0., bool COMPRESSED=false)
31{
32 ArrayXXd res(L,L);
33 res.setZero();
34
35 res = create_1D_OBC(L,lambda1,lambda2);
36
37 res(0,L-1) = lambda1;
38 res(L-1,0) = lambda1;
39
40 res(0,L-2) = lambda2;
41 res(L-2,0) = lambda2;
42 res(1,L-1) = lambda2;
43 res(L-1,1) = lambda2;
44
45 if (COMPRESSED and lambda2 == 0.)
46 {
47 res.setZero();
48
49 res(0,1) = lambda1; res(1,0) = lambda1;
50 res(L-2,L-1) = lambda1; res(L-1,L-2) = lambda1;
51 for (size_t l=0; l<L-2; l++)
52 {
53 res(l,l+2) = lambda1;
54 res(l+2,l) = lambda1;
55 }
56 }
57 else if (COMPRESSED and lambda2 != 0.)
58 {
59 auto res_ = compress_CuthillMcKee(res,true);
60 res = res_;
61 }
62 return res;
63}
64
65ArrayXXd extend_to_thermal (const ArrayXXd &tFull, double factor)
66{
67 int L = tFull.rows();
68 ArrayXXd res(2*L,2*L); res.setZero();
69
70 for (int i=0; i<L; ++i)
71 for (int j=0; j<L; ++j)
72 {
73 if (abs(tFull(i,j)) > 1e-10)
74 {
75 res(2*i,2*j) = tFull(i,j);
76 res(2*i+1,2*j+1) = factor*tFull(i,j);
77 }
78 }
79 return res;
80}
81
82ArrayXXd create_1D_AB (size_t L, double lambda1A=1., double lambda1B=1., double lambda2A=0., double lambda2B=0.)
83{
84 ArrayXXd res(L,L);
85 res.setZero();
86
87 for (int i=0; i<L; i+=2)
88 {
89 if (i+1<L) res(i, i+1) = lambda1A;
90 if (i+2<L) res(i+1,i+2) = lambda1B;
91 }
92
93 for (int i=0; i<L; i+=2)
94 {
95 if (i+2<L) res(i ,i+2) = lambda2A;
96 if (i+3<L) res(i+1,i+3) = lambda2B;
97 }
98
99 res += res.transpose().eval();
100
101 return res;
102}
103
104ArrayXXd create_1D_PBC_AB (size_t L, double lambda1A=1., double lambda1B=1., double lambda2A=0., double lambda2B=0., bool COMPRESSED=true)
105{
106 ArrayXXd res(L,L);
107 res.setZero();
108
109 for (int i=0; i<L; i+=2)
110 {
111 res(i, (i+1)%L) = lambda1A;
112 if (i+1<L) res(i+1,(i+2)%L) = lambda1B;
113 }
114
115 for (int i=0; i<L; i+=2)
116 {
117 res(i, (i+2)%L) = lambda2A;
118 if (i+1<L) res(i+1,(i+3)%L) = lambda2B;
119 }
120
121 res += res.transpose().eval();
122
123// cout << "res=" << endl << res << endl;
124
125 if (COMPRESSED)
126 {
127 auto res_ = compress_CuthillMcKee(res,true);
128 res = res_;
129 }
130
131 return res;
132}
133
134ArrayXXd hopping_square (int Lx, int Ly, bool PBCx=false, bool PBCy=true, double lambda=1.)
135{
136 ArrayXXd res(Lx*Ly,Lx*Ly); res.setZero();
137 // vertical
138 for (int x=0; x<Lx; ++x)
139 {
140 int i0 = Ly*x;
141 for (int i=i0; i<=i0+Ly-2; ++i)
142 {
143 res(i,i+1) = lambda;
144 }
145 // y-periodic part for all x:
146 if (PBCy) res(i0,i0+Ly-1) = lambda;
147 }
148 // horizontal
149 for (int y=0; y<Ly; ++y)
150 {
151 for (int x=0; x<Lx-1; ++x)
152 {
153 int i = Ly*x+y;
154 res(i,i+Ly) = lambda;
155 }
156 }
157 // x-periodic part for all y:
158 if (PBCx)
159 {
160 for (int i=0; i<Ly; ++i)
161 {
162 res(i,Ly*Lx-Ly+i) = lambda;
163 }
164 }
165 res += res.transpose().eval();
166 return res;
167}
168
169// split into batches of Ly for A, and Ly/2 for B (depleted sites)
170void split_kagomeYC_AB (int Lx, int Ly, std::vector<std::vector<int>> &A, std::vector<std::vector<int>> &B)
171{
172 int Nevn = Lx/2;
173 int Nodd = Lx/2;
174 int L = Ly*Nevn+Ly/2*Nodd;
175
176 vector<int> input(L);
177 for (int j=0; j<L; ++j) input[j] = j;
178
179 int i = 0;
180
181 while (i+Ly <= L)
182 {
183 std::vector<int> batchA;
184 std::vector<int> batchB;
185
186 for (int j=0; j<Ly; j++)
187 {
188 batchA.push_back(input[i+j]);
189 }
190
191 A.push_back(batchA);
192 i += Ly;
193
194 if (i + Ly/2 <= L)
195 {
196 for (int j=0; j<Ly/2; j++)
197 {
198 batchB.push_back(input[i+j]);
199 }
200 B.push_back(batchB);
201 i += Ly/2;
202 }
203 }
204}
205
206// split into batches of Ly/2 for A (depleted sites), and Ly for B
207void split_kagomeYC_BAB (int Lx, int Ly, std::vector<std::vector<int>> &A, std::vector<std::vector<int>> &B)
208{
209 int Nevn = Lx/2+1;
210 int Nodd = Lx/2;
211 int L = Ly/2*Nevn+Ly*Nodd;
212
213 std::vector<int> input(L);
214 for (int j=0; j<L; ++j) input[j] = j;
215
216 int i = 0;
217
218 while (i + Ly/2 <= L)
219 {
220 std::vector<int> batchA;
221
222 for (int j=0; j<Ly/2; j++)
223 {
224 batchA.push_back(input[i+j]);
225 }
226
227 A.push_back(batchA);
228 i += Ly/2;
229
230 if (i + Ly <= L)
231 {
232 std::vector<int> batchB;
233
234 for (int j=0; j<Ly; j++)
235 {
236 batchB.push_back(input[i+j]);
237 }
238
239 B.push_back(batchB);
240 i += Ly;
241 }
242 }
243}
244
245int find_x_kagomeYC (int index, int L, const std::vector<std::vector<int>> &A, const std::vector<std::vector<int>> &B)
246{
247 //int Nevn = Lx/2;
248 //int Nodd = Lx/2;
249 //int L = Ly*Nevn+Ly/2*Nodd;
250
251 MatrixXi res(L,3);
252
253 bool EVN = false;
254 bool ODD = false;
255 int i0 = -1;
256
257 for (int i=0; i<A.size(); ++i)
258 {
259 auto it = find(A[i].begin(), A[i].end(), index);
260 if (it != A[i].end())
261 {
262 EVN = true;
263 i0 = i;
264 break;
265 }
266 }
267
268 if (!EVN)
269 {
270 for (int i=0; i<B.size(); ++i)
271 {
272 auto it = find(B[i].begin(), B[i].end(), index);
273 if (it != B[i].end())
274 {
275 ODD = true;
276 i0 = i;
277 break;
278 }
279 }
280 }
281
282 assert(EVN or ODD);
283
284 return (EVN)? 2*i0 : 2*i0+1;
285}
286
287ArrayXXd hopping_kagomeYC_AB (int Lx, int Ly, bool PBCx=false, bool PBCy=true, double lambda=1.)
288{
289 assert(PBCx==false and PBCy==true);
290
291 int Nevn = Lx/2;
292 int Nodd = Lx/2;
293 int L = Ly*Nevn+Ly/2*Nodd;
294
295 ArrayXXd res(L,L); res.setZero();
296
297 vector<vector<int>> i_evn;
298 vector<vector<int>> i_odd;
299
300 split_kagomeYC_BAB(Lx, Ly, i_evn, i_odd);
301
302// lout << "EVN:" << endl;
303//
304// for (int i=0; i<i_evn.size(); ++i)
305// {
306// for (int j=0; j<i_evn[i].size(); ++j)
307// {
308// lout << i_evn[i][j] << endl;
309// }
310// lout << "----" << endl;
311// }
312//
313// lout << "ODD:" << endl;
314//
315// for (int i=0; i<i_odd.size(); ++i)
316// {
317// for (int j=0; j<i_odd[i].size(); ++j)
318// {
319// lout << i_odd[i][j] << endl;
320// }
321// lout << "----" << endl;
322// }
323
324 // vertical
325 for (int x=0; x<i_evn.size(); ++x)
326 for (int i=0; i<i_evn[x].size(); ++i)
327 {
328 int k = i_evn[x][i];
329 int l = i_evn[x][(i+1)%Ly];
330 //lout << "vertical evn bond: " << k << ", " << l << endl;
331 res(min(k,l),max(k,l)) = lambda;
332 }
333
334 // horizontal, even x
335 for (int x=0; x<i_evn.size(); ++x)
336 for (int i=1; i<i_evn[x].size(); i+=2)
337 {
338 int k = i_evn[x][i];
339 int l = i_odd[x][(i-1)/2];
340 //lout << "horizontal evn bond: " << k << ", " << l << endl;
341 res(min(k,l),max(k,l)) = lambda;
342 }
343
344 // horizontal, odd x
345 for (int x=0; x<i_odd.size()-1; ++x)
346 for (int i=0; i<i_odd[x].size(); i+=1)
347 {
348 int k = i_odd[x][i];
349 int l = i_evn[x+1][2*i+1];
350 //lout << "horizontal odd bond: " << k << ", " << l << endl;
351 res(min(k,l),max(k,l)) = lambda;
352 }
353
354 // diagonal, even x
355 for (int x=0; x<i_evn.size(); ++x)
356 for (int i=0; i<i_evn[x].size(); i+=2)
357 {
358 int k = i_evn[x][i];
359 int l = i_odd[x][i/2];
360 //lout << "diagonal evn bond: " << k << ", " << l << endl;
361 res(min(k,l),max(k,l)) = lambda;
362 }
363
364 // diagonal, odd x
365 for (int x=0; x<i_odd.size()-1; ++x)
366 for (int i=0; i<i_odd[x].size(); ++i)
367 {
368 int k = i_odd[x][i];
369 int l = i_evn[x+1][(2*i+2)%Ly];
370 //lout << "diagonal odd bond: " << k << ", " << l << endl;
371 res(min(k,l),max(k,l)) = lambda;
372 }
373
374 res += res.transpose().eval();
375
376 //lout << res << endl;
377
378 return res;
379}
380
381ArrayXXd hopping_kagomeYC_BAB (int Lx, int Ly, bool PBCx=false, bool PBCy=true, double lambda=1., bool VERBOSE=false)
382{
383 assert(PBCx==false and PBCy==true);
384
385 int Nevn = Lx/2+1;
386 int Nodd = Lx/2;
387 int L = Ly/2*Nevn+Ly*Nodd;
388
389 ArrayXXd res(L,L); res.setZero();
390
391 vector<vector<int>> i_evn;
392 vector<vector<int>> i_odd;
393
394 split_kagomeYC_BAB(Lx, Ly, i_evn, i_odd);
395
396 if (VERBOSE)
397 {
398 lout << "EVN:" << endl;
399 for (int i=0; i<i_evn.size(); ++i)
400 {
401 for (int j=0; j<i_evn[i].size(); ++j)
402 {
403 lout << i_evn[i][j] << endl;
404 }
405 lout << "----" << endl;
406 }
407
408 lout << "ODD:" << endl;
409 for (int i=0; i<i_odd.size(); ++i)
410 {
411 for (int j=0; j<i_odd[i].size(); ++j)
412 {
413 lout << i_odd[i][j] << endl;
414 }
415 lout << "----" << endl;
416 }
417 }
418
419 // vertical, odd x
420 for (int x=0; x<i_odd.size(); ++x)
421 for (int i=0; i<i_odd[x].size(); ++i)
422 {
423 int k = i_odd[x][i];
424 int l = i_odd[x][(i+1)%Ly];
425 if (VERBOSE) lout << "vertical odd bond: " << min(k,l) << ", " << max(k,l) << endl;
426 res(min(k,l),max(k,l)) = lambda;
427 }
428
429 // horizontal, evn x
430 for (int x=0; x<i_evn.size()-1; ++x)
431 for (int i=0; i<i_evn[x].size(); i+=1)
432 {
433 int k = i_evn[x][i];
434 int l = i_odd[x][2*i+1];
435 if (VERBOSE) lout << "horizontal evn bond: " << min(k,l) << ", " << max(k,l) << endl;
436 res(min(k,l),max(k,l)) = lambda;
437 }
438
439 // horizontal, odd x
440 for (int x=0; x<i_odd.size(); ++x)
441 for (int i=1; i<i_odd[x].size(); i+=2)
442 {
443 int k = i_odd[x][i];
444 int l = i_evn[x+1][(i-1)/2];
445 if (VERBOSE) lout << "horizontal odd bond: " << min(k,l) << ", " << max(k,l) << endl;
446 res(min(k,l),max(k,l)) = lambda;
447 }
448
449 // diagonal, evn x
450 for (int x=0; x<i_evn.size()-1; ++x)
451 for (int i=0; i<i_evn[x].size(); ++i)
452 {
453 int k = i_evn[x][i];
454 int l = i_odd[x][(2*i+2)%Ly];
455 if (VERBOSE) lout << "diagonal evn bond: " << min(k,l) << ", " << max(k,l) << endl;
456 res(min(k,l),max(k,l)) = lambda;
457 }
458
459 // diagonal, odd x
460 for (int x=0; x<i_odd.size(); ++x)
461 for (int i=0; i<i_odd[x].size(); i+=2)
462 {
463 int k = i_odd[x][i];
464 int l = i_evn[x+1][i/2];
465 if (VERBOSE) lout << "diagonal odd bond: " << min(k,l) << ", " << max(k,l) << endl;
466 res(min(k,l),max(k,l)) = lambda;
467 }
468
469 res += res.transpose().eval();
470
471 //lout << res << endl;
472
473 return res;
474}
475
476void split_kagomeXC (int Lx, int Ly, std::vector<std::vector<int>> &A, std::vector<std::vector<int>> &B, std::vector<std::vector<int>> &C, std::vector<std::vector<int>> &D)
477{
478 int N = Ly/4;
479 int L = N*Lx + N*(Lx/2+1) + N*Lx + N*Lx/2;
480
481 std::vector<int> input(L);
482 for (int j=0; j<L; ++j) input[j] = j;
483
484 int currentIndex = 0;
485 while(currentIndex < L)
486 {
487 // Batch of Lx into A
488 vector<int> subA;
489 for(int i=currentIndex; i<currentIndex+Lx && i<L; i++) {
490 subA.push_back(input[i]);
491 }
492 A.push_back(subA);
493 currentIndex += Lx;
494 if(currentIndex >= L) break;
495
496 // Batch of Lx/2+1 into B
497 vector<int> subB;
498 for(int i=currentIndex; i<currentIndex+Lx/2+1 && i<L; i++) {
499 subB.push_back(input[i]);
500 }
501 B.push_back(subB);
502 currentIndex += Lx/2+1;
503 if(currentIndex >= L) break;
504
505 // Another batch of Lx into C
506 vector<int> subC;
507 for(int i=currentIndex; i<currentIndex+Lx && i<L; i++) {
508 subC.push_back(input[i]);
509 }
510 C.push_back(subC);
511 currentIndex += Lx;
512 if(currentIndex >= L) break;
513
514 // Batch of Lx/2 into D
515 vector<int> subD;
516 for(int i=currentIndex; i<currentIndex+Lx/2 && i<L; i++) {
517 subD.push_back(input[i]);
518 }
519 D.push_back(subD);
520 currentIndex += Lx/2;
521 }
522
523// for (int i=0; i<Ly/4; ++i)
524// {
525// lout << "A:" << endl;
526// for (int j=0; j<iA[i].size(); ++j)
527// {
528// lout << iA[i][j] << ", ";
529// }
530// lout << endl;
531// lout << "B:" << endl;
532// for (int j=0; j<iB[i].size(); ++j)
533// {
534// lout << iB[i][j] << ", ";
535// }
536// lout << endl;
537// lout << "C:" << endl;
538// for (int j=0; j<iC[i].size(); ++j)
539// {
540// lout << iC[i][j] << ", ";
541// }
542// lout << endl;
543// lout << "D:" << endl;
544// for (int j=0; j<iD[i].size(); ++j)
545// {
546// lout << iD[i][j] << ", ";
547// }
548// lout << endl;
549// }
550}
551
552ArrayXXd hopping_kagomeXC (int Lx, int Ly, bool PBCx=false, bool PBCy=true, double lambda=1., bool VERBOSE=false)
553{
554 int N = Ly/4;
555 int L = N*Lx + N*(Lx/2+1) + N*Lx + N*Lx/2;
556
557 ArrayXXd res(L,L); res.setZero();
558
559 vector<vector<int>> iA, iB, iC, iD;
560 split_kagomeXC(Lx, Ly, iA, iB, iC, iD);
561
562 for (int y=0; y<N; ++y)
563 {
564 if (VERBOSE) lout << "y=" << y << endl;
565
566 // horizontal bonds
567 for (int i=0; i<iA[y].size()-1; ++i)
568 {
569 int k = iA[y][i];
570 int l = iA[y][i+1];
571 res(min(k,l),max(k,l)) = lambda;
572 if (VERBOSE) lout << "A-A: " << min(k,l) << ", " << max(k,l) << endl;
573 }
574 for (int i=0; i<iC[y].size()-1; ++i)
575 {
576 int k = iC[y][i];
577 int l = iC[y][i+1];
578 res(min(k,l),max(k,l)) = lambda;
579 if (VERBOSE) lout << "C-C: " << min(k,l) << ", " << max(k,l) << endl;
580 }
581
582 // vertical bonds
583 for (int i=0; i<iA[y].size(); ++i)
584 {
585 if (i==0)
586 {
587 int k = iA[y][0];
588 int l = iB[y][0];
589 res(min(k,l),max(k,l)) = lambda;
590 if (VERBOSE) lout << "A-B first: " << min(k,l) << ", " << max(k,l) << endl;
591 }
592 else if (i==Lx-1)
593 {
594 int k = iA[y][Lx-1];
595 int l = iB[y][Lx/2];
596 res(min(k,l),max(k,l)) = lambda;
597 if (VERBOSE) lout << "A-B last: " << min(k,l) << ", " << max(k,l) << endl;
598 }
599 else
600 {
601 if (i%2==0)
602 {
603 int k = iA[y][i];
604 int l = iB[y][i/2];
605 res(min(k,l),max(k,l)) = lambda;
606 if (VERBOSE) lout << "A-B: " << min(k,l) << ", " << max(k,l) << endl;
607 }
608 else
609 {
610 int k = iA[y][i];
611 int l = iB[y][(i+1)/2];
612 res(min(k,l),max(k,l)) = lambda;
613 if (VERBOSE) lout << "A-B: " << min(k,l) << ", " << max(k,l) << endl;
614 }
615 }
616 }
617
618 for (int i=0; i<iB[y].size(); ++i)
619 {
620 if (i==0)
621 {
622 int k = iB[y][0];
623 int l = iC[y][0];
624 res(min(k,l),max(k,l)) = lambda;
625 if (VERBOSE) lout << "B-C first: " << min(k,l) << ", " << max(k,l) << endl;
626 }
627 else if (i==Lx/2)
628 {
629 int k = iB[y][Lx/2];
630 int l = iC[y][Lx-1];
631 res(min(k,l),max(k,l)) = lambda;
632 if (VERBOSE) lout << "B-C last: " << min(k,l) << ", " << max(k,l) << endl;
633 }
634 else
635 {
636 if (i%2==1)
637 {
638 int k = iB[y][i];
639 int l = iC[y][2*i];
640 res(min(k,l),max(k,l)) = lambda;
641 if (VERBOSE) lout << "B-C: " << min(k,l) << ", " << max(k,l) << endl;
642
643 k = iB[y][i];
644 l = iC[y][2*i-1];
645 res(min(k,l),max(k,l)) = lambda;
646 if (VERBOSE) lout << "B-C: " << min(k,l) << ", " << max(k,l) << endl;
647 }
648 else
649 {
650 int k = iB[y][i];
651 int l = iC[y][2*i];
652 res(min(k,l),max(k,l)) = lambda;
653 if (VERBOSE) lout << "B-C: " << min(k,l) << ", " << max(k,l) << endl;
654
655 k = iB[y][i];
656 l = iC[y][2*i-1];
657 res(min(k,l),max(k,l)) = lambda;
658 if (VERBOSE) lout << "B-C: " << min(k,l) << ", " << max(k,l) << endl;
659 }
660 }
661 }
662
663 for (int i=0; i<iD[y].size(); ++i)
664 {
665 int k = iC[y][2*i];
666 int l = iD[y][i];
667 res(min(k,l),max(k,l)) = lambda;
668 if (VERBOSE) lout << "C-D: " << min(k,l) << ", " << max(k,l) << endl;
669
670 k = iC[y][2*i+1];
671 l = iD[y][i];
672 res(min(k,l),max(k,l)) = lambda;
673 if (VERBOSE) lout << "C-D: " << min(k,l) << ", " << max(k,l) << endl;
674 }
675
676 // includes periodic y-bonds
677 for (int i=0; i<iD[y].size(); ++i)
678 {
679 int k = iA[(y+1)%iA.size()][2*i];
680 int l = iD[y][i];
681 res(min(k,l),max(k,l)) = lambda;
682 if (VERBOSE) lout << "D-A: " << min(k,l) << ", " << max(k,l) << endl;
683
684 k = iA[(y+1)%iA.size()][2*i+1];
685 l = iD[y][i];
686 res(min(k,l),max(k,l)) = lambda;
687 if (VERBOSE) lout << "D-A: " << min(k,l) << ", " << max(k,l) << endl;
688 }
689 }
690
691 res += res.transpose().eval();
692 return res;
693}
694
695ArrayXXd hopping_triangularYC (int Lx, int Ly, bool PBCx=false, bool PBCy=true, double lambda=1.)
696{
697 ArrayXXd res = hopping_square(Lx,Ly,PBCx,PBCy,lambda);
698 res.matrix().triangularView<Eigen::Upper>().setZero();
699 // diagonal
700 for (int x=0; x<Lx-1; ++x)
701 {
702 for (int y=1; y<Ly; ++y)
703 {
704 int i = y+Ly*x;
705 res(i,i+Ly-1) = lambda;
706 }
707 // y-periodic part for all x:
708 if (PBCy) res(Ly*x,Ly*x+2*Ly-1) = lambda;
709 }
710 // x-periodic part for all y:
711 if (PBCx)
712 {
713 for (int i=0; i<Ly-1; ++i)
714 {
715 res(i,Ly*Lx-Ly+i+1) = lambda;
716 }
717 res(Ly-1,Ly*Lx-Ly) = lambda;
718 }
719 res += res.transpose().eval();
720 return res;
721}
722
723void add_triangle (int i, int j, int k, ArrayXXd &target, double lambda=1.)
724{
725 assert(i<j and j<k);
726 target(i,j) = lambda;
727 target(j,k) = lambda;
728 target(i,k) = lambda;
729}
730
731// Leung, Elser PRB 47, 9 (1992)
732ArrayXXd hopping_kagome36d (double lambda=1.)
733{
734 ArrayXXd res(36,36);
735 res.setZero();
736
737 add_triangle(0,3,4,res,lambda);
738 add_triangle(1,5,6,res,lambda);
739 add_triangle(2,3,7,res,lambda);
740 add_triangle(4,5,8,res,lambda);
741 add_triangle(6,9,21,res,lambda);
742 add_triangle(7,11,12,res,lambda);
743 add_triangle(8,13,14,res,lambda);
744 add_triangle(9,15,16,res,lambda);
745 add_triangle(10,11,17,res,lambda);
746 add_triangle(12,13,18,res,lambda);
747 add_triangle(14,15,19,res,lambda);
748 add_triangle(16,20,31,res,lambda);
749 add_triangle(17,21,22,res,lambda);
750 add_triangle(18,23,24,res,lambda);
751 add_triangle(19,25,26,res,lambda);
752 add_triangle(2,20,27,res,lambda);
753 add_triangle(22,23,28,res,lambda);
754 add_triangle(24,25,29,res,lambda);
755 add_triangle(26,27,30,res,lambda);
756 add_triangle(28,31,32,res,lambda);
757 add_triangle(29,33,35,res,lambda);
758 add_triangle(10,30,35,res,lambda);
759 add_triangle(0,32,33,res,lambda);
760 add_triangle(1,34,35,res,lambda);
761
762 res += res.transpose().eval();
763
764 auto res_ = compress_CuthillMcKee(res,true);
765 res = res_;
766
767 return res;
768}
769
770ArrayXXd triangularFlake (int L, double lambda=1.)
771{
772 ArrayXXd res(L,L); res.setZero();
773
774 if (L == 28 or L == 36 or L == 45 or L == 55 or L == 66)
775 {
776 add_triangle(0,1,2,res,lambda);
777
778 add_triangle(1,3,4,res,lambda);
779 add_triangle(1,2,4,res,lambda);
780 add_triangle(2,4,5,res,lambda);
781
782 add_triangle(3,6,7,res,lambda);
783 add_triangle(3,4,7,res,lambda);
784 add_triangle(4,7,8,res,lambda);
785 add_triangle(4,5,8,res,lambda);
786 add_triangle(5,8,9,res,lambda);
787
788 add_triangle(6,10,11,res,lambda);
789 add_triangle(6,7,11,res,lambda);
790 add_triangle(7,11,12,res,lambda);
791 add_triangle(7,8,12,res,lambda);
792 add_triangle(8,12,13,res,lambda);
793 add_triangle(8,9,13,res,lambda);
794 add_triangle(9,13,14,res,lambda);
795
796 add_triangle(10,15,16,res,lambda);
797 add_triangle(10,11,16,res,lambda);
798 add_triangle(11,16,17,res,lambda);
799 add_triangle(11,12,17,res,lambda);
800 add_triangle(12,17,18,res,lambda);
801 add_triangle(12,13,18,res,lambda);
802 add_triangle(13,18,19,res,lambda);
803 add_triangle(13,14,19,res,lambda);
804 add_triangle(14,19,20,res,lambda);
805
806 add_triangle(15,21,22,res,lambda);
807 add_triangle(15,16,22,res,lambda);
808 add_triangle(16,22,23,res,lambda);
809 add_triangle(16,17,23,res,lambda);
810 add_triangle(17,23,24,res,lambda);
811 add_triangle(17,18,24,res,lambda);
812 add_triangle(18,24,25,res,lambda);
813 add_triangle(18,19,25,res,lambda);
814 add_triangle(19,25,26,res,lambda);
815 add_triangle(19,20,26,res,lambda);
816 add_triangle(20,26,27,res,lambda);
817
818 if (L >= 36)
819 {
820 add_triangle(21,28,29,res,lambda);
821 add_triangle(21,22,29,res,lambda);
822 add_triangle(22,29,30,res,lambda);
823 add_triangle(22,23,30,res,lambda);
824 add_triangle(23,30,31,res,lambda);
825 add_triangle(23,24,31,res,lambda);
826 add_triangle(24,31,32,res,lambda);
827 add_triangle(24,25,32,res,lambda);
828 add_triangle(25,32,33,res,lambda);
829 add_triangle(25,26,33,res,lambda);
830 add_triangle(26,33,34,res,lambda);
831 add_triangle(26,27,34,res,lambda);
832 add_triangle(27,34,35,res,lambda);
833 }
834
835 if (L >= 45)
836 {
837 add_triangle(28,36,37,res,lambda);
838 add_triangle(29,37,38,res,lambda);
839 add_triangle(30,38,39,res,lambda);
840 add_triangle(31,39,40,res,lambda);
841 add_triangle(32,40,41,res,lambda);
842 add_triangle(33,41,42,res,lambda);
843 add_triangle(34,42,43,res,lambda);
844 add_triangle(35,43,44,res,lambda);
845
846 add_triangle(28,29,37,res,lambda);
847 add_triangle(29,30,38,res,lambda);
848 add_triangle(30,31,39,res,lambda);
849 add_triangle(31,32,40,res,lambda);
850 add_triangle(32,33,41,res,lambda);
851 add_triangle(33,34,42,res,lambda);
852 add_triangle(34,35,43,res,lambda);
853
854 if (L >= 55)
855 {
856 add_triangle(36,45,46,res,lambda);
857 add_triangle(37,38,47,res,lambda);
858 add_triangle(38,39,48,res,lambda);
859 add_triangle(39,40,49,res,lambda);
860 add_triangle(40,41,50,res,lambda);
861 add_triangle(41,42,51,res,lambda);
862 add_triangle(42,43,52,res,lambda);
863 add_triangle(43,44,53,res,lambda);
864
865 add_triangle(36,45,46,res,lambda);
866 add_triangle(37,46,47,res,lambda);
867 add_triangle(28,36,37,res,lambda);
868 add_triangle(38,47,48,res,lambda);
869 add_triangle(39,48,49,res,lambda);
870 add_triangle(40,49,50,res,lambda);
871 add_triangle(41,50,51,res,lambda);
872 add_triangle(42,51,52,res,lambda);
873 add_triangle(43,52,53,res,lambda);
874 add_triangle(44,53,54,res,lambda);
875
876 if (L >= 66)
877 {
878 add_triangle(45,46,56,res,lambda);
879 add_triangle(46,47,57,res,lambda);
880 add_triangle(47,48,58,res,lambda);
881 add_triangle(48,49,59,res,lambda);
882 add_triangle(49,50,60,res,lambda);
883 add_triangle(50,51,61,res,lambda);
884 add_triangle(51,52,62,res,lambda);
885 add_triangle(52,53,63,res,lambda);
886 add_triangle(53,54,64,res,lambda);
887
888 add_triangle(45,55,56,res,lambda);
889 add_triangle(46,56,57,res,lambda);
890 add_triangle(47,57,58,res,lambda);
891 add_triangle(48,58,59,res,lambda);
892 add_triangle(49,59,60,res,lambda);
893 add_triangle(50,60,61,res,lambda);
894 add_triangle(51,61,62,res,lambda);
895 add_triangle(52,62,63,res,lambda);
896 add_triangle(53,63,64,res,lambda);
897 add_triangle(54,64,65,res,lambda);
898 }
899 }
900 }
901 }
902
903 res += res.transpose().eval();
904
905 auto res_ = compress_CuthillMcKee(res,true);
906 res = res_;
907
908 return res;
909}
910
911ArrayXXd hexagonalFlake (int L, double lambda=1.)
912{
913 ArrayXXd res(L,L); res.setZero();
914
915 if (L == 48)
916 {
917 add_triangle(0,4,5,res,lambda);
918 add_triangle(0,1,5,res,lambda);
919 add_triangle(1,5,6,res,lambda);
920 add_triangle(1,2,6,res,lambda);
921 add_triangle(2,6,7,res,lambda);
922 add_triangle(2,3,7,res,lambda);
923 add_triangle(3,7,8,res,lambda);
924
925 add_triangle(4,9,10,res,lambda);
926 add_triangle(4,5,10,res,lambda);
927 add_triangle(5,10,11,res,lambda);
928 add_triangle(5,6,11,res,lambda);
929 add_triangle(6,11,12,res,lambda);
930 add_triangle(6,7,12,res,lambda);
931 add_triangle(7,12,13,res,lambda);
932 add_triangle(7,8,13,res,lambda);
933 add_triangle(8,13,14,res,lambda);
934
935 add_triangle(9,15,16,res,lambda);
936 add_triangle(9,10,16,res,lambda);
937 add_triangle(10,16,17,res,lambda);
938 add_triangle(10,11,17,res,lambda);
939 add_triangle(11,17,18,res,lambda);
940 add_triangle(11,12,18,res,lambda);
941 add_triangle(12,18,19,res,lambda);
942 add_triangle(12,13,19,res,lambda);
943 add_triangle(13,19,20,res,lambda);
944 add_triangle(13,14,20,res,lambda);
945 add_triangle(14,20,21,res,lambda);
946
947 add_triangle(15,16,23,res,lambda);
948 add_triangle(16,17,24,res,lambda);
949 add_triangle(17,18,25,res,lambda);
950 add_triangle(18,19,26,res,lambda);
951 add_triangle(19,20,27,res,lambda);
952 add_triangle(20,21,28,res,lambda);
953 add_triangle(15,22,23,res,lambda);
954 add_triangle(16,23,24,res,lambda);
955 add_triangle(17,24,25,res,lambda);
956 add_triangle(18,25,26,res,lambda);
957 add_triangle(19,26,27,res,lambda);
958 add_triangle(20,27,28,res,lambda);
959 add_triangle(21,28,29,res,lambda);
960
961 add_triangle(22,23,30,res,lambda);
962 add_triangle(23,24,31,res,lambda);
963 add_triangle(24,25,32,res,lambda);
964 add_triangle(25,26,33,res,lambda);
965 add_triangle(26,27,34,res,lambda);
966 add_triangle(27,28,35,res,lambda);
967 add_triangle(28,29,36,res,lambda);
968 add_triangle(23,30,31,res,lambda);
969 add_triangle(24,31,32,res,lambda);
970 add_triangle(25,32,33,res,lambda);
971 add_triangle(26,33,34,res,lambda);
972 add_triangle(27,34,35,res,lambda);
973 add_triangle(28,35,36,res,lambda);
974
975 add_triangle(30,31,37,res,lambda);
976 add_triangle(31,32,38,res,lambda);
977 add_triangle(32,33,39,res,lambda);
978 add_triangle(33,34,40,res,lambda);
979 add_triangle(34,35,41,res,lambda);
980 add_triangle(35,36,42,res,lambda);
981 add_triangle(31,37,38,res,lambda);
982 add_triangle(32,38,39,res,lambda);
983 add_triangle(33,39,40,res,lambda);
984 add_triangle(34,40,41,res,lambda);
985 add_triangle(35,41,42,res,lambda);
986
987 add_triangle(37,38,43,res,lambda);
988 add_triangle(38,39,44,res,lambda);
989 add_triangle(39,40,45,res,lambda);
990 add_triangle(40,41,46,res,lambda);
991 add_triangle(41,42,47,res,lambda);
992 add_triangle(38,43,44,res,lambda);
993 add_triangle(39,44,45,res,lambda);
994 add_triangle(40,45,46,res,lambda);
995 add_triangle(41,46,47,res,lambda);
996 }
997
998 res += res.transpose().eval();
999
1000 auto res_ = compress_CuthillMcKee(res,true);
1001 res = res_;
1002
1003 return res;
1004}
1005
1006ArrayXXd hopping_Archimedean (string vertex_conf, int VARIANT=0, double lambda1=1., double lambda2=1.)
1007{
1008 ArrayXXd res;
1009
1010 if (vertex_conf == "3.5.3.5") // icosidodecahedron
1011 {
1012 int L=30;
1013 res.resize(L,L); res.setZero();
1014
1015 if (VARIANT==1 or VARIANT==0) // my naive counting
1016 {
1017 res(0,1) = lambda1;
1018 res(1,2) = lambda1;
1019 res(2,3) = lambda1;
1020 res(3,4) = lambda1;
1021 res(0,4) = lambda1;
1022
1023 res(0,5) = lambda1;
1024 res(0,6) = lambda1;
1025 res(1,6) = lambda1;
1026 res(1,7) = lambda1;
1027 res(2,7) = lambda1;
1028 res(2,8) = lambda1;
1029 res(3,8) = lambda1;
1030 res(3,9) = lambda1;
1031 res(4,5) = lambda1;
1032 res(4,9) = lambda1;
1033
1034 res(5,11) = lambda1;
1035 res(5,12) = lambda1;
1036 res(6,13) = lambda1;
1037 res(6,14) = lambda1;
1038 res(7,15) = lambda1;
1039 res(7,16) = lambda1;
1040 res(8,17) = lambda1;
1041 res(8,18) = lambda1;
1042 res(9,10) = lambda1;
1043 res(9,19) = lambda1;
1044
1045 res(12,13) = lambda1;
1046 res(13,14) = lambda1;
1047 res(14,15) = lambda1;
1048 res(15,16) = lambda1;
1049 res(16,17) = lambda1;
1050 res(17,18) = lambda1;
1051 res(18,19) = lambda1;
1052 res(10,19) = lambda1;
1053 res(10,11) = lambda1;
1054 res(11,12) = lambda1;
1055
1056 res(12,22) = lambda1;
1057 res(13,22) = lambda1;
1058 res(14,23) = lambda1;
1059 res(15,23) = lambda1;
1060 res(16,24) = lambda1;
1061 res(17,24) = lambda1;
1062 res(18,20) = lambda1;
1063 res(19,20) = lambda1;
1064 res(10,21) = lambda1;
1065 res(11,21) = lambda1;
1066
1067 res(20,28) = lambda1;
1068 res(20,29) = lambda1;
1069 res(21,27) = lambda1;
1070 res(21,28) = lambda1;
1071 res(22,26) = lambda1;
1072 res(22,27) = lambda1;
1073 res(23,25) = lambda1;
1074 res(23,26) = lambda1;
1075 res(24,25) = lambda1;
1076 res(24,29) = lambda1;
1077
1078 res(25,26) = lambda1;
1079 res(26,27) = lambda1;
1080 res(27,28) = lambda1;
1081 res(28,29) = lambda1;
1082 res(25,29) = lambda1;
1083 }
1084 else if (VARIANT==2) // According to Exler, Schnack (2003)
1085 {
1086 res(1-1,3-1) = lambda1;
1087 res(1-1,4-1) = lambda1;
1088 res(1-1,28-1) = lambda1;
1089 res(1-1,29-1) = lambda1;
1090
1091 res(2-1,5-1) = lambda1;
1092 res(2-1,8-1) = lambda1;
1093 res(2-1,26-1) = lambda1;
1094 res(2-1,29-1) = lambda1;
1095
1096 res(3-1,4-1) = lambda1;
1097 res(3-1,6-1) = lambda1;
1098 res(3-1,30-1) = lambda1;
1099
1100 res(4-1,5-1) = lambda1;
1101 res(4-1,7-1) = lambda1;
1102
1103 res(5-1,7-1) = lambda1;
1104 res(5-1,8-1) = lambda1;
1105
1106 res(6-1,30-1) = lambda1;
1107 res(6-1,9-1) = lambda1;
1108 res(6-1,12-1) = lambda1;
1109
1110 res(7-1,9-1) = lambda1;
1111 res(7-1,10-1) = lambda1;
1112
1113 res(8-1,11-1) = lambda1;
1114 res(8-1,14-1) = lambda1;
1115
1116 res(9-1,10-1) = lambda1;
1117 res(9-1,12-1) = lambda1;
1118
1119 res(10-1,11-1) = lambda1;
1120 res(10-1,13-1) = lambda1;
1121
1122 res(11-1,13-1) = lambda1;
1123 res(11-1,14-1) = lambda1;
1124
1125 res(12-1,15-1) = lambda1;
1126 res(12-1,18-1) = lambda1;
1127
1128 res(13-1,15-1) = lambda1;
1129 res(13-1,16-1) = lambda1;
1130
1131 res(14-1,17-1) = lambda1;
1132 res(14-1,20-1) = lambda1;
1133
1134 res(15-1,16-1) = lambda1;
1135 res(15-1,18-1) = lambda1;
1136
1137 res(16-1,17-1) = lambda1;
1138 res(16-1,19-1) = lambda1;
1139
1140 res(17-1,19-1) = lambda1;
1141 res(17-1,20-1) = lambda1;
1142
1143 res(18-1,21-1) = lambda1;
1144 res(18-1,24-1) = lambda1;
1145
1146 res(19-1,21-1) = lambda1;
1147 res(19-1,22-1) = lambda1;
1148
1149 res(20-1,23-1) = lambda1;
1150 res(20-1,26-1) = lambda1;
1151
1152 res(21-1,22-1) = lambda1;
1153 res(21-1,24-1) = lambda1;
1154
1155 res(22-1,23-1) = lambda1;
1156 res(22-1,25-1) = lambda1;
1157
1158 res(23-1,25-1) = lambda1;
1159 res(23-1,26-1) = lambda1;
1160
1161 res(24-1,27-1) = lambda1;
1162 res(24-1,30-1) = lambda1;
1163
1164 res(25-1,27-1) = lambda1;
1165 res(25-1,28-1) = lambda1;
1166
1167 res(26-1,29-1) = lambda1;
1168
1169 res(27-1,28-1) = lambda1;
1170 res(27-1,30-1) = lambda1;
1171
1172 res(28-1,29-1) = lambda1;
1173 }
1174 }
1175 else if (vertex_conf == "3^4.5") // snub dodecahedron
1176 {
1177 int L=60;
1178 res.resize(L,L); res.setZero();
1179
1180 for (int i=0; i<=58; ++i) res(i,i+1) = lambda1;
1181
1182 res( 0, 4) = lambda1;
1183 res( 0, 6) = lambda1;
1184 res( 0, 7) = lambda1;
1185 res( 0, 8) = lambda1;
1186
1187 res( 1, 8) = lambda1;
1188 res( 1, 9) = lambda1;
1189 res( 1,10) = lambda1;
1190
1191 res( 2,10) = lambda1;
1192 res( 2,11) = lambda1;
1193 res( 2,12) = lambda1;
1194
1195 res( 3,12) = lambda1;
1196 res( 3,13) = lambda1;
1197 res( 3,14) = lambda1;
1198
1199 res( 4, 6) = lambda1;
1200 res( 4,14) = lambda1;
1201
1202 res( 5,14) = lambda1;
1203 res( 5,15) = lambda1;
1204 res( 5,16) = lambda1;
1205
1206 res( 6,18) = lambda1;
1207
1208 res( 7,18) = lambda1;
1209 res( 7,19) = lambda1;
1210
1211 res( 8,21) = lambda1;
1212
1213 res( 9,21) = lambda1;
1214 res( 9,22) = lambda1;
1215
1216 res(10,24) = lambda1;
1217
1218 res(11,24) = lambda1;
1219 res(11,25) = lambda1;
1220
1221 res(12,27) = lambda1;
1222
1223 res(13,27) = lambda1;
1224 res(13,28) = lambda1;
1225
1226 res(15,29) = lambda1;
1227 res(15,30) = lambda1;
1228
1229 res(16,30) = lambda1;
1230 res(16,31) = lambda1;
1231
1232 res(17,31) = lambda1;
1233 res(17,32) = lambda1;
1234 res(17,33) = lambda1;
1235
1236 res(18,33) = lambda1;
1237
1238 res(19,33) = lambda1;
1239 res(19,34) = lambda1;
1240
1241 res(20,34) = lambda1;
1242 res(20,35) = lambda1;
1243 res(20,36) = lambda1;
1244
1245 res(21,36) = lambda1;
1246
1247 res(22,36) = lambda1;
1248 res(22,37) = lambda1;
1249
1250 res(23,37) = lambda1;
1251 res(23,38) = lambda1;
1252 res(23,39) = lambda1;
1253
1254 res(24,39) = lambda1;
1255
1256 res(25,39) = lambda1;
1257 res(25,40) = lambda1;
1258
1259 res(26,40) = lambda1;
1260 res(26,41) = lambda1;
1261 res(26,42) = lambda1;
1262
1263 res(27,42) = lambda1;
1264
1265 res(28,42) = lambda1;
1266 res(28,43) = lambda1;
1267
1268 res(29,43) = lambda1;
1269 res(29,44) = lambda1;
1270
1271 res(30,44) = lambda1;
1272
1273 res(31,46) = lambda1;
1274
1275 res(32,46) = lambda1;
1276 res(32,47) = lambda1;
1277
1278 res(34,48) = lambda1;
1279
1280 res(35,48) = lambda1;
1281 res(35,49) = lambda1;
1282
1283 res(37,50) = lambda1;
1284
1285 res(38,50) = lambda1;
1286 res(38,51) = lambda1;
1287
1288 res(40,52) = lambda1;
1289
1290 res(41,52) = lambda1;
1291 res(41,53) = lambda1;
1292
1293 res(43,54) = lambda1;
1294
1295 res(44,54) = lambda1;
1296
1297 res(45,54) = lambda1;
1298 res(45,55) = lambda1;
1299 res(45,56) = lambda1;
1300
1301 res(46,56) = lambda1;
1302
1303 res(47,56) = lambda1;
1304 res(47,57) = lambda1;
1305
1306 res(48,57) = lambda1;
1307
1308 res(49,57) = lambda1;
1309 res(49,58) = lambda1;
1310
1311 res(50,58) = lambda1;
1312
1313 res(51,58) = lambda1;
1314 res(51,59) = lambda1;
1315
1316 res(52,59) = lambda1;
1317
1318 res(53,55) = lambda1;
1319 res(53,59) = lambda1;
1320
1321 res(55,59) = lambda1;
1322 }
1323 else if (vertex_conf == "4.6^2") // truncated octahedron
1324 {
1325 int L=24;
1326 res.resize(L,L); res.setZero();
1327
1328 for (int i=0; i<=4; ++i) res(i,i+1) = lambda1;
1329 res(0,5) = lambda1;
1330
1331 res(5,6) = lambda1;
1332 res(0,9) = lambda1;
1333 res(1,10) = lambda1;
1334 res(2,13) = lambda1;
1335 res(3,14) = lambda1;
1336 res(4,17) = lambda1;
1337
1338 for (int i=6; i<=16; ++i) res(i,i+1) = lambda1;
1339 res(6,17) = lambda1;
1340
1341 res(7,19) = lambda1;
1342 res(8,20) = lambda1;
1343 res(11,21) = lambda1;
1344 res(12,22) = lambda1;
1345 res(15,23) = lambda1;
1346 res(16,18) = lambda1;
1347
1348 for (int i=18; i<=22; ++i) res(i,i+1) = lambda1;
1349 res(18,23) = lambda1;
1350 }
1351 else if (vertex_conf == "3.8^2") // truncated cube
1352 {
1353 int L=24;
1354 res.resize(L,L); res.setZero();
1355
1356 for (int i=0; i<=6; ++i) res(i,i+1) = lambda1;
1357 res(0,7) = lambda1;
1358
1359 res(0,8) = lambda1;
1360 res(7,8) = lambda1;
1361 res(1,9) = lambda1;
1362 res(2,9) = lambda1;
1363 res(3,10) = lambda1;
1364 res(4,10) = lambda1;
1365 res(5,11) = lambda1;
1366 res(6,11) = lambda1;
1367
1368 res(8,13) = lambda1;
1369 res(9,14) = lambda1;
1370 res(10,15) = lambda1;
1371 res(11,12) = lambda1;
1372
1373 res(12,17) = lambda1;
1374 res(12,18) = lambda1;
1375 res(13,19) = lambda1;
1376 res(13,20) = lambda1;
1377 res(14,21) = lambda1;
1378 res(14,22) = lambda1;
1379 res(15,16) = lambda1;
1380 res(15,23) = lambda1;
1381
1382 for (int i=15; i<=22; ++i) res(i,i+1) = lambda1;
1383 res(15,23) = lambda1;
1384
1385 res(19,20) = lambda1;
1386 res(16,23) = lambda1;
1387 }
1388 else if (vertex_conf == "3.4.3.4") // cuboctahedron
1389 {
1390 int L=12;
1391 res.resize(L,L); res.setZero();
1392
1393 for (int i=0; i<=2; ++i) res(i,i+1) = lambda1;
1394 res(0,3) = lambda1;
1395
1396 res(0,4) = lambda1;
1397 res(3,4) = lambda1;
1398 res(0,5) = lambda1;
1399 res(1,5) = lambda1;
1400 res(1,6) = lambda1;
1401 res(2,6) = lambda1;
1402 res(2,7) = lambda1;
1403 res(3,7) = lambda1;
1404
1405 res(4,8) = lambda1;
1406 res(7,8) = lambda1;
1407 res(4,9) = lambda1;
1408 res(5,9) = lambda1;
1409 res(5,10) = lambda1;
1410 res(6,10) = lambda1;
1411 res(6,11) = lambda1;
1412 res(7,11) = lambda1;
1413
1414 for (int i=8; i<=10; ++i) res(i,i+1) = lambda1;
1415 res(8,11) = lambda1;
1416 }
1417 else if (vertex_conf == "3.6^2") // truncated tetrahedron = C12
1418 {
1419 int L=12;
1420 res.resize(L,L); res.setZero();
1421
1422 // A
1423 res(0,1) = lambda1;
1424 res(0,3) = lambda1;
1425 res(1,3) = lambda1;
1426
1427 // B
1428 res(2,6) = lambda1;
1429 res(2,5) = lambda1;
1430 res(5,6) = lambda1;
1431
1432 // C
1433 res(4,9) = lambda1;
1434 res(4,8) = lambda1;
1435 res(8,9) = lambda1;
1436
1437 // D
1438 res(7,10) = lambda1;
1439 res(7,11) = lambda1;
1440 res(10,11) = lambda1;
1441
1442 res(3,7) = lambda2; // A-D
1443 res(6,10) = lambda2; // B-D
1444 res(9,11) = lambda2; // C-D <--
1445
1446 res(0,2) = lambda2; // A-B <--
1447 res(1,4) = lambda2; // A-C
1448 res(5,8) = lambda2; // B-C
1449 }
1450
1451 res += res.transpose().eval();
1452
1453 if (VARIANT==0 and vertex_conf != "3.6^2" and vertex_conf != "3.4.3.4")
1454 {
1455 auto res_ = compress_CuthillMcKee(res,true);
1456 res = res_;
1457 }
1458
1459 return res;
1460}
1461
1462ArrayXXd hopping_fullerene_C40Td (int VARIANT=0, double lambda1=1., double lambda2=1.)
1463{
1464 int L = 40;
1465 ArrayXXd res(L,L); res.setZero();
1466
1467 res( 0 , 2 ) = lambda1;
1468 res( 0 , 26 ) = lambda1;
1469 res( 0 , 1 ) = lambda1;
1470 res( 1 , 9 ) = lambda1;
1471 res( 1 , 28 ) = lambda1;
1472 res( 2 , 3 ) = lambda1;
1473 res( 2 , 5 ) = lambda1;
1474 res( 3 , 4 ) = lambda1;
1475 res( 3 , 15 ) = lambda1;
1476 res( 4 , 7 ) = lambda1;
1477 res( 4 , 16 ) = lambda1;
1478 res( 5 , 6 ) = lambda1;
1479 res( 5 , 9 ) = lambda1;
1480 res( 6 , 7 ) = lambda1;
1481 res( 6 , 10 ) = lambda1;
1482 res( 7 , 18 ) = lambda1;
1483 res( 8 , 9 ) = lambda1;
1484 res( 8 , 12 ) = lambda1;
1485 res( 8 , 11 ) = lambda1;
1486 res( 10 , 11 ) = lambda1;
1487 res( 10 , 37 ) = lambda1;
1488 res( 11 , 38 ) = lambda1;
1489 res( 12 , 39 ) = lambda1;
1490 res( 12 , 29 ) = lambda1;
1491 res( 13 , 15 ) = lambda1;
1492 res( 13 , 19 ) = lambda1;
1493 res( 13 , 14 ) = lambda1;
1494 res( 14 , 16 ) = lambda1;
1495 res( 14 , 23 ) = lambda1;
1496 res( 15 , 26 ) = lambda1;
1497 res( 16 , 17 ) = lambda1;
1498 res( 17 , 18 ) = lambda1;
1499 res( 17 , 21 ) = lambda1;
1500 res( 18 , 37 ) = lambda1;
1501 res( 19 , 20 ) = lambda1;
1502 res( 19 , 24 ) = lambda1;
1503 res( 20 , 23 ) = lambda1;
1504 res( 20 , 31 ) = lambda1;
1505 res( 21 , 22 ) = lambda1;
1506 res( 21 , 35 ) = lambda1;
1507 res( 22 , 23 ) = lambda1;
1508 res( 22 , 32 ) = lambda1;
1509 res( 24 , 25 ) = lambda1;
1510 res( 24 , 27 ) = lambda1;
1511 res( 25 , 26 ) = lambda1;
1512 res( 25 , 28 ) = lambda1;
1513 res( 27 , 30 ) = lambda1;
1514 res( 27 , 31 ) = lambda1;
1515 res( 28 , 29 ) = lambda1;
1516 res( 29 , 30 ) = lambda1;
1517 res( 30 , 34 ) = lambda1;
1518 res( 31 , 32 ) = lambda1;
1519 res( 32 , 33 ) = lambda1;
1520 res( 33 , 34 ) = lambda1;
1521 res( 33 , 36 ) = lambda1;
1522 res( 34 , 39 ) = lambda1;
1523 res( 35 , 37 ) = lambda1;
1524 res( 35 , 36 ) = lambda1;
1525 res( 36 , 38 ) = lambda1;
1526 res( 38 , 39 ) = lambda1;
1527
1528 res += res.transpose().eval();
1529
1530 if (VARIANT==0)
1531 {
1532 auto res_ = compress_CuthillMcKee(res,true);
1533 res = res_;
1534 }
1535
1536 return res;
1537}
1538
1539// reference: PRB 93, 165406 (2016), Appendix C
1540ArrayXXd hopping_fullerene (int L=60, int VARIANT=0, double lambda1=1., double lambda2=1.)
1541{
1542 ArrayXXd res(L,L); res.setZero();
1543
1544 // lambda1: pentagon bond, lambda2: hexagon bond
1545 if (L == 60)
1546 {
1547 if (VARIANT==1) // inwards spiral
1548 {
1549 for (int i=0; i<=3; ++i) res(i,i+1) = lambda1;
1550 res(0,4) = lambda1;
1551
1552 res(0,8) = lambda2; //1
1553 res(1,11) = lambda2; //2
1554 res(2,14) = lambda2; //3
1555 res(3,17) = lambda2; //4
1556 res(4,5) = lambda2; //5
1557
1558 for (int i=5; i<=18; ++i) res(i,i+1) = lambda1;
1559 res(6,7) = lambda2; //6
1560 res(9,10) = lambda2; //7
1561 res(12,13) = lambda2; //8
1562 res(15,16) = lambda2; //9
1563 res(18,19) = lambda2; //10
1564 res(5,19) = lambda1;
1565
1566 res(7,24) = lambda1;
1567 res(9,25) = lambda1;
1568 res(10,28) = lambda1;
1569 res(12,29) = lambda1;
1570 res(13,32) = lambda1;
1571 res(15,33) = lambda1;
1572 res(16,36) = lambda1;
1573 res(18,37) = lambda1;
1574 res(19,20) = lambda1;
1575 res(6,21) = lambda1;
1576
1577 for (int i=20; i<=38; ++i) res(i,i+1) = lambda1;
1578 res(21,22) = lambda2; //11
1579 res(23,24) = lambda2; //12
1580 res(25,26) = lambda2; //13
1581 res(27,28) = lambda2; //14
1582 res(29,30) = lambda2; //15
1583 res(31,32) = lambda2; //16
1584 res(33,34) = lambda2; //17
1585 res(35,36) = lambda2; //18
1586 res(37,38) = lambda2; //19
1587 res(20,39) = lambda2; //20
1588
1589 res(26,44) = lambda1;
1590 res(27,46) = lambda1;
1591 res(30,47) = lambda1;
1592 res(31,49) = lambda1;
1593 res(34,50) = lambda1;
1594 res(35,52) = lambda1;
1595 res(38,53) = lambda1;
1596 res(39,40) = lambda1;
1597 res(22,41) = lambda1;
1598 res(23,43) = lambda1;
1599
1600 for (int i=40; i<=54; ++i) res(i,i+1) = lambda1;
1601 res(40,41) = lambda2; //21
1602 res(43,44) = lambda2; //22
1603 res(46,47) = lambda2; //23
1604 res(49,50) = lambda2; //24
1605 res(52,53) = lambda2; //25
1606 res(40,54) = lambda1;
1607
1608 res(45,57) = lambda2; //26
1609 res(48,58) = lambda2; //27
1610 res(51,59) = lambda2; //28
1611 res(54,55) = lambda2; //29
1612 res(42,56) = lambda2; //30
1613
1614 for (int i=55; i<=58; ++i) res(i,i+1) = lambda1;
1615 res(55,59) = lambda1;
1616 }
1617 // reference: https://www.qmul.ac.uk/sbcs/iupac/fullerene2/311.html
1618 // also in: Phys. Rev. B 72, 064453 (2005)
1619 else if (VARIANT==2 or VARIANT==0) // Kivelson
1620 {
1621 res(0,4) = lambda1;
1622 res(3,4) = lambda1;
1623 res(2,3) = lambda1;
1624 res(1,2) = lambda1;
1625 res(0,1) = lambda1;
1626
1627 res(0,8) = lambda1;
1628 res(4,5) = lambda1;
1629 res(3,17) = lambda1;
1630 res(2,14) = lambda1;
1631 res(1,11) = lambda1;
1632
1633 res(7,8) = lambda1;
1634 res(6,7) = lambda1;
1635 res(5,6) = lambda1;
1636 res(5,19) = lambda1;
1637 res(18,19) = lambda1;
1638 res(17,18) = lambda1;
1639 res(16,17) = lambda1;
1640 res(15,16) = lambda1;
1641 res(14,15) = lambda1;
1642 res(13,14) = lambda1;
1643 res(12,13) = lambda1;
1644 res(11,12) = lambda1;
1645 res(10,11) = lambda1;
1646 res(9,10) = lambda1;
1647 res(8,9) = lambda1;
1648
1649 res(7,22) = lambda1;
1650 res(6,21) = lambda1;
1651 res(19,20) = lambda1;
1652 res(18,29) = lambda1;
1653 res(16,28) = lambda1;
1654 res(15,27) = lambda1;
1655 res(13,26) = lambda1;
1656 res(12,25) = lambda1;
1657 res(10,24) = lambda1;
1658 res(9,23) = lambda1;
1659
1660 res(22,33) = lambda1;
1661 res(32,33) = lambda1;
1662 res(21,32) = lambda1;
1663 res(20,21) = lambda1;
1664 res(20,31) = lambda1;
1665 res(30,31) = lambda1;
1666 res(29,30) = lambda1;
1667 res(28,29) = lambda1;
1668 res(28,39) = lambda1;
1669 res(38,39) = lambda1;
1670 res(27,38) = lambda1;
1671 res(26,27) = lambda1;
1672 res(26,37) = lambda1;
1673 res(36,37) = lambda1;
1674 res(25,36) = lambda1;
1675 res(24,25) = lambda1;
1676 res(24,35) = lambda1;
1677 res(34,35) = lambda1;
1678 res(23,34) = lambda1;
1679 res(22,23) = lambda1;
1680
1681 res(33,46) = lambda1;
1682 res(32,44) = lambda1;
1683 res(31,43) = lambda1;
1684 res(30,41) = lambda1;
1685 res(39,40) = lambda1;
1686 res(38,53) = lambda1;
1687 res(37,52) = lambda1;
1688 res(36,50) = lambda1;
1689 res(35,49) = lambda1;
1690 res(34,47) = lambda1;
1691
1692 res(45,46) = lambda1;
1693 res(44,45) = lambda1;
1694 res(43,44) = lambda1;
1695 res(42,43) = lambda1;
1696 res(41,42) = lambda1;
1697 res(40,41) = lambda1;
1698 res(40,54) = lambda1;
1699 res(53,54) = lambda1;
1700 res(52,53) = lambda1;
1701 res(51,52) = lambda1;
1702 res(50,51) = lambda1;
1703 res(49,50) = lambda1;
1704 res(48,49) = lambda1;
1705 res(47,48) = lambda1;
1706 res(46,47) = lambda1;
1707
1708 res(45,57) = lambda1;
1709 res(42,56) = lambda1;
1710 res(54,55) = lambda1;
1711 res(51,59) = lambda1;
1712 res(48,58) = lambda1;
1713 res(56,57) = lambda1;
1714 res(55,56) = lambda1;
1715 res(55,59) = lambda1;
1716 res(58,59) = lambda1;
1717 res(57,58) = lambda1;
1718 }
1719 }
1720 else if (L == 20)
1721 {
1722// res(11,12) = lambda1;
1723// res(12,13) = lambda1;
1724// res(13,14) = lambda1;
1725// res(5,14) = lambda1;
1726// res(5,6) = lambda1;
1727// res(6,7) = lambda1;
1728// res(7,8) = lambda1;
1729// res(8,9) = lambda1;
1730// res(9,10) = lambda1;
1731// res(10,11) = lambda1;
1732//
1733// res(0,1) = lambda1;
1734// res(1,2) = lambda1;
1735// res(2,3) = lambda1;
1736// res(3,4) = lambda1;
1737// res(0,4) = lambda1;
1738//
1739// res(2,11) = lambda1;
1740// res(3,13) = lambda1;
1741// res(4,5) = lambda1;
1742// res(0,7) = lambda1;
1743// res(1,9) = lambda1;
1744// res(2,11) = lambda1;
1745//
1746// res(15,16) = lambda1;
1747// res(16,17) = lambda1;
1748// res(17,18) = lambda1;
1749// res(18,19) = lambda1;
1750// res(15,19) = lambda1;
1751//
1752// res(12,19) = lambda1;
1753// res(14,15) = lambda1;
1754// res(6,16) = lambda1;
1755// res(8,17) = lambda1;
1756// res(10,18) = lambda1;
1757
1758 // better numbering (inwards spiral):
1759
1760 res(0,1) = lambda1;
1761 res(1,2) = lambda1;
1762 res(2,3) = lambda1;
1763 res(3,4) = lambda1;
1764 res(0,4) = lambda1;
1765
1766 res(0,7) = lambda1;
1767 res(1,9) = lambda1;
1768 res(2,11) = lambda1;
1769 res(3,13) = lambda1;
1770 res(4,5) = lambda1;
1771
1772 res(5,6) = lambda1;
1773 res(6,7) = lambda1;
1774 res(7,8) = lambda1;
1775 res(8,9) = lambda1;
1776 res(9,10) = lambda1;
1777 res(10,11) = lambda1;
1778 res(11,12) = lambda1;
1779 res(12,13) = lambda1;
1780 res(13,14) = lambda1;
1781 res(5,14) = lambda1;
1782
1783 res(6,16) = lambda1;
1784 res(8,17) = lambda1;
1785 res(10,18) = lambda1;
1786 res(12,19) = lambda1;
1787 res(14,15) = lambda1;
1788
1789 res(15,16) = lambda1;
1790 res(16,17) = lambda1;
1791 res(17,18) = lambda1;
1792 res(18,19) = lambda1;
1793 res(15,19) = lambda1;
1794 }
1795 else if (L==40) // symmetry D_5d(I)
1796 {
1797 for (int i=0; i<=3; ++i) res(i,i+1) = lambda1;
1798 res(0,4) = lambda1;
1799
1800 res(0,8) = lambda1;
1801 res(1,11) = lambda1;
1802 res(2,14) = lambda1;
1803 res(3,17) = lambda1;
1804 res(4,5) = lambda1;
1805
1806 for (int i=5; i<=18; ++i) res(i,i+1) = lambda1;
1807 res(5,19) = lambda1;
1808
1809 res(6,21) = lambda1;
1810 res(7,23) = lambda1;
1811 res(9,24) = lambda1;
1812 res(10,26) = lambda1;
1813 res(12,27) = lambda1;
1814 res(13,29) = lambda1;
1815 res(15,30) = lambda1;
1816 res(16,32) = lambda1;
1817 res(18,33) = lambda1;
1818 res(19,20) = lambda1;
1819
1820 for (int i=20; i<=33; ++i) res(i,i+1) = lambda1;
1821 res(20,34) = lambda1;
1822
1823 res(22,36) = lambda1;
1824 res(25,37) = lambda1;
1825 res(28,38) = lambda1;
1826 res(31,39) = lambda1;
1827 res(34,35) = lambda1;
1828
1829 for (int i=35; i<=38; ++i) res(i,i+1) = lambda1;
1830 res(35,39) = lambda1;
1831 }
1832 else if (L==36) // symmetry D_6h // https://nanotube.msu.edu/fullerene/fullerene.php?C=36
1833 {
1834 for (int i=0; i<=4; ++i) res(i,i+1) = lambda1;
1835 res(0,5) = lambda1;
1836
1837 res(0,8) = lambda1;
1838 res(1,10) = lambda1;
1839 res(2,12) = lambda1;
1840 res(3,14) = lambda1;
1841 res(4,16) = lambda1;
1842 res(5,6) = lambda1;
1843
1844 for (int i=6; i<=16; ++i) res(i,i+1) = lambda1;
1845 res(6,17) = lambda1;
1846
1847 res(7,20) = lambda1;
1848 res(9,22) = lambda1;
1849 res(11,24) = lambda1;
1850 res(13,26) = lambda1;
1851 res(15,28) = lambda1;
1852 res(17,18) = lambda1;
1853
1854 for (int i=18; i<=28; ++i) res(i,i+1) = lambda1;
1855 res(18,29) = lambda1;
1856
1857 res(19,31) = lambda1;
1858 res(21,32) = lambda1;
1859 res(23,33) = lambda1;
1860 res(25,34) = lambda1;
1861 res(27,35) = lambda1;
1862 res(29,30) = lambda1;
1863
1864 for (int i=30; i<=34; ++i) res(i,i+1) = lambda1;
1865 res(30,35) = lambda1;
1866 }
1867 else if (L==30) // symmetry D5h // https://nanotube.msu.edu/fullerene/fullerene.php?C=30
1868 {
1869 for (int i=0; i<=3; ++i) res(i,i+1) = lambda1;
1870 res(0,4) = lambda1;
1871
1872 res(0,7) = lambda1;
1873 res(1,9) = lambda1;
1874 res(2,11) = lambda1;
1875 res(3,13) = lambda1;
1876 res(4,5) = lambda1;
1877
1878 for (int i=5; i<=13; ++i) res(i,i+1) = lambda1;
1879 res(5,14) = lambda1;
1880
1881 res(6,17) = lambda1;
1882 res(8,19) = lambda1;
1883 res(10,21) = lambda1;
1884 res(12,23) = lambda1;
1885 res(14,15) = lambda1;
1886
1887 for (int i=15; i<=23; ++i) res(i,i+1) = lambda1;
1888 res(15,24) = lambda1;
1889
1890 res(16,26) = lambda1;
1891 res(18,27) = lambda1;
1892 res(20,28) = lambda1;
1893 res(22,29) = lambda1;
1894 res(24,25) = lambda1;
1895
1896 for (int i=25; i<=28; ++i) res(i,i+1) = lambda1;
1897 res(25,29) = lambda1;
1898 }
1899 else if (L==28) // symmetry Td // https://nanotube.msu.edu/fullerene/fullerene.php?C=28
1900 {
1901 res( 0 , 4 ) = lambda1;
1902 res( 0 , 20 ) = lambda1;
1903 res( 0 , 3 ) = lambda1;
1904 res( 1 , 19 ) = lambda1;
1905 res( 1 , 17 ) = lambda1;
1906 res( 1 , 2 ) = lambda1;
1907 res( 2 , 7 ) = lambda1;
1908 res( 2 , 3 ) = lambda1;
1909 res( 3 , 5 ) = lambda1;
1910 res( 4 , 8 ) = lambda1;
1911 res( 4 , 6 ) = lambda1;
1912 res( 5 , 23 ) = lambda1;
1913 res( 5 , 6 ) = lambda1;
1914 res( 6 , 14 ) = lambda1;
1915 res( 7 , 23 ) = lambda1;
1916 res( 7 , 27 ) = lambda1;
1917 res( 8 , 11 ) = lambda1;
1918 res( 8 , 9 ) = lambda1;
1919 res( 9 , 20 ) = lambda1;
1920 res( 9 , 10 ) = lambda1;
1921 res( 10 , 16 ) = lambda1;
1922 res( 10 , 13 ) = lambda1;
1923 res( 11 , 14 ) = lambda1;
1924 res( 11 , 12 ) = lambda1;
1925 res( 12 , 22 ) = lambda1;
1926 res( 12 , 13 ) = lambda1;
1927 res( 13 , 18 ) = lambda1;
1928 res( 14 , 25 ) = lambda1;
1929 res( 15 , 18 ) = lambda1;
1930 res( 15 , 16 ) = lambda1;
1931 res( 15 , 17 ) = lambda1;
1932 res( 16 , 19 ) = lambda1;
1933 res( 17 , 27 ) = lambda1;
1934 res( 18 , 21 ) = lambda1;
1935 res( 19 , 20 ) = lambda1;
1936 res( 21 , 22 ) = lambda1;
1937 res( 21 , 26 ) = lambda1;
1938 res( 22 , 25 ) = lambda1;
1939 res( 23 , 24 ) = lambda1;
1940 res( 24 , 25 ) = lambda1;
1941 res( 24 , 26 ) = lambda1;
1942 res( 26 , 27 ) = lambda1;
1943 }
1944 else if (L==26)
1945 {
1946 for (int i=0; i<=3; ++i) res(i,i+1) = lambda1;
1947 res(0,4) = lambda1;
1948
1949 res(0,7) = lambda1;
1950 res(1,10) = lambda1;
1951 res(2,12) = lambda1;
1952 res(3,14) = lambda1;
1953 res(4,5) = lambda1;
1954
1955 for (int i=5; i<=14; ++i) res(i,i+1) = lambda1;
1956 res(5,15) = lambda1;
1957
1958 res(6,18) = lambda1;
1959 res(8,19) = lambda1;
1960 res(9,21) = lambda1;
1961 res(11,22) = lambda1;
1962 res(13,25) = lambda1;
1963 res(15,16) = lambda1;
1964
1965 for (int i=16; i<=23; ++i) res(i,i+1) = lambda1;
1966
1967 res(16,25) = lambda1;
1968 res(17,24) = lambda1;
1969 res(20,24) = lambda1;
1970 res(23,25) = lambda1;
1971 }
1972 else if (L==24)
1973 {
1974 for (int i=0; i<=4; ++i) res(i,i+1) = lambda1;
1975 res(0,5) = lambda1;
1976
1977 res(0,8) = lambda1;
1978 res(1,10) = lambda1;
1979 res(2,12) = lambda1;
1980 res(3,14) = lambda1;
1981 res(4,16) = lambda1;
1982 res(5,6) = lambda1;
1983
1984 for (int i=6; i<=16; ++i) res(i,i+1) = lambda1;
1985 res(6,17) = lambda1;
1986
1987 res(17,18) = lambda1;
1988 res(7,19) = lambda1;
1989 res(9,20) = lambda1;
1990 res(11,21) = lambda1;
1991 res(13,22) = lambda1;
1992 res(15,23) = lambda1;
1993
1994 for (int i=18; i<=22; ++i) res(i,i+1) = lambda1;
1995 res(18,23) = lambda1;
1996 }
1997 else if (L==12)
1998 {
1999 return hopping_Archimedean("3.6^2",VARIANT,lambda1,lambda2);
2000 }
2001
2002 res += res.transpose().eval();
2003
2004 if (VARIANT==0)
2005 {
2006 auto res_ = compress_CuthillMcKee(res,true);
2007 res = res_;
2008 }
2009
2010 return res;
2011}
2012
2013ArrayXXd hopping_Platonic (int L, int VARIANT=0, double lambda1=1.)
2014{
2015 ArrayXXd res(L,L); res.setZero();
2016
2017 if (L == 4) // tetrahedron
2018 {
2019 res(0,1) = lambda1;
2020 res(1,2) = lambda1;
2021 res(0,2) = lambda1;
2022 res(0,3) = lambda1;
2023 res(1,3) = lambda1;
2024 res(2,3) = lambda1;
2025 }
2026 if (L == 5) // bipyramid
2027 {
2028 res(0,1) = lambda1;
2029 res(0,2) = lambda1;
2030 res(0,3) = lambda1;
2031 res(1,2) = lambda1;
2032 res(1,3) = lambda1;
2033 res(2,3) = lambda1;
2034 res(1,4) = lambda1;
2035 res(2,4) = lambda1;
2036 res(3,4) = lambda1;
2037 }
2038 else if (L == 6) // octahedron
2039 {
2040 res(0,1) = lambda1;
2041 res(1,2) = lambda1;
2042 res(2,3) = lambda1;
2043 res(0,3) = lambda1;
2044
2045 res(0,4) = lambda1;
2046 res(1,4) = lambda1;
2047 res(2,4) = lambda1;
2048 res(3,4) = lambda1;
2049
2050 res(0,5) = lambda1;
2051 res(1,5) = lambda1;
2052 res(2,5) = lambda1;
2053 res(3,5) = lambda1;
2054 }
2055 else if (L == 8) // cube
2056 {
2057 res(0,1) = lambda1;
2058 res(1,2) = lambda1;
2059 res(2,3) = lambda1;
2060 res(0,3) = lambda1;
2061
2062 res(4,5) = lambda1;
2063 res(5,6) = lambda1;
2064 res(6,7) = lambda1;
2065 res(4,7) = lambda1;
2066
2067 res(0,4) = lambda1;
2068 res(1,5) = lambda1;
2069 res(2,6) = lambda1;
2070 res(3,7) = lambda1;
2071 }
2072 // reference: Phys. Rev. B 72, 064453 (2005)
2073 else if (L == 12) // icosahedron
2074 {
2075// res(0,1) = lambda1;
2076// res(1,2) = lambda1;
2077// res(0,2) = lambda1;
2078//
2079// res(1,5) = lambda1;
2080// res(2,5) = lambda1;
2081// res(2,6) = lambda1;
2082// res(2,3) = lambda1;
2083// res(0,3) = lambda1;
2084// res(0,7) = lambda1;
2085// res(0,4) = lambda1;
2086// res(1,4) = lambda1;
2087// res(1,8) = lambda1;
2088//
2089// res(4,8) = lambda1;
2090// res(5,8) = lambda1;
2091// res(5,6) = lambda1;
2092// res(3,6) = lambda1;
2093// res(3,7) = lambda1;
2094// res(4,7) = lambda1;
2095//
2096// res(8,10) = lambda1;
2097// res(5,10) = lambda1;
2098// res(6,10) = lambda1;
2099// res(6,11) = lambda1;
2100// res(3,11) = lambda1;
2101// res(7,11) = lambda1;
2102// res(7,9) = lambda1;
2103// res(4,9) = lambda1;
2104// res(9,8) = lambda1;
2105//
2106// res(9,10) = lambda1;
2107// res(10,11) = lambda1;
2108// res(9,11) = lambda1;
2109
2110 // better numbering (inwards spiral):
2111
2112 res(0,1) = lambda1;
2113 res(1,2) = lambda1;
2114 res(0,2) = lambda1;
2115
2116 res(0,3) = lambda1;
2117 res(0,4) = lambda1;
2118 res(0,5) = lambda1;
2119
2120 res(1,5) = lambda1;
2121 res(1,6) = lambda1;
2122 res(1,7) = lambda1;
2123
2124 res(2,3) = lambda1;
2125 res(2,7) = lambda1;
2126 res(2,8) = lambda1;
2127
2128 res(3,4) = lambda1;
2129 res(4,5) = lambda1;
2130 res(5,6) = lambda1;
2131 res(6,7) = lambda1;
2132 res(7,8) = lambda1;
2133 res(3,8) = lambda1;
2134
2135 res(3,9) = lambda1;
2136 res(4,9) = lambda1;
2137 res(4,10) = lambda1;
2138 res(5,10) = lambda1;
2139 res(6,10) = lambda1;
2140 res(6,11) = lambda1;
2141 res(7,11) = lambda1;
2142 res(8,11) = lambda1;
2143 res(8,9) = lambda1;
2144
2145 res(9,10) = lambda1;
2146 res(10,11) = lambda1;
2147 res(9,11) = lambda1;
2148 }
2149 else if (L==20) // dodecahedron
2150 {
2151 res = hopping_fullerene(L, VARIANT, lambda1, lambda1);
2152 }
2153 else if (L==8) // cube
2154 {
2155 res(0,1) = lambda1;
2156 res(1,2) = lambda1;
2157 res(2,3) = lambda1;
2158 res(0,3) = lambda1;
2159
2160 res(0,5) = lambda1;
2161 res(1,6) = lambda1;
2162 res(2,7) = lambda1;
2163 res(3,4) = lambda1;
2164
2165 res(4,5) = lambda1;
2166 res(5,6) = lambda1;
2167 res(6,7) = lambda1;
2168 res(4,7) = lambda1;
2169 }
2170
2171 res += res.transpose().eval();
2172
2173 // not required for small Platonic solids
2174// if (VARIANT==0)
2175// {
2176// auto res_ = compress_CuthillMcKee(res);
2177// res = res_;
2178// }
2179
2180 return res;
2181}
2182
2183ArrayXXd hopping_triangular (int L, int VARIANT=0, double lambda1=1.)
2184{
2185 ArrayXXd res(L,L); res.setZero();
2186
2187 if (L==3)
2188 {
2189 res(0,1) = lambda1;
2190 res(1,2) = lambda1;
2191 res(0,2) = lambda1;
2192 }
2193 else if (L==4) // two edge-shared triangles
2194 {
2195 res(0,1) = lambda1;
2196 res(1,2) = lambda1;
2197 res(2,3) = lambda1;
2198 res(0,3) = lambda1;
2199 res(0,2) = lambda1;
2200 }
2201 else if (L==5)
2202 {
2203 res(0,1) = lambda1;
2204 res(1,2) = lambda1;
2205 res(0,2) = lambda1;
2206
2207 res(2,3) = lambda1;
2208 res(2,4) = lambda1;
2209 res(3,4) = lambda1;
2210 }
2211 else if (L==6)
2212 {
2213 res(0,1) = lambda1;
2214 res(1,2) = lambda1;
2215 res(0,2) = lambda1;
2216
2217 res(2,3) = lambda1;
2218 res(2,4) = lambda1;
2219 res(3,4) = lambda1;
2220
2221 res(1,4) = lambda1;
2222 res(1,5) = lambda1;
2223 res(4,5) = lambda1;
2224 }
2225
2226 res += res.transpose().eval();
2227 return res;
2228}
2229
2230void add_tetrahedron (int i, int j, int k, int l, vector<pair<size_t,size_t>> &target)
2231{
2232// std::cout << "ijkl=" << i << ", " << j << ", " << k << ", " << l << std::endl;
2233//
2234 (i<j)? target.push_back(pair<size_t,size_t>(i,j)) : target.push_back(pair<size_t,size_t>(j,i));
2235 (i<k)? target.push_back(pair<size_t,size_t>(i,k)) : target.push_back(pair<size_t,size_t>(k,i));
2236 (i<l)? target.push_back(pair<size_t,size_t>(i,l)) : target.push_back(pair<size_t,size_t>(l,i));
2237 (j<k)? target.push_back(pair<size_t,size_t>(j,k)) : target.push_back(pair<size_t,size_t>(k,j));
2238 (j<l)? target.push_back(pair<size_t,size_t>(j,l)) : target.push_back(pair<size_t,size_t>(l,j));
2239 (k<l)? target.push_back(pair<size_t,size_t>(k,l)) : target.push_back(pair<size_t,size_t>(l,k));
2240//
2241// for (int n=0; n<6; ++n)
2242// {
2243// std::cout << "tetrahedron: " << target[target.size()-6+n].first << ", " << target[target.size()-6+n].second << std::endl;
2244// }
2245// std::cout << std::endl;
2246}
2247
2248void add_edge (int i, int j, vector<pair<size_t,size_t>> &target)
2249{
2250 (i<j)? target.push_back(pair<size_t,size_t>(i,j)) : target.push_back(pair<size_t,size_t>(j,i));
2251}
2252
2253ArrayXXd hopping_sodaliteCage (int L=60, int VARIANT=0, double lambda1=1.)
2254{
2255 std::vector<std::pair<std::size_t, std::size_t>> edges;
2256
2257 ArrayXXd res(L,L); res.setZero();
2258
2259 if (L==60)
2260 {
2261 if (VARIANT==0)
2262 {
2263 // dmax=16 with Cuthill-McKee algorithm
2264 add_tetrahedron(0,1,3,4,edges); // dmax=4
2265 add_tetrahedron(0,2,5,6,edges); // dmax=6
2266 add_tetrahedron(5,13,14,15,edges); // dmax=10
2267 add_tetrahedron(14,20,27,28,edges); // dmax=14
2268 add_tetrahedron(8,19,20,21,edges); // dmax=13
2269 add_tetrahedron(3,7,8,9,edges); // dmax=6
2270
2271 add_tetrahedron(9,11,22,23,edges); // dmax=14
2272 add_tetrahedron(4,10,11,12,edges); // dmax=8
2273 add_tetrahedron(6,16,17,18,edges); // dmax=12
2274 add_tetrahedron(15,18,29,30,edges); // dmax=15
2275 add_tetrahedron(28,35,42,43,edges); // dmax=15
2276 add_tetrahedron(21,33,34,35,edges); // dmax=14
2277
2278 add_tetrahedron(34,37,49,50,edges); // dmax=16
2279 add_tetrahedron(23,36,37,38,edges); // dmax=15
2280 add_tetrahedron(12,24,25,26,edges); // dmax=14
2281 add_tetrahedron(17,25,31,32,edges); // dmax=15
2282 add_tetrahedron(30,44,45,46,edges); // dmax=16
2283 add_tetrahedron(43,45,53,54,edges); // dmax=11
2284
2285 add_tetrahedron(54,56,58,59,edges); // dmax=5
2286 add_tetrahedron(50,52,57,58,edges); // dmax=8
2287 add_tetrahedron(38,41,51,52,edges); // dmax=14
2288 add_tetrahedron(26,39,40,41,edges); // dmax=15
2289 add_tetrahedron(32,40,47,48,edges); // dmax=16
2290 add_tetrahedron(46,48,55,56,edges); // dmax=10
2291 }
2292 else if (VARIANT==1)
2293 {
2294 // dmax=25 without optimization
2295 add_tetrahedron(0,5,6,13,edges); // dmax=13
2296 add_tetrahedron(0,1,7,14,edges); // dmax=14
2297 add_tetrahedron(1,2,8,15,edges); // dmax=14
2298 add_tetrahedron(2,3,9,16,edges); // dmax=14
2299 add_tetrahedron(3,4,10,17,edges); // dmax=14
2300 add_tetrahedron(4,5,11,12,edges); // dmax=8
2301
2302 add_tetrahedron(12,20,21,33,edges); // dmax=21
2303 add_tetrahedron(13,21,22,34,edges); // dmax=21
2304 add_tetrahedron(14,24,25,37,edges); // dmax=23
2305 add_tetrahedron(15,25,26,38,edges); // dmax=23
2306 add_tetrahedron(16,28,29,41,edges); // dmax=25
2307 add_tetrahedron(17,18,29,30,edges); // dmax=13
2308
2309 add_tetrahedron(18,19,31,43,edges); // dmax=25
2310 add_tetrahedron(19,20,32,44,edges); // dmax=25
2311 add_tetrahedron(22,23,35,45,edges); // dmax=23
2312 add_tetrahedron(23,24,36,46,edges); // dmax=23
2313 add_tetrahedron(26,27,39,47,edges); // dmax=21
2314 add_tetrahedron(27,28,40,42,edges); // dmax=15
2315
2316 add_tetrahedron(42,57,58,59,edges); // dmax=17
2317 add_tetrahedron(43,55,56,57,edges); // dmax=14
2318 add_tetrahedron(44,53,54,55,edges); // dmax=11
2319 add_tetrahedron(45,51,52,53,edges); // dmax=8
2320 add_tetrahedron(46,49,50,51,edges); // dmax=5
2321 add_tetrahedron(47,48,49,59,edges); // dmax=12
2322 }
2323 }
2324 else if (L==50)
2325 {
2326 if (VARIANT==0)
2327 {
2328 add_tetrahedron(0,2,3,6,edges);
2329 add_tetrahedron(3,8,7,9,edges);
2330 add_tetrahedron(8,11,19,20,edges);
2331 add_tetrahedron(4,10,11,12,edges);
2332 add_tetrahedron(0,1,4,5,edges);
2333
2334 add_tetrahedron(6,16,17,18,edges);
2335 add_tetrahedron(9,21,22,23,edges);
2336 add_tetrahedron(20,33,34,35,edges);
2337 add_tetrahedron(12,24,25,26,edges);
2338 add_tetrahedron(5,13,14,15,edges);
2339
2340 add_tetrahedron(15,17,29,30,edges);
2341 add_tetrahedron(18,22,31,32,edges);
2342 add_tetrahedron(23,34,36,37,edges);
2343 add_tetrahedron(25,35,38,39,edges);
2344 add_tetrahedron(14,26,27,28,edges);
2345
2346 add_tetrahedron(30,42,43,44,edges);
2347 add_tetrahedron(32,44,45,46,edges);
2348 add_tetrahedron(37,46,47,48,edges);
2349 add_tetrahedron(39,41,48,49,edges);
2350 add_tetrahedron(28,40,41,42,edges);
2351 }
2352 else if (VARIANT==1)
2353 {
2354 add_tetrahedron(0,1,6,12,edges);
2355 add_tetrahedron(1,2,7,13,edges);
2356 add_tetrahedron(2,3,8,14,edges);
2357 add_tetrahedron(3,4,9,10,edges);
2358 add_tetrahedron(0,4,5,11,edges);
2359
2360 add_tetrahedron(12,21,31,32,edges);
2361 add_tetrahedron(13,23,33,34,edges);
2362 add_tetrahedron(14,15,25,26,edges);
2363 add_tetrahedron(10,17,27,28,edges);
2364 add_tetrahedron(11,19,29,30,edges);
2365
2366 add_tetrahedron(20,30,31,38,edges);
2367 add_tetrahedron(22,32,33,39,edges);
2368 add_tetrahedron(24,25,34,35,edges);
2369 add_tetrahedron(16,26,27,36,edges);
2370 add_tetrahedron(18,28,29,37,edges);
2371
2372 add_tetrahedron(38,44,45,49,edges);
2373 add_tetrahedron(39,40,45,46,edges);
2374 add_tetrahedron(35,41,46,47,edges);
2375 add_tetrahedron(36,42,47,48,edges);
2376 add_tetrahedron(37,43,48,49,edges);
2377 }
2378 }
2379 else if (L==20)
2380 {
2381 if (VARIANT==0)
2382 {
2383 add_tetrahedron(0,1,4,5,edges);
2384 add_tetrahedron(0,2,3,6,edges);
2385 add_tetrahedron(3,7,8,9,edges);
2386 add_tetrahedron(8,4,10,11,edges);
2387
2388 add_tetrahedron(5,13,12,14,edges);
2389 add_tetrahedron(6,14,15,16,edges);
2390 add_tetrahedron(9,16,17,18,edges);
2391 add_tetrahedron(11,13,18,19,edges);
2392 }
2393 else if (VARIANT==1)
2394 {
2395 add_tetrahedron(0,3,4,12,edges);
2396 add_tetrahedron(0,1,5,13,edges);
2397 add_tetrahedron(1,2,6,14,edges);
2398 add_tetrahedron(2,3,7,15,edges);
2399
2400 add_tetrahedron(4,8,9,16,edges);
2401 add_tetrahedron(5,9,10,17,edges);
2402 add_tetrahedron(6,10,11,18,edges);
2403 add_tetrahedron(7,8,11,19,edges);
2404 }
2405 }
2406 else if (L==16)
2407 {
2408 add_tetrahedron(1,2,3,4,edges);
2409 add_tetrahedron(4,6,7,8,edges);
2410 add_tetrahedron(8,9,10,11,edges);
2411 add_tetrahedron(3,9,14,15,edges);
2412
2413 add_edge(0,1,edges);
2414 add_edge(5,6,edges);
2415 add_edge(11,12,edges);
2416 add_edge(13,14,edges);
2417 }
2418 else if (L==28)
2419 {
2420 add_tetrahedron(0,1,2,3,edges);
2421 add_tetrahedron(3,4,5,6,edges);
2422 add_tetrahedron(7,8,9,10,edges);
2423 add_tetrahedron(10,11,12,13,edges);
2424 add_tetrahedron(14,15,16,17,edges);
2425 add_tetrahedron(17,18,19,20,edges);
2426 add_tetrahedron(21,22,23,24,edges);
2427 add_tetrahedron(24,25,26,27,edges);
2428
2429 add_edge(0,26,edges);
2430 add_edge(1,27,edges);
2431 add_edge(5,8,edges);
2432 add_edge(6,7,edges);
2433 add_edge(12,15,edges);
2434 add_edge(13,14,edges);
2435 add_edge(19,22,edges);
2436 add_edge(20,21,edges);
2437 }
2438
2439 for (int e=0; e<edges.size(); ++e)
2440 {
2441 int i = edges[e].first;
2442 int j = edges[e].second;
2443 res(i,j) = lambda1;
2444 }
2445
2446 res += res.transpose().eval();
2447
2448 if (VARIANT==0)
2449 {
2450 compress_CuthillMcKee(res,true);
2451 }
2452
2453 return res;
2454}
2455
2456ArrayXXd hopping_Mn32 (double lambda_cap=1., double lambda_corner=0., double lambda_edge=1., int VARIANT=0)
2457{
2458 std::vector<std::pair<std::size_t, std::size_t>> edges;
2459
2460 ArrayXXd res(32,32); res.setZero();
2461
2462 // outer circle:
2463 add_tetrahedron(0,1,2,3,edges);
2464 add_tetrahedron(4,5,6,7,edges);
2465 add_tetrahedron(8,9,10,11,edges);
2466 add_tetrahedron(12,13,14,15,edges);
2467 // inner circle:
2468 add_tetrahedron(16,17,18,19,edges);
2469 add_tetrahedron(20,21,22,23,edges);
2470 add_tetrahedron(24,25,26,27,edges);
2471 add_tetrahedron(28,29,30,31,edges);
2472
2473 vector<int> caps = {3,7,11,15, 19,23,27,31};
2474
2475 for (int e=0; e<edges.size(); ++e)
2476 {
2477 int i = edges[e].first;
2478 int j = edges[e].second;
2479
2480 auto it_i = find(caps.begin(), caps.end(), i);
2481 auto it_j = find(caps.begin(), caps.end(), j);
2482
2483 if (it_i!=caps.end() or it_j!=caps.end())
2484 {
2485 res(i,j) = lambda_cap;
2486 }
2487 else
2488 {
2489 res(i,j) = lambda_corner;
2490 }
2491 }
2492
2493 edges.clear();
2494
2495 // outer circle:
2496 edges.push_back(pair<size_t,size_t>(1,4));
2497 edges.push_back(pair<size_t,size_t>(5,8));
2498 edges.push_back(pair<size_t,size_t>(9,12));
2499 edges.push_back(pair<size_t,size_t>(0,13));
2500 // connection:
2501 edges.push_back(pair<size_t,size_t>(2,16));
2502 edges.push_back(pair<size_t,size_t>(6,21));
2503 edges.push_back(pair<size_t,size_t>(10,26));
2504 edges.push_back(pair<size_t,size_t>(14,29));
2505 // inner circle:
2506 edges.push_back(pair<size_t,size_t>(17,20));
2507 edges.push_back(pair<size_t,size_t>(22,24));
2508 edges.push_back(pair<size_t,size_t>(25,28));
2509 edges.push_back(pair<size_t,size_t>(18,30));
2510
2511 for (int e=0; e<edges.size(); ++e)
2512 {
2513 int i = edges[e].first;
2514 int j = edges[e].second;
2515 res(i,j) = lambda_edge;
2516 }
2517
2518 res += res.transpose().eval();
2519
2520 if (VARIANT==0)
2521 {
2522 compress_CuthillMcKee(res,true);
2523 }
2524
2525 return res;
2526}
2527
2528pair<ArrayXXd,vector<SUB_LATTICE> > hopping_PPV (int L, int VARIANT=0, double t0=1., double tsingle=1., double tdouble=1., string BC="")
2529{
2530 ArrayXXd res(L,L); res.setZero();
2531 assert(L==8);
2532 vector<SUB_LATTICE> G;
2533
2534 if (BC == "INIT_HEX")
2535 {
2536 res(0,1) = t0;
2537 res(0,2) = t0;
2538 res(1,3) = t0;
2539 res(2,4) = t0;
2540 res(3,5) = t0;
2541 res(4,5) = t0;
2542 res(5,6) = tsingle; // tvinyl(single)
2543 res(6,7) = tdouble; // tvinyl(double)
2544 // intermolecular hopping is is tsingle
2545
2546 // AABA|BBAB
2547 G.push_back(static_cast<SUB_LATTICE>(1)); //0
2548 G.push_back(static_cast<SUB_LATTICE>(-1)); //1
2549 G.push_back(static_cast<SUB_LATTICE>(-1)); //2
2550 G.push_back(static_cast<SUB_LATTICE>(1)); //3
2551 G.push_back(static_cast<SUB_LATTICE>(1)); //4
2552 G.push_back(static_cast<SUB_LATTICE>(-1)); //5
2553 G.push_back(static_cast<SUB_LATTICE>(1)); //6
2554 G.push_back(static_cast<SUB_LATTICE>(-1)); //7
2555 }
2556 else // PPV
2557 {
2558 res(0,1) = tsingle;
2559 res(1,2) = t0;
2560 res(1,3) = t0;
2561 res(2,4) = t0;
2562 res(3,5) = t0;
2563 res(4,6) = t0;
2564 res(5,6) = t0;
2565 res(6,7) = tsingle; // tvinyl(single)
2566 // intermolecular hopping is is tdouble
2567
2568 G.push_back(static_cast<SUB_LATTICE>(1)); //0
2569 G.push_back(static_cast<SUB_LATTICE>(-1)); //1
2570 G.push_back(static_cast<SUB_LATTICE>(1)); //2
2571 G.push_back(static_cast<SUB_LATTICE>(1)); //3
2572 G.push_back(static_cast<SUB_LATTICE>(-1)); //4
2573 G.push_back(static_cast<SUB_LATTICE>(-1)); //5
2574 G.push_back(static_cast<SUB_LATTICE>(1)); //6
2575 G.push_back(static_cast<SUB_LATTICE>(-1)); //7
2576 }
2577
2578 res += res.transpose().eval();
2579
2580 if (VARIANT==0)
2581 {
2582 CuthillMcKeeCompressor CMK(res,true);
2583 CMK.apply_compression(res);
2584
2585 vector<SUB_LATTICE> G_ = G;
2586 for (int i=0; i<L; ++i)
2587 {
2588 G[CMK.get_transform()[i]] = G_[i];
2589 }
2590 }
2591
2592 pair<ArrayXXd,vector<SUB_LATTICE> > ret(res,G);
2593 return ret;
2594}
2595
2596pair<ArrayXXd,vector<SUB_LATTICE> > hopping_triangulene (int L, int VARIANT=0, double lambda=1., double lambda2=1., string BC="")
2597{
2598 ArrayXXd res(L,L); res.setZero();
2599 assert(L==13 or L==8 or L==10 or L==22 or L==4); // 33, 46, 61 // 4: simplified model, 8: simplified, larger unit cell with different coupling
2600 // (6+)7+9+11+13
2601
2602 vector<SUB_LATTICE> G;
2603
2604 if (L==4)
2605 {
2606 res(0,2) = lambda;
2607 res(1,2) = lambda;
2608 res(2,3) = lambda;
2609
2610 G.push_back(static_cast<SUB_LATTICE>(1)); //0
2611 G.push_back(static_cast<SUB_LATTICE>(1)); //1
2612 G.push_back(static_cast<SUB_LATTICE>(-1)); //2
2613 G.push_back(static_cast<SUB_LATTICE>(1)); //3
2614 }
2615 else if (L==8)
2616 {
2617 if (BC == "CUT_HEX")
2618 {
2619 res(0,2) = lambda;
2620 res(1,2) = lambda;
2621 res(2,3) = lambda;
2622 res(3,4) = lambda2; // tinter
2623 res(4,5) = lambda;
2624 res(5,6) = lambda;
2625 res(5,7) = lambda;
2626 // intermolecular hopping is thex (2x)
2627
2628 // AABA|BBAB
2629 G.push_back(static_cast<SUB_LATTICE>(1)); //0
2630 G.push_back(static_cast<SUB_LATTICE>(1)); //1
2631 G.push_back(static_cast<SUB_LATTICE>(-1)); //2
2632 G.push_back(static_cast<SUB_LATTICE>(1)); //3
2633 G.push_back(static_cast<SUB_LATTICE>(-1)); //4
2634 G.push_back(static_cast<SUB_LATTICE>(1)); //5
2635 G.push_back(static_cast<SUB_LATTICE>(-1)); //6
2636 G.push_back(static_cast<SUB_LATTICE>(-1)); //7
2637 }
2638 else
2639 {
2640 res(0,1) = lambda;
2641 res(1,2) = lambda;
2642 res(1,3) = lambda;
2643 res(2,4) = lambda2; // thex
2644 res(3,5) = lambda2; // thex
2645 res(4,6) = lambda;
2646 res(5,6) = lambda;
2647 res(6,7) = lambda;
2648 // intermolecular hopping is tinter
2649
2650 G.push_back(static_cast<SUB_LATTICE>(1)); //0
2651 G.push_back(static_cast<SUB_LATTICE>(-1)); //1
2652 G.push_back(static_cast<SUB_LATTICE>(1)); //2
2653 G.push_back(static_cast<SUB_LATTICE>(1)); //3
2654 G.push_back(static_cast<SUB_LATTICE>(-1)); //4
2655 G.push_back(static_cast<SUB_LATTICE>(-1)); //5
2656 G.push_back(static_cast<SUB_LATTICE>(1)); //6
2657 G.push_back(static_cast<SUB_LATTICE>(-1)); //7
2658 }
2659 }
2660 else if (L==10) //porphyrine coarse-grained
2661 {
2662 res(0,1) = lambda;
2663
2664 res(1,2) = lambda;
2665 res(1,3) = lambda;
2666
2667 res(2,4) = lambda;
2668 res(4,6) = lambda;
2669
2670 res(3,5) = lambda;
2671 res(5,7) = lambda;
2672
2673 res(6,8) = lambda;
2674 res(7,8) = lambda;
2675
2676 res(8,9) = lambda;
2677
2678 G.push_back(static_cast<SUB_LATTICE>(1)); //0
2679 G.push_back(static_cast<SUB_LATTICE>(-1)); //1
2680 G.push_back(static_cast<SUB_LATTICE>(1)); //2
2681 G.push_back(static_cast<SUB_LATTICE>(1)); //3
2682
2683 G.push_back(static_cast<SUB_LATTICE>(-1)); //4
2684 G.push_back(static_cast<SUB_LATTICE>(-1)); //5
2685
2686 G.push_back(static_cast<SUB_LATTICE>(1)); //6
2687 G.push_back(static_cast<SUB_LATTICE>(1)); //7
2688 G.push_back(static_cast<SUB_LATTICE>(-1)); //8
2689 G.push_back(static_cast<SUB_LATTICE>(1)); //9
2690 }
2691 else if (L==13)
2692 {
2693 for (int i=0; i<=10; ++i) res(i,i+1) = lambda;
2694 res(0,11) = lambda;
2695
2696 res(6,12) = lambda;
2697 res(2,12) = lambda;
2698 res(10,12) = lambda;
2699 }
2700 else if (L==22) // for triangulene proper, lambda2 is the NNN hopping t3
2701 {
2702 res(11,15) = lambda;
2703 res(15,19) = lambda;
2704
2705 res(8,11) = lambda;
2706 res(16,19) = lambda;
2707
2708 res(5,8) = lambda;
2709 res(8,12) = lambda;
2710 res(12,16) = lambda;
2711 res(16,20) = lambda;
2712
2713 res(3,5) = lambda;
2714 res(9,12) = lambda;
2715 res(17,20) = lambda;
2716
2717 res(1,3) = lambda;
2718 res(3,6) = lambda;
2719 res(6,9) = lambda;
2720 res(9,13) = lambda;
2721 res(13,17) = lambda;
2722 res(17,21) = lambda;
2723
2724 res(0,1) = lambda;
2725 res(4,6) = lambda;
2726 res(10,13) = lambda;
2727 res(18,21) = lambda;
2728
2729 res(0,2) = lambda;
2730 res(2,4) = lambda;
2731 res(4,7) = lambda;
2732 res(7,10) = lambda;
2733 res(10,14) = lambda;
2734 res(14,18) = lambda;
2735
2736 G.push_back(static_cast<SUB_LATTICE>(1)); //0
2737 G.push_back(static_cast<SUB_LATTICE>(-1)); //1
2738 G.push_back(static_cast<SUB_LATTICE>(-1)); //2
2739 G.push_back(static_cast<SUB_LATTICE>(1)); //3
2740 G.push_back(static_cast<SUB_LATTICE>(1)); //4
2741 G.push_back(static_cast<SUB_LATTICE>(-1)); //5
2742 G.push_back(static_cast<SUB_LATTICE>(-1)); //6
2743 G.push_back(static_cast<SUB_LATTICE>(-1)); //7
2744 G.push_back(static_cast<SUB_LATTICE>(1)); //8
2745 G.push_back(static_cast<SUB_LATTICE>(1)); //9
2746 G.push_back(static_cast<SUB_LATTICE>(1)); //10
2747 G.push_back(static_cast<SUB_LATTICE>(-1)); //11
2748 G.push_back(static_cast<SUB_LATTICE>(-1)); //12
2749 G.push_back(static_cast<SUB_LATTICE>(-1)); //13
2750 G.push_back(static_cast<SUB_LATTICE>(-1)); //14
2751 G.push_back(static_cast<SUB_LATTICE>(1)); //15
2752 G.push_back(static_cast<SUB_LATTICE>(1)); //16
2753 G.push_back(static_cast<SUB_LATTICE>(1)); //17
2754 G.push_back(static_cast<SUB_LATTICE>(1)); //18
2755 G.push_back(static_cast<SUB_LATTICE>(-1)); //19
2756 G.push_back(static_cast<SUB_LATTICE>(-1)); //20
2757 G.push_back(static_cast<SUB_LATTICE>(-1)); //21
2758
2759 res(12,15) = lambda2;
2760 res(11,16) = lambda2;
2761 res(8,19) = lambda2;
2762
2763 res(5,9) = lambda2;
2764 res(3,12) = lambda2;
2765 res(6,8) = lambda2;
2766
2767 res(13,16) = lambda2;
2768 res(12,17) = lambda2;
2769 res(9,20) = lambda2;
2770
2771 res(2,3) = lambda2;
2772 res(1,4) = lambda2;
2773 res(0,6) = lambda2;
2774
2775 res(6,10) = lambda2;
2776 res(7,9) = lambda2;
2777 res(4,13) = lambda2;
2778
2779 res(10,21) = lambda2;
2780 res(14,17) = lambda2;
2781 res(13,18) = lambda2;
2782 }
2783
2784 res += res.transpose().eval();
2785
2786 //cout << "VARIANT=" << VARIANT << endl;
2787 if (VARIANT==0)
2788 {
2789 //auto res_ = compress_CuthillMcKee(res,true);
2790 CuthillMcKeeCompressor CMK(res,true);
2791 CMK.apply_compression(res);
2792
2793 vector<SUB_LATTICE> G_ = G;
2794 for (int i=0; i<L; ++i)
2795 {
2796 G[CMK.get_transform()[i]] = G_[i];
2797 }
2798
2799// for (int i=0; i<L; ++i)
2800// {
2801// lout << "i=" << i << ", G[i]=" << G[i]<< ", orig.G[i]=" << G_[i] << endl;
2802// }
2803 }
2804
2805 pair<ArrayXXd,vector<SUB_LATTICE> > ret(res,G);
2806 return ret;
2807}
2808
2809pair<ArrayXXd,vector<SUB_LATTICE> > hopping_triangulene_dimer (int VARIANT=0, double lambda=1., double lambda2=1., string BC="")
2810{
2811 auto [T,G0] = hopping_triangulene(22,1,lambda,lambda2);
2812 ArrayXXd res(44,44); res = 0.;
2813 res.topLeftCorner(22,22) = T;
2814 res.bottomRightCorner(22,22) = T;
2815
2816 ArrayXXd coupling(44,44); coupling = 0.;
2817 coupling(2,1+22) = lambda;
2818 coupling(7,5+22) = lambda;
2819 coupling(14,11+22) = lambda;
2820
2821 coupling(2,1+22) = lambda2;
2822 coupling(2,5+22) = lambda2;
2823
2824 coupling(4,3+22) = lambda2;
2825
2826 coupling(7,1+22) = lambda2;
2827 coupling(7,11+22) = lambda2;
2828
2829 coupling(10,8+22) = lambda2;
2830
2831 coupling(14,5+22) = lambda2;
2832 coupling(14,11+22) = lambda2;
2833
2834 coupling += coupling.transpose().eval();
2835
2836 res += coupling;
2837
2838 vector<SUB_LATTICE> G;
2839 for (int i=0; i<22; ++i) G.push_back(G0[i]);
2840 for (int i=0; i<22; ++i) G.push_back(flip_sublattice(G0[i]));
2841
2842 if (VARIANT==0)
2843 {
2844 CuthillMcKeeCompressor CMK(res,true);
2845 CMK.apply_compression(res);
2846
2847 vector<SUB_LATTICE> G_ = G;
2848 for (int i=0; i<44; ++i)
2849 {
2850 G[CMK.get_transform()[i]] = G_[i];
2851 }
2852 }
2853
2854 pair<ArrayXXd,vector<SUB_LATTICE> > ret(res,G);
2855 return ret;
2856}
2857
2858ArrayXXd hopping_coronene (int L, int VARIANT=0, double lambda=1.)
2859{
2860 ArrayXXd res(L,L); res.setZero();
2861 assert(L==24);
2862
2863 if (L==24)
2864 {
2865 res(0,2) = lambda;
2866 res(0,3) = lambda;
2867 res(1,3) = lambda;
2868 res(1,4) = lambda;
2869 res(2,5) = lambda;
2870 res(3,6) = lambda;
2871 res(4,7) = lambda;
2872 res(5,8) = lambda;
2873 res(5,9) = lambda;
2874 res(6,9) = lambda;
2875 res(6,10) = lambda;
2876 res(7,10) = lambda;
2877 res(7,11) = lambda;
2878 res(8,12) = lambda;
2879 res(9,13) = lambda;
2880 res(10,14) = lambda;
2881 res(11,15) = lambda;
2882 res(12,16) = lambda;
2883 res(13,16) = lambda;
2884 res(13,17) = lambda;
2885 res(14,17) = lambda;
2886 res(14,18) = lambda;
2887 res(15,18) = lambda;
2888 res(16,19) = lambda;
2889 res(17,20) = lambda;
2890 res(18,21) = lambda;
2891 res(19,22) = lambda;
2892 res(20,22) = lambda;
2893 res(20,23) = lambda;
2894 res(21,23) = lambda;
2895 }
2896
2897 res += res.transpose().eval();
2898
2899 if (VARIANT==0)
2900 {
2901 compress_CuthillMcKee(res,true);
2902 }
2903
2904 return res;
2905}
2906
2907ArrayXXd hopping_corannulene (int L, int VARIANT=0, double lambda=1.)
2908{
2909 ArrayXXd res(L,L); res.setZero();
2910 assert(L==20);
2911
2912 if (L==20)
2913 {
2914 for (int i=0; i<=13; ++i) res(i,i+1) = lambda;
2915 res(0,14) = lambda;
2916
2917 res(13,15) = lambda;
2918 res(10,16) = lambda;
2919 res(7,17) = lambda;
2920 res(4,18) = lambda;
2921 res(1,19) = lambda;
2922
2923 for (int i=15; i<=18; ++i) res(i,i+1) = lambda;
2924 res(15,19) = lambda;
2925 }
2926
2927 res += res.transpose().eval();
2928
2929 if (VARIANT==0)
2930 {
2931 compress_CuthillMcKee(res,true);
2932 }
2933
2934 return res;
2935}
2936
2937ArrayXXd hopping_square_plaquette (int L, int VARIANT=0, double lambda=1.)
2938{
2939 ArrayXXd res(L,L); res.setZero();
2940
2941 if (L==16)
2942 {
2943 for (int i=0; i<3; ++i) res(i,i+1) = lambda;
2944 for (int i=4; i<7; ++i) res(i,i+1) = lambda;
2945 for (int i=8; i<11; ++i) res(i,i+1) = lambda;
2946 for (int i=12; i<15; ++i) res(i,i+1) = lambda;
2947
2948 for (int i=0; i<=2; ++i) res(4*i,4*i+4) = lambda;
2949 for (int i=0; i<=2; ++i) res(1+4*i,1+4*i+4) = lambda;
2950 for (int i=0; i<=2; ++i) res(2+4*i,2+4*i+4) = lambda;
2951 for (int i=0; i<=2; ++i) res(3+4*i,3+4*i+4) = lambda;
2952
2953// res(0,3) = lambda;
2954// res(4,7) = lambda;
2955// res(8,11) = lambda;
2956// res(12,15) = lambda;
2957//
2958// res(0,12) = lambda;
2959// res(1,13) = lambda;
2960// res(2,14) = lambda;
2961// res(3,15) = lambda;
2962 }
2963 if (L==20)
2964 {
2965 for (int i=0; i<3; ++i) res(i,i+1) = lambda;
2966 for (int i=4; i<7; ++i) res(i,i+1) = lambda;
2967 for (int i=8; i<11; ++i) res(i,i+1) = lambda;
2968 for (int i=12; i<15; ++i) res(i,i+1) = lambda;
2969 for (int i=16; i<19; ++i) res(i,i+1) = lambda;
2970
2971 for (int i=0; i<=3; ++i) res(4*i,4*i+4) = lambda;
2972 for (int i=0; i<=3; ++i) res(1+4*i,1+4*i+4) = lambda;
2973 for (int i=0; i<=3; ++i) res(2+4*i,2+4*i+4) = lambda;
2974 for (int i=0; i<=3; ++i) res(3+4*i,3+4*i+4) = lambda;
2975
2976 res(0,3) = lambda;
2977 res(4,7) = lambda;
2978 res(8,11) = lambda;
2979 res(12,15) = lambda;
2980 res(16,19) = lambda;
2981
2982 res(0,16) = lambda;
2983 res(1,17) = lambda;
2984 res(2,18) = lambda;
2985 res(3,19) = lambda;
2986 }
2987
2988 res += res.transpose().eval();
2989
2990 if (VARIANT==0)
2991 {
2992 compress_CuthillMcKee(res,true);
2993 }
2994
2995 return res;
2996}
2997
2998// Calculates distance matrix from adjacency matrix of a graph
2999// The distance matrix has 1 for nearest neighbours, 2 for next-nearest neighbours etc.
3000ArrayXXi calc_distanceMatrix (ArrayXXd adjacencyMatrix)
3001{
3002 int L = adjacencyMatrix.rows();
3003 assert(adjacencyMatrix.cols() == L);
3004
3005 ArrayXXi dist0(L,L); dist0 = 0;
3006 dist0.matrix().diagonal().setConstant(1);
3007
3008 ArrayXXi dist1(L,L); dist1 = 0;
3009 for (int i=0; i<L; ++i)
3010 for (int j=0; j<L; ++j)
3011 {
3012 if (abs(adjacencyMatrix(i,j)) > 0.)
3013 {
3014 dist1(i,j) = 1;
3015 }
3016 }
3017
3018 ArrayXXi next = dist1;
3019
3020 // each entry is a matrix with 1 for the given distance, 0 otherwise
3021 vector<ArrayXXi> dist;
3022 dist.push_back(dist0);
3023 dist.push_back(dist1);
3024
3025 while (next.sum() != 0)
3026 {
3027 // calculate d-th power of adjacency matrix
3028 next = next.matrix()*dist1.matrix();
3029
3030 // remove already known distances (multiply array-wise with 0)
3031 for (int d=0; d<dist.size(); ++d) next = next*(1-dist[d]);
3032 // set non-zero entries to 1
3033 for (int i=0; i<L; ++i) for (int j=0; j<L; ++j) if (next(i,j)>0) next(i,j) = 1;
3034
3035 dist.push_back(next);
3036 }
3037
3038 ArrayXXi res(L,L); res = 0;
3039 for (int d=0; d<dist.size(); ++d) res = res+d*dist[d];
3040 return res;
3041}
3042
3043ArrayXXi hopping_MnRing (int Ncells, double J1, double J2, double J3, double J4=0., double J5=0., double J6=0., double J7=0.)
3044{
3045 int L = 7*Ncells;
3046 ArrayXXi res(L,L);
3047
3048 for (int ic=0; ic<Ncells; ++ic)
3049 {
3050 int icell = ic*7;
3051
3052 res((0+icell)%L,(1+icell)%L) = 1.;
3053 res((1+icell)%L,(2+icell)%L) = 1.;
3054
3055 res((0+icell)%L,(3+icell)%L) = J2/J1;
3056 res((1+icell)%L,(3+icell)%L) = J2/J1;
3057 res((1+icell)%L,(4+icell)%L) = J2/J1;
3058 res((2+icell)%L,(4+icell)%L) = J2/J1;
3059
3060 res((3+icell)%L,(4+icell)%L) = J3/J1;
3061 res((3+icell)%L,(5+icell)%L) = J3/J1;
3062 res((3+icell)%L,(6+icell)%L) = J3/J1;
3063 res((4+icell)%L,(5+icell)%L) = J3/J1;
3064 res((5+icell)%L,(6+icell)%L) = J3/J1;
3065
3066 res((5+icell)%L,(7+icell)%L) = J2/J1;
3067 res((5+icell)%L,(8+icell)%L) = J2/J1;
3068 res((6+icell)%L,(8+icell)%L) = J2/J1;
3069 res((6+icell)%L,(9+icell)%L) = J2/J1;
3070
3071 res((0+icell)%L,(4+icell)%L) = J4/J1;
3072 res((2+icell)%L,(3+icell)%L) = J4/J1;
3073 res((1+icell)%L,(6+icell)%L) = J4/J1;
3074
3075 res((5+icell)%L,(9+icell)%L) = J4/J1;
3076 res((6+icell)%L,(7+icell)%L) = J4/J1;
3077 res((4+icell)%L,(8+icell)%L) = J4/J1;
3078 }
3079
3080 res += res.transpose().eval();
3081 compress_CuthillMcKee(res,true);
3082
3083 return res;
3084}
3085
3086void push_back_KondoUnpacked (vector<Param> &params, size_t L, double J, double t, size_t D, bool START_WITH_SPIN=true)
3087{
3088 int SPIN_PARITY = (START_WITH_SPIN==true)? 0:1;
3089 for (size_t l=0; l<2*L; ++l)
3090 {
3091 // spin site
3092 if (l%2 == SPIN_PARITY)
3093 {
3094 params.push_back({"D",D,l});
3095 params.push_back({"LyF",0ul,l});
3096 if (START_WITH_SPIN==true)
3097 {
3098 params.push_back({"Inext",J,l});
3099 params.push_back({"Iprev",0.,l});
3100 }
3101 else
3102 {
3103 params.push_back({"Inext",0.,l});
3104 params.push_back({"Iprev",J,l});
3105 }
3106 params.push_back({"tPrime",0.,l});
3107 }
3108 // fermionic site
3109 else
3110 {
3111 params.push_back({"D",1ul,l});
3112 params.push_back({"LyF",1ul,l});
3113 params.push_back({"tPrime",t,l});
3114 params.push_back({"Inext",0.,l});
3115 params.push_back({"Iprev",0.,l});
3116 }
3117 }
3118}
3119
3120//Eigen::ArrayXXd coupling2d (double coupl_x, double coupl_y, size_t Lx, size_t Ly, bool PERIODIC_Y=false)
3121//{
3122// Eigen::ArrayXXd out(Lx*Ly,Lx*Ly); out.setZero();
3123// for (int ix=0; ix<Lx; ix++)
3124// for (int iy=0; iy<Ly; iy++)
3125// for (int jx=0; jx<Lx; jx++)
3126// for (int jy=0; jy<Ly; jy++)
3127// {
3128// if (abs(ix-jx) == 1 and (iy==jy))
3129// {
3130// out(ix*Ly+iy, jx*Ly+jy) += coupl_x;
3131// }
3132// else if (abs(iy-jy) == 1 and (ix==jx))
3133// {
3134// out(ix*Ly+iy, jx*Ly+jy) += coupl_y;
3135// }
3136// else if (abs(iy-jy) == Ly-1 and (ix==jx) and PERIODIC_Y == true)
3137// {
3138// out(ix*Ly+iy, jx*Ly+jy) += coupl_y;
3139// }
3140// }
3141// return out;
3142//}
3143
3144//Eigen::ArrayXXd coupling2d (double coupl, size_t Lx, size_t Ly, bool PERIODIC_Y=false)
3145//{
3146// return coupling2d (coupl,coupl,Lx,Ly,PERIODIC_Y);
3147//}
3148
3149//Eigen::ArrayXXd coupling2d_snake (double coupl_x, double coupl_y, size_t Lx, size_t Ly, bool PERIODIC_Y=false)
3150//{
3151// Eigen::ArrayXXd out(Lx*Ly,Lx*Ly); out.setZero();
3152//
3153// // Mirrors the y coordinate to create a snake.
3154// auto mirror = [&Ly] (int iy) -> int
3155// {
3156// vector<int> v(Ly);
3157// iota (begin(v),end(v),0);
3158// reverse(v.begin(),v.end());
3159// return v[iy];
3160// };
3161//
3162// for (int ix=0; ix<Lx; ix++)
3163// for (int iy=0; iy<Ly; iy++)
3164// for (int jx=0; jx<Lx; jx++)
3165// for (int jy=0; jy<Ly; jy++)
3166// {
3167// // mirror even y only
3168// int iy_ = (ix%2==0)? iy : mirror(iy);
3169// int jy_ = (jx%2==0)? jy : mirror(jy);
3170// int index_i = iy_+Ly*ix;
3171// int index_j = jy_+Ly*jx;
3172//
3177//
3178// if (abs(ix-jx) == 1 and (iy==jy))
3179// {
3180// out(index_i,index_j) += coupl_x;
3181// }
3182// else if (abs(iy-jy) == 1 and (ix==jx))
3183// {
3184// out(index_i,index_j) += coupl_y;
3185// }
3186// else if (abs(iy-jy) == Ly-1 and (ix==jx) and PERIODIC_Y == true)
3187// {
3188// out(index_i,index_j) += coupl_y;
3189// }
3190// }
3191// return out;
3192//}
3193
3194//Eigen::ArrayXXd coupling2d_snake (double coupl, size_t Lx, size_t Ly, bool PERIODIC_Y=false)
3195//{
3196// return coupling2d_snake (coupl,coupl,Lx,Ly,PERIODIC_Y);
3197//}
3198
3200//double sigma_hop (const Eigen::ArrayXXd &thop)
3201//{
3202// double res = 0.;
3203//
3204// Eigen::ArrayXd x(thop.rows()); x = 0;
3205//
3206// for (int i=0; i<thop.rows(); ++i)
3207// for (int j=0; j<thop.cols(); ++j)
3208// {
3209// if (thop(i,j) != 0.) x(i) += abs(j-i);
3210// }
3211//
3212// double avg = x.sum()/x.rows();
3213// double var = ((x-avg)*(x-avg)).sum();
3214//
3215// return sqrt(var);
3216//}
3217
3218// Warning: changed arguments to L,Ly,maxPower from Ly,maxPower on 2024-01-18
3219vector<Param> Tinf_params_fermions (size_t L, size_t Ly, size_t maxPower=1ul)
3220{
3221 vector<Param> res;
3222 res.push_back({"Ly",Ly});
3223 res.push_back({"maxPower",maxPower});
3224 res.push_back({"OPEN_BC",true});
3225 if (Ly == 2ul)
3226 {
3227 res.push_back({"t",0.});
3228 res.push_back({"tRung",1.});
3229 }
3230 else
3231 {
3232 res.push_back({"t",1.,0});
3233 res.push_back({"t",0.,1});
3234 }
3235 vector<SUB_LATTICE> G;
3236 for (int l=0; l<L; l+=4)
3237 {
3238 G.push_back(static_cast<SUB_LATTICE>(1));
3239 G.push_back(static_cast<SUB_LATTICE>(-1));
3240 G.push_back(static_cast<SUB_LATTICE>(-1));
3241 G.push_back(static_cast<SUB_LATTICE>(1));
3242 }
3243 res.push_back({"G",G});
3244 //for (size_t l=0; l<L; ++l) {lout << G[l];}
3245 //lout << endl;
3246 return res;
3247}
3248
3249vector<Param> Tinf_params_spins (size_t L, size_t Ly, size_t maxPower=1ul, size_t DA=2ul, size_t DB=2ul, bool SOFT=false)
3250{
3251 vector<Param> res;
3252 res.push_back({"Ly",Ly});
3253 res.push_back({"maxPower",maxPower});
3254 if (SOFT and DA==3ul and DB==3ul)
3255 {
3256 size_t D = 3ul;
3257 if (Ly==2)
3258 {
3259 for (size_t l=1; l<L-1; ++l)
3260 {
3261 res.push_back({"D",D,l});
3262 }
3263 res.push_back({"D",2ul,0});
3264 res.push_back({"D",2ul,L-1});
3265 }
3266 else if (Ly==1)
3267 {
3268 for (size_t l=2; l<2*L-2; ++l)
3269 {
3270 res.push_back({"D",D,l});
3271 }
3272 res.push_back({"D",2ul,0});
3273 res.push_back({"D",2ul,1});
3274 res.push_back({"D",2ul,L-2});
3275 res.push_back({"D",2ul,L-1});
3276 }
3277 }
3278 else
3279 {
3280 if (Ly == 2ul)
3281 {
3282 for (size_t l=0; l<L; l+=2)
3283 {
3284 res.push_back({"D",DA,l});
3285 res.push_back({"D",DB,l+1});
3286 }
3287 }
3288 else
3289 {
3290 for (size_t l=0; l<L; l+=4)
3291 {
3292 res.push_back({"D",DA,l});
3293 res.push_back({"D",DA,l+1});
3294 res.push_back({"D",DB,l+2});
3295 res.push_back({"D",DB,l+3});
3296 }
3297 }
3298 }
3299 if (Ly == 2ul)
3300 {
3301 res.push_back({"J",0.});
3302 res.push_back({"Jrung",1.});
3303 }
3304 else
3305 {
3306 res.push_back({"J",1.,0});
3307 res.push_back({"J",0.,1});
3308 }
3309 return res;
3310}
3311
3312inline double conjIfImag (double x) {return x;}
3313inline std::complex<double> conjIfImag (std::complex<double> x) {return conj(x);}
3314
3315template<typename Scalar>
3316Array<Scalar,Dynamic,Dynamic> hopping_PAM (int L, Scalar tfc, Scalar tcc, Scalar tff, Scalar tx, Scalar ty, bool PBC=false)
3317{
3318 Array<Scalar,Dynamic,Dynamic> t1site(2,2); t1site = 0;
3319 t1site(0,1) = tfc;
3320
3321 // L: Anzahl der physikalischen fc-Sites
3322 Array<Scalar,Dynamic,Dynamic> res(2*L,2*L); res = 0;
3323
3324 for (int l=0; l<L; ++l)
3325 {
3326 res.block(2*l,2*l, 2,2) = t1site;
3327 }
3328
3329 for (int l=0; l<L-1; ++l)
3330 {
3331 res(2*l, 2*l+2) = tcc;
3332 res(2*l+1, 2*l+3) = tff;
3333 res(2*l+1, 2*l+2) = tx;
3334 res(2*l, 2*l+3) = conjIfImag(ty);
3335 }
3336
3337// cout << "OBC:" << endl;
3338// cout << endl << res << endl << endl;
3339 if (PBC)
3340 {
3341 for (int l=0; l<2*L-4; l+=2)
3342 {
3343 // increase hopping range
3344 res(l, l+4) = res(l, l+2); // d=2 -> d=4
3345 res(l, l+5) = res(l, l+3); // d=3 -> d=5
3346 res(l+1,l+4) = res(l+1,l+2); // d=1 -> d=3
3347 res(l+1,l+5) = res(l+1,l+3); // d=2 -> d=4
3348 // delete NN hopping except for edges
3349 if (l>=2)
3350 {
3351 res(l, l+2) = 0.;
3352 res(l, l+3) = 0.;
3353 res(l+1,l+2) = 0.;
3354 res(l+1,l+3) = 0.;
3355 }
3356 }
3357// res(0,2*L-2) = res(0,2);
3358// res(0,2*L-1) = res(0,3);
3359// res(1,2*L-2) = res(1,2);
3360// res(1,2*L-1) = res(1,3);
3361 }
3362// cout << "PBC:" << endl;
3363// cout << endl << res << endl << endl;
3364
3365 res.matrix() += res.matrix().adjoint().eval();
3366
3367 return res;
3368}
3369
3370template<typename Scalar>
3371Array<Scalar,Dynamic,Dynamic> hopping_PAM_T (int L, Scalar tfc, Scalar tcc, Scalar tff, Scalar tx, Scalar ty, bool ANCILLA_HOPPING=false, double bugfix=1e-7)
3372{
3373 Array<Scalar,Dynamic,Dynamic> res_tmp = hopping_PAM(L/2,tfc,tcc,tff,tx,ty);
3374 Array<Scalar,Dynamic,Dynamic> res(2*L,2*L); res = 0;
3375
3376 for (int i=0; i<L; ++i)
3377 for (int j=0; j<L; ++j)
3378 {
3379 res(2*i,2*j) = res_tmp(i,j);
3380 if (ANCILLA_HOPPING)
3381 {
3382 res(2*i+1,2*j+1) = res_tmp(i,j);
3383 }
3384 }
3385
3386 for (int i=0; i<2*L-1; ++i)
3387 {
3388 res(i,i+1) += bugfix;
3389 res(i+1,i) += bugfix;
3390 }
3391
3392// cout << res.real() << endl;
3393// cout << endl;
3394// cout << res.imag() << endl;
3395 return res;
3396}
3397
3398ArrayXXd hopping_spinChain (int L, double JA, double JB, double JpA, double JpB, bool ANCILLA_HOPPING=false, bool PBC=false)
3399{
3400 ArrayXXd res_tmp(L,L);
3401
3402 if (PBC)
3403 {
3404 lout << termcolor::yellow << "Warning: ignoring JB, JpB for PBC!" << termcolor::reset << endl;
3405 res_tmp = create_1D_PBC(L,JA,JpA);
3406 }
3407 else
3408 {
3409 res_tmp.setZero();
3410 res_tmp.matrix().diagonal<1> ()(Eigen::seq(0,Eigen::last,2)).setConstant(JA);
3411 res_tmp.matrix().diagonal<1> ()(Eigen::seq(1,Eigen::last,2)).setConstant(JB);
3412 res_tmp.matrix().diagonal<2> ()(Eigen::seq(0,Eigen::last,2)).setConstant(JpA);
3413 res_tmp.matrix().diagonal<2> ()(Eigen::seq(1,Eigen::last,2)).setConstant(JpB);
3414 res_tmp += res_tmp.transpose().eval();
3415 }
3416 return res_tmp;
3417}
3418
3419ArrayXXd hopping_spinChain_T (int L, double JA, double JB, double JpA, double JpB, bool ANCILLA_HOPPING=false, double bugfix=1e-7, bool PBC=false)
3420{
3421 ArrayXXd res_tmp;
3422 if (PBC)
3423 {
3424 lout << termcolor::yellow << "Warning: ignoring JB, JpB for PBC!" << termcolor::reset << endl;
3425 res_tmp = create_1D_PBC(L,JA,JpA);
3426 }
3427 else
3428 {
3429 res_tmp = hopping_spinChain(L,JA,JB,JpA,JpB,false);
3430 }
3431
3432 ArrayXXd res(2*L,2*L); res = 0;
3433
3434 for (int i=0; i<L; ++i)
3435 for (int j=0; j<L; ++j)
3436 {
3437 res(2*i,2*j) = res_tmp(i,j);
3438 if (ANCILLA_HOPPING)
3439 {
3440 res(2*i+1,2*j+1) = res_tmp(i,j);
3441 }
3442 }
3443
3444 for (int i=0; i<2*L-1; ++i)
3445 {
3446 res(i,i+1) += bugfix;
3447 res(i+1,i) += bugfix;
3448 }
3449
3450 return res;
3451}
3452
3453ArrayXXd hopping_ladder (int L, double tPara=1., double tPerp=1., double tPrime=0., double tPPrime=0., bool PBC=false, bool BABA=false)
3454{
3455 ArrayXXd res(L,L);
3456 res.setZero();
3457 if (!PBC)
3458 {
3459 if (!BABA)
3460 {
3461 for (int l=0; l<L; ++l)
3462 {
3463 if (l+1<L) res(l,l+1) = (l%2==0)? tPerp:0.;
3464 if (l+2<L) res(l,l+2) = tPara;
3465
3466 if (l%2==0 and l+3<L) res(l,l+3) = tPrime;
3467 if (l%2==1 and l+1<L) res(l,l+1) = tPrime;
3468
3469 if (l+4<L) res(l,l+4) = tPPrime;
3470 }
3471 }
3472 else
3473 {
3474 for (int l=0; l<L; ++l)
3475 {
3476 if (l+1<L) res(l,l+1) = (l%2==1)? tPerp:0.;
3477 if (l+2<L) res(l,l+2) = tPara;
3478
3479 if (l%2==0 and l+1<L) res(l,l+1) = tPrime;
3480 if (l%2==1 and l+3<L) res(l,l+3) = tPrime;
3481
3482// // different enumeration, doesn't seem to be correct:
3483// if (l+1<L) res(l,l+1) = (l%4==0 or l%4==2)? tPara:0.;
3484// if (l+3<L) res(l,l+3) = (l%4==1 or l%4==3)? tPara:0.;
3485// if (l+1<L) res(l,l+1) = (l%4==1 or l%4==3)? tPerp:0.;
3486// if (l+2<L) res(l,l+2) = tPrime;
3487 }
3488 }
3489 }
3490 else
3491 {
3492 for (int l=0; l<L; ++l)
3493 {
3494 res(l,(l+1)%L) = (l%2==0)? tPerp:0.;
3495 res(l,(l+2)%L) = tPara;
3496
3497 if (l%2==0) res(l,(l+3)%L) = tPrime;
3498 if (l%2==1) res(l,(l+1)%L) = tPrime;
3499
3500 res(l,(l+4)%L) = tPPrime;
3501 }
3502 }
3503 res += res.transpose().eval();
3504 return res;
3505}
3506
3507tuple<double,double,double> params_bilineraBiquadratic_beta (double beta, double J_input=1.)
3508{
3509 double J = J_input+0.5*beta;
3510 double R = -0.5*beta;
3511 double offset = -4./3.*beta;
3512 return make_tuple(J,R,offset);
3513}
3514
3515// returns J, R and offset for bilinear-biquadratic Hamiltonian with J1=1 and J2=beta
3516// H = J*S_i*S_{i+1} - beta*(S_i*S_{i+1})^2
3517// = (J+beta/2)*S_i*S_{i+1} - beta/2 Q_i*Q_{i+1}
3518tuple<double,double,double> params_bilineraBiquadratic_beta (boost::rational<int> beta_rational = boost::rational<int>(-1,3), double J_input=1.)
3519{
3520// double beta = boost::rational_cast<double>(beta_rational);
3521// double J = J_input+0.5*beta;
3522// double R = -0.5*beta;
3523// double offset = -4./3.*beta;
3524// return make_tuple(J,R,offset);
3525 return params_bilineraBiquadratic_beta(boost::rational_cast<double>(beta_rational), J_input);
3526}
3527
3528// H = cos(theta)*S_i*S_{i+1} + sin(theta)*(S_i*S_{i+1})^2
3529tuple<double,double,double> params_bilineraBiquadratic_theta (double theta)
3530{
3531 double beta = -sin(theta);
3532 double J = cos(theta)+0.5*beta;
3533 double R = -0.5*beta;
3534 double offset = -4./3.*beta;
3535 return make_tuple(J,R,offset);
3536}
3537
3538#endif
setZero
void setZero(PivotVector< Symmetry, Scalar_ > &V)
Definition: DmrgPivotVector.h:399
ODD
@ ODD
Definition: DmrgTypedefs.h:128
BC
BC
Definition: DmrgTypedefs.h:161
SUB_LATTICE
SUB_LATTICE
Definition: DmrgTypedefs.h:130
B
@ B
Definition: DmrgTypedefs.h:130
A
@ A
Definition: DmrgTypedefs.h:130
flip_sublattice
SUB_LATTICE flip_sublattice(SUB_LATTICE sublat)
Definition: DmrgTypedefs.h:141
split_kagomeYC_AB
void split_kagomeYC_AB(int Lx, int Ly, std::vector< std::vector< int > > &A, std::vector< std::vector< int > > &B)
Definition: ParamCollection.h:170
hopping_square
ArrayXXd hopping_square(int Lx, int Ly, bool PBCx=false, bool PBCy=true, double lambda=1.)
Definition: ParamCollection.h:134
hopping_PAM_T
Array< Scalar, Dynamic, Dynamic > hopping_PAM_T(int L, Scalar tfc, Scalar tcc, Scalar tff, Scalar tx, Scalar ty, bool ANCILLA_HOPPING=false, double bugfix=1e-7)
Definition: ParamCollection.h:3371
hopping_square_plaquette
ArrayXXd hopping_square_plaquette(int L, int VARIANT=0, double lambda=1.)
Definition: ParamCollection.h:2937
hopping_corannulene
ArrayXXd hopping_corannulene(int L, int VARIANT=0, double lambda=1.)
Definition: ParamCollection.h:2907
create_1D_PBC_AB
ArrayXXd create_1D_PBC_AB(size_t L, double lambda1A=1., double lambda1B=1., double lambda2A=0., double lambda2B=0., bool COMPRESSED=true)
Definition: ParamCollection.h:104
hopping_triangular
ArrayXXd hopping_triangular(int L, int VARIANT=0, double lambda1=1.)
Definition: ParamCollection.h:2183
hopping_coronene
ArrayXXd hopping_coronene(int L, int VARIANT=0, double lambda=1.)
Definition: ParamCollection.h:2858
hopping_spinChain_T
ArrayXXd hopping_spinChain_T(int L, double JA, double JB, double JpA, double JpB, bool ANCILLA_HOPPING=false, double bugfix=1e-7, bool PBC=false)
Definition: ParamCollection.h:3419
create_1D_OBC
ArrayXXd create_1D_OBC(size_t L, double lambda1=1., double lambda2=0.)
Definition: ParamCollection.h:15
hopping_triangulene_dimer
pair< ArrayXXd, vector< SUB_LATTICE > > hopping_triangulene_dimer(int VARIANT=0, double lambda=1., double lambda2=1., string BC="")
Definition: ParamCollection.h:2809
add_triangle
void add_triangle(int i, int j, int k, ArrayXXd &target, double lambda=1.)
Definition: ParamCollection.h:723
hopping_kagome36d
ArrayXXd hopping_kagome36d(double lambda=1.)
Definition: ParamCollection.h:732
hopping_kagomeXC
ArrayXXd hopping_kagomeXC(int Lx, int Ly, bool PBCx=false, bool PBCy=true, double lambda=1., bool VERBOSE=false)
Definition: ParamCollection.h:552
hopping_spinChain
ArrayXXd hopping_spinChain(int L, double JA, double JB, double JpA, double JpB, bool ANCILLA_HOPPING=false, bool PBC=false)
Definition: ParamCollection.h:3398
add_edge
void add_edge(int i, int j, vector< pair< size_t, size_t > > &target)
Definition: ParamCollection.h:2248
create_1D_PBC
ArrayXXd create_1D_PBC(size_t L, double lambda1=1., double lambda2=0., bool COMPRESSED=false)
Definition: ParamCollection.h:30
params_bilineraBiquadratic_beta
tuple< double, double, double > params_bilineraBiquadratic_beta(double beta, double J_input=1.)
Definition: ParamCollection.h:3507
hopping_triangulene
pair< ArrayXXd, vector< SUB_LATTICE > > hopping_triangulene(int L, int VARIANT=0, double lambda=1., double lambda2=1., string BC="")
Definition: ParamCollection.h:2596
push_back_KondoUnpacked
void push_back_KondoUnpacked(vector< Param > &params, size_t L, double J, double t, size_t D, bool START_WITH_SPIN=true)
Definition: ParamCollection.h:3086
add_tetrahedron
void add_tetrahedron(int i, int j, int k, int l, vector< pair< size_t, size_t > > &target)
Definition: ParamCollection.h:2230
hopping_PPV
pair< ArrayXXd, vector< SUB_LATTICE > > hopping_PPV(int L, int VARIANT=0, double t0=1., double tsingle=1., double tdouble=1., string BC="")
Definition: ParamCollection.h:2528
hopping_fullerene
ArrayXXd hopping_fullerene(int L=60, int VARIANT=0, double lambda1=1., double lambda2=1.)
Definition: ParamCollection.h:1540
hopping_ladder
ArrayXXd hopping_ladder(int L, double tPara=1., double tPerp=1., double tPrime=0., double tPPrime=0., bool PBC=false, bool BABA=false)
Definition: ParamCollection.h:3453
hopping_MnRing
ArrayXXi hopping_MnRing(int Ncells, double J1, double J2, double J3, double J4=0., double J5=0., double J6=0., double J7=0.)
Definition: ParamCollection.h:3043
Tinf_params_spins
vector< Param > Tinf_params_spins(size_t L, size_t Ly, size_t maxPower=1ul, size_t DA=2ul, size_t DB=2ul, bool SOFT=false)
Definition: ParamCollection.h:3249
params_bilineraBiquadratic_theta
tuple< double, double, double > params_bilineraBiquadratic_theta(double theta)
Definition: ParamCollection.h:3529
hopping_sodaliteCage
ArrayXXd hopping_sodaliteCage(int L=60, int VARIANT=0, double lambda1=1.)
Definition: ParamCollection.h:2253
hopping_fullerene_C40Td
ArrayXXd hopping_fullerene_C40Td(int VARIANT=0, double lambda1=1., double lambda2=1.)
Definition: ParamCollection.h:1462
hopping_kagomeYC_BAB
ArrayXXd hopping_kagomeYC_BAB(int Lx, int Ly, bool PBCx=false, bool PBCy=true, double lambda=1., bool VERBOSE=false)
Definition: ParamCollection.h:381
split_kagomeYC_BAB
void split_kagomeYC_BAB(int Lx, int Ly, std::vector< std::vector< int > > &A, std::vector< std::vector< int > > &B)
Definition: ParamCollection.h:207
Tinf_params_fermions
vector< Param > Tinf_params_fermions(size_t L, size_t Ly, size_t maxPower=1ul)
Definition: ParamCollection.h:3219
extend_to_thermal
ArrayXXd extend_to_thermal(const ArrayXXd &tFull, double factor)
Definition: ParamCollection.h:65
hopping_PAM
Array< Scalar, Dynamic, Dynamic > hopping_PAM(int L, Scalar tfc, Scalar tcc, Scalar tff, Scalar tx, Scalar ty, bool PBC=false)
Definition: ParamCollection.h:3316
find_x_kagomeYC
int find_x_kagomeYC(int index, int L, const std::vector< std::vector< int > > &A, const std::vector< std::vector< int > > &B)
Definition: ParamCollection.h:245
calc_distanceMatrix
ArrayXXi calc_distanceMatrix(ArrayXXd adjacencyMatrix)
Definition: ParamCollection.h:3000
hopping_kagomeYC_AB
ArrayXXd hopping_kagomeYC_AB(int Lx, int Ly, bool PBCx=false, bool PBCy=true, double lambda=1.)
Definition: ParamCollection.h:287
create_1D_AB
ArrayXXd create_1D_AB(size_t L, double lambda1A=1., double lambda1B=1., double lambda2A=0., double lambda2B=0.)
Definition: ParamCollection.h:82
hopping_Platonic
ArrayXXd hopping_Platonic(int L, int VARIANT=0, double lambda1=1.)
Definition: ParamCollection.h:2013
conjIfImag
double conjIfImag(double x)
Definition: ParamCollection.h:3312
split_kagomeXC
void split_kagomeXC(int Lx, int Ly, std::vector< std::vector< int > > &A, std::vector< std::vector< int > > &B, std::vector< std::vector< int > > &C, std::vector< std::vector< int > > &D)
Definition: ParamCollection.h:476
hopping_triangularYC
ArrayXXd hopping_triangularYC(int Lx, int Ly, bool PBCx=false, bool PBCy=true, double lambda=1.)
Definition: ParamCollection.h:695
hopping_Archimedean
ArrayXXd hopping_Archimedean(string vertex_conf, int VARIANT=0, double lambda1=1., double lambda2=1.)
Definition: ParamCollection.h:1006
hexagonalFlake
ArrayXXd hexagonalFlake(int L, double lambda=1.)
Definition: ParamCollection.h:911
triangularFlake
ArrayXXd triangularFlake(int L, double lambda=1.)
Definition: ParamCollection.h:770
hopping_Mn32
ArrayXXd hopping_Mn32(double lambda_cap=1., double lambda_corner=0., double lambda_edge=1., int VARIANT=0)
Definition: ParamCollection.h:2456