HΦ  3.2.0
diagonalcalc.c
Go to the documentation of this file.
1 /* HPhi - Quantum Lattice Model Simulator */
2 /* Copyright (C) 2015 The University of Tokyo */
3 
4 /* This program is free software: you can redistribute it and/or modify */
5 /* it under the terms of the GNU General Public License as published by */
6 /* the Free Software Foundation, either version 3 of the License, or */
7 /* (at your option) any later version. */
8 
9 /* This program is distributed in the hope that it will be useful, */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
12 /* GNU General Public License for more details. */
13 
14 /* You should have received a copy of the GNU General Public License */
15 /* along with this program. If not, see <http://www.gnu.org/licenses/>. */
16 
36 #include <bitcalc.h>
37 #include "FileIO.h"
38 #include "diagonalcalc.h"
39 #include "mltplySpinCore.h"
40 #include "wrapperMPI.h"
41 
42 
43 int SetDiagonalTETransfer(
44  long unsigned int isite1,
45  double dtmp_V,
46  long unsigned int spin,
47  struct BindStruct *X,
48  double complex *tmp_v0,
49  double complex *tmp_v1
50 );
51 
52 int SetDiagonalTEInterAll(
53  long unsigned int isite1,
54  long unsigned int isite2,
55  long unsigned int isigma1,
56  long unsigned int isigma2,
57  double dtmp_V,
58  struct BindStruct *X,
59  double complex *tmp_v0,
60  double complex *tmp_v1
61 );
62 
63 int SetDiagonalTEChemi(
64  long unsigned int isite1,
65  long unsigned int spin,
66  double dtmp_V,
67  struct BindStruct *X,
68  double complex *tmp_v0,
69  double complex *tmp_v1
70 );
71 
84 int diagonalcalc
85 (
86  struct BindStruct *X
87  ){
88 
89  FILE *fp;
90  long unsigned int i,j;
91  long unsigned int isite1,isite2;
92  long unsigned int spin;
93  double tmp_V;
94 
95  /*[s] For InterAll*/
96  long unsigned int A_spin,B_spin;
97  /*[e] For InterAll*/
98  long unsigned int i_max=X->Check.idim_max;
99 
100  fprintf(stdoutMPI, "%s", cProStartCalcDiag);
101  TimeKeeper(X, cFileNameTimeKeep, cDiagonalCalcStart, "a");
102 
103 #pragma omp parallel for default(none) private(j) shared(list_Diagonal) firstprivate(i_max)
104  for(j = 1;j <= i_max; j++){
105  list_Diagonal[j]=0.0;
106  }
107 
108  if(X->Def.NCoulombIntra>0){
109  if(childfopenMPI(cFileNameCheckCoulombIntra, "w", &fp)!=0){
110  return -1;
111  }
112  for(i = 0; i < X->Def.NCoulombIntra; i++){
113  isite1 = X->Def.CoulombIntra[i][0]+1;
114  tmp_V = X->Def.ParaCoulombIntra[i];
115  fprintf(fp,"i=%ld isite1=%ld tmp_V=%lf \n",i,isite1,tmp_V);
116  SetDiagonalCoulombIntra(isite1, tmp_V, X);
117  }
118  fclose(fp);
119  }
120 
121  if(X->Def.EDNChemi>0){
122  if(childfopenMPI(cFileNameCheckChemi,"w", &fp)!=0){
123  return -1;
124  }
125  for(i = 0; i < X->Def.EDNChemi; i++){
126  isite1 = X->Def.EDChemi[i]+1;
127  spin = X->Def.EDSpinChemi[i];
128  tmp_V = -X->Def.EDParaChemi[i];
129  fprintf(fp,"i=%ld spin=%ld isite1=%ld tmp_V=%lf \n",i,spin,isite1,tmp_V);
130  if(SetDiagonalChemi(isite1, tmp_V,spin, X) !=0){
131  return -1;
132  }
133  }
134  fclose(fp);
135  }
136 
137  if(X->Def.NCoulombInter>0){
138  if(childfopenMPI(cFileNameCheckInterU,"w", &fp)!=0){
139  return -1;
140  }
141  for(i = 0; i < X->Def.NCoulombInter; i++){
142  isite1 = X->Def.CoulombInter[i][0]+1;
143  isite2 = X->Def.CoulombInter[i][1]+1;
144  tmp_V = X->Def.ParaCoulombInter[i];
145  fprintf(fp,"i=%ld isite1=%ld isite2=%ld tmp_V=%lf \n",i,isite1,isite2,tmp_V);
146  if(SetDiagonalCoulombInter(isite1, isite2, tmp_V, X) !=0){
147  return -1;
148  }
149  }
150  fclose(fp);
151  }
152  if(X->Def.NHundCoupling>0){
153  if(childfopenMPI(cFileNameCheckHund,"w", &fp) !=0){
154  return -1;
155  }
156  for(i = 0; i < X->Def.NHundCoupling; i++){
157  isite1 = X->Def.HundCoupling[i][0]+1;
158  isite2 = X->Def.HundCoupling[i][1]+1;
159  tmp_V = -X->Def.ParaHundCoupling[i];
160  if(SetDiagonalHund(isite1, isite2, tmp_V, X) !=0){
161  return -1;
162  }
163  fprintf(fp,"i=%ld isite1=%ld isite2=%ld tmp_V=%lf \n",i,isite1,isite2,tmp_V);
164  }
165  fclose(fp);
166  }
167 
168  if(X->Def.NInterAll_Diagonal>0){
169  if(childfopenMPI(cFileNameCheckInterAll,"w", &fp) !=0){
170  return -1;
171  }
172  for(i = 0; i < X->Def.NInterAll_Diagonal; i++){
173  isite1=X->Def.InterAll_Diagonal[i][0]+1;
174  A_spin=X->Def.InterAll_Diagonal[i][1];
175  isite2=X->Def.InterAll_Diagonal[i][2]+1;
176  B_spin=X->Def.InterAll_Diagonal[i][3];
177  tmp_V = X->Def.ParaInterAll_Diagonal[i];
178  fprintf(fp,"i=%ld isite1=%ld A_spin=%ld isite2=%ld B_spin=%ld tmp_V=%lf \n", i, isite1, A_spin, isite2, B_spin, tmp_V);
179  SetDiagonalInterAll(isite1, isite2, A_spin, B_spin, tmp_V, X);
180  }
181  fclose(fp);
182  }
183 
184  TimeKeeper(X, cFileNameTimeKeep, cDiagonalCalcFinish, "a");
185  fprintf(stdoutMPI, "%s", cProEndCalcDiag);
186  return 0;
187 }
188 
196 int diagonalcalcForTE
197  (
198  const int _istep,
199  struct BindStruct *X,
200  double complex *tmp_v0,
201  double complex *tmp_v1
202  ) {
203 
204  long unsigned int i;
205  long unsigned int isite1, isite2;
206  long unsigned int A_spin, B_spin;
207  double tmp_V;
208 
209  if (X->Def.NTETransferDiagonal[_istep] > 0) {
210  for (i = 0; i < X->Def.NTETransferDiagonal[_istep]; i++) {
211  isite1 = X->Def.TETransferDiagonal[_istep][i][0] + 1;
212  A_spin = X->Def.TETransferDiagonal[_istep][i][1];
213  tmp_V = X->Def.ParaTETransferDiagonal[_istep][i];
214  SetDiagonalTETransfer(isite1, tmp_V, A_spin, X, tmp_v0, tmp_v1);
215  }
216  }
217  else if (X->Def.NTEInterAllDiagonal[_istep] >0) {
218  for (i = 0; i < X->Def.NTEInterAllDiagonal[_istep]; i++) {
219  //Assume n_{1\sigma_1} n_{2\sigma_2}
220  isite1 = X->Def.TEInterAllDiagonal[_istep][i][0] + 1;
221  A_spin = X->Def.TEInterAllDiagonal[_istep][i][1];
222  isite2 = X->Def.TEInterAllDiagonal[_istep][i][2] + 1;
223  B_spin = X->Def.TEInterAllDiagonal[_istep][i][3];
224  tmp_V = X->Def.ParaTEInterAllDiagonal[_istep][i];
225 
226  if (SetDiagonalTEInterAll(isite1, isite2, A_spin, B_spin, tmp_V, X, tmp_v0, tmp_v1) != 0) {
227  return -1;
228  }
229  }
230 
231  if (X->Def.NTEChemi[_istep] > 0) {
232  for(i=0; i< X->Def.NTEChemi[_istep]; i++) {
233  isite1 = X->Def.TEChemi[_istep][i] + 1;
234  A_spin = X->Def.SpinTEChemi[_istep][i];
235  tmp_V = -X->Def.ParaTEChemi[_istep][i];
236  if (SetDiagonalTEChemi(isite1, A_spin, tmp_V, X, tmp_v0, tmp_v1) != 0) {
237  return -1;
238  }
239  }
240  }
241  }
242  return 0;
243 }
244 
245 
260 int SetDiagonalCoulombIntra
261 (
262  long unsigned int isite1,
263  double dtmp_V,
264  struct BindStruct *X
265  ){
266  long unsigned int is;
267  long unsigned int ibit;
268  long unsigned int is1_up, is1_down;
269 
270  long unsigned int j;
271  long unsigned int i_max=X->Check.idim_max;
272 
273  /*
274  When isite1 is in the inter process region
275  */
276  if (isite1 > X->Def.Nsite){
277 
278  switch (X->Def.iCalcModel) {
279 
280  case HubbardGC:
281  case KondoGC:
282  case Hubbard:
283  case Kondo:
284 
285  is1_up = X->Def.Tpow[2 * isite1 - 2];
286  is1_down = X->Def.Tpow[2 * isite1 - 1];
287  is = is1_up + is1_down;
288  ibit = (unsigned long int)myrank & is;
289  if (ibit == is) {
290 #pragma omp parallel for default(none) shared(list_Diagonal) \
291  firstprivate(i_max, dtmp_V) private(j)
292  for (j = 1; j <= i_max; j++) list_Diagonal[j] += dtmp_V;
293  }
294 
295  break; /*case HubbardGC, KondoGC, Hubbard, Kondo:*/
296 
297  case Spin:
298  case SpinGC:
299  /*
300  They do not have the Coulomb term
301  */
302  break;
303 
304  default:
305  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
306  return -1;
307  //break;
308 
309  }/*switch (X->Def.iCalcModel)*/
310 
311  return 0;
312 
313  }/*if (isite1 >= X->Def.Nsite*/
314  else{
315  switch (X->Def.iCalcModel){
316  case HubbardGC:
317  is1_up = X->Def.Tpow[2*isite1-2];
318  is1_down = X->Def.Tpow[2*isite1-1];
319  is=is1_up+is1_down;
320 #pragma omp parallel for default(none) shared(list_Diagonal, list_1) firstprivate(i_max, is, dtmp_V) private(ibit)
321  for(j = 1;j <= i_max;j++){
322  ibit=(j-1)&is;
323  if(ibit==is){
324  list_Diagonal[j]+=dtmp_V;
325  }
326  }
327 
328  break;
329  case KondoGC:
330  case Hubbard:
331  case Kondo:
332  is1_up = X->Def.Tpow[2*isite1-2];
333  is1_down = X->Def.Tpow[2*isite1-1];
334  is=is1_up+is1_down;
335 #pragma omp parallel for default(none) shared(list_Diagonal, list_1) firstprivate(i_max, is, dtmp_V) private(ibit)
336  for(j = 1;j <= i_max;j++){
337  ibit=list_1[j]&is;
338  if(ibit==is){
339  list_Diagonal[j]+=dtmp_V;
340  }
341  }
342  break;
343 
344  case Spin:
345  case SpinGC:
346  break;
347 
348  default:
349  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
350  return -1;
351  //break;
352  }
353  }
354  return 0;
355 }
356 
357 
373 int SetDiagonalChemi
374 (
375  long unsigned int isite1,
376  double dtmp_V,
377  long unsigned int spin,
378  struct BindStruct *X
379  ){
380  long unsigned int is1_up;
381  long unsigned int ibit1_up;
382  long unsigned int num1;
383  long unsigned int isigma1 =spin;
384  long unsigned int is1,ibit1;
385 
386  long unsigned int j;
387  long unsigned int i_max=X->Check.idim_max;
388 
389  /*
390  When isite1 is in the inter process region
391  */
392  if (isite1 > X->Def.Nsite){
393 
394  switch (X->Def.iCalcModel) {
395 
396  case HubbardGC:
397  case KondoGC:
398  case Hubbard:
399  case Kondo:
400 
401  if (spin == 0) {
402  is1 = X->Def.Tpow[2 * isite1 - 2];
403  }
404  else {
405  is1 = X->Def.Tpow[2 * isite1 - 1];
406  }
407  ibit1 = (unsigned long int)myrank & is1;
408  num1 = ibit1 / is1;
409 #pragma omp parallel for default(none) shared(list_Diagonal) \
410  firstprivate(i_max, dtmp_V, num1) private(j)
411  for (j = 1; j <= i_max; j++) list_Diagonal[j] += num1*dtmp_V;
412 
413  break;/*case HubbardGC, case KondoGC, Hubbard, Kondo:*/
414 
415  case SpinGC:
416  case Spin:
417 
418  if (X->Def.iFlgGeneralSpin == FALSE) {
419  is1_up = X->Def.Tpow[isite1 - 1];
420  ibit1_up = (((unsigned long int)myrank& is1_up) / is1_up) ^ (1 - spin);
421 #pragma omp parallel for default(none) shared(list_Diagonal) \
422 firstprivate(i_max, dtmp_V, ibit1_up) private(j)
423  for (j = 1; j <= i_max; j++) list_Diagonal[j] += dtmp_V * ibit1_up;
424  } /*if (X->Def.iFlgGeneralSpin == FALSE)*/
425  else /*if (X->Def.iFlgGeneralSpin == TRUE)*/ {
426  num1 = BitCheckGeneral((unsigned long int)myrank,
427  isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
428  if (num1 != 0) {
429 #pragma omp parallel for default(none) shared(list_Diagonal) \
430 firstprivate(i_max, dtmp_V) private(j)
431  for (j = 1; j <= i_max; j++) list_Diagonal[j] += dtmp_V;
432  }/*if (num1 != 0)*/
433  }/*if (X->Def.iFlgGeneralSpin == TRUE)*/
434  break;/*case SpinGC, Spin:*/
435 
436  default:
437  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
438  return -1;
439 
440  } /*switch (X->Def.iCalcModel)*/
441 
442  return 0;
443 
444  }/*if (isite1 >= X->Def.Nsite*/
445 
446  switch (X->Def.iCalcModel){
447  case HubbardGC:
448  if(spin==0){
449  is1 = X->Def.Tpow[2*isite1-2];
450  }else{
451  is1 = X->Def.Tpow[2*isite1-1];
452  }
453 #pragma omp parallel for default(none) shared(list_1, list_Diagonal) firstprivate(i_max, dtmp_V, is1) private(num1, ibit1)
454  for(j = 1;j <= i_max;j++){
455 
456  ibit1 = (j-1)&is1;
457  num1 = ibit1/is1;
458  //fprintf(stdoutMPI, "DEBUG: spin=%ld is1=%ld: isite1=%ld j=%ld num1=%ld \n",spin,is1,isite1,j,num1);
459 
460  list_Diagonal[j]+=num1*dtmp_V;
461  }
462  break;
463  case KondoGC:
464  case Hubbard:
465  case Kondo:
466  if(spin==0){
467  is1 = X->Def.Tpow[2*isite1-2];
468  }else{
469  is1 = X->Def.Tpow[2*isite1-1];
470  }
471 
472 #pragma omp parallel for default(none) shared(list_1, list_Diagonal) firstprivate(i_max, dtmp_V, is1) private(num1, ibit1)
473  for(j = 1;j <= i_max;j++){
474 
475  ibit1 = list_1[j]&is1;
476  num1 = ibit1/is1;
477  list_Diagonal[j]+=num1*dtmp_V;
478  }
479  break;
480 
481  case SpinGC:
482  if(X->Def.iFlgGeneralSpin==FALSE){
483  is1_up = X->Def.Tpow[isite1-1];
484 #pragma omp parallel for default(none) shared(list_1, list_Diagonal) firstprivate(i_max, dtmp_V, is1_up, spin) private(num1, ibit1_up)
485  for(j = 1;j <= i_max;j++){
486  ibit1_up=(((j-1)& is1_up)/is1_up)^(1-spin);
487  list_Diagonal[j] += dtmp_V * ibit1_up;
488  }
489  }
490  else{
491 #pragma omp parallel for default(none) shared(list_Diagonal) firstprivate(i_max, dtmp_V, isite1, isigma1, X) private(j, num1)
492  for(j = 1;j <= i_max; j++){
493  num1=BitCheckGeneral (j-1, isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
494  if(num1 != 0){
495  list_Diagonal[j] += dtmp_V;
496  }
497  }
498  }
499  break;
500 
501  case Spin:
502  if(X->Def.iFlgGeneralSpin==FALSE){
503  is1_up = X->Def.Tpow[isite1-1];
504 #pragma omp parallel for default(none) shared(list_1, list_Diagonal) firstprivate(i_max, dtmp_V, is1_up, spin) private(num1, ibit1_up)
505  for(j = 1;j <= i_max;j++){
506  ibit1_up=((list_1[j]& is1_up)/is1_up)^(1-spin);
507  list_Diagonal[j] += dtmp_V * ibit1_up;
508  }
509  }
510  else{
511 #pragma omp parallel for default(none) shared(list_Diagonal, list_1) firstprivate(i_max, dtmp_V, isite1, isigma1, X) private(j, num1)
512  for(j = 1;j <= i_max; j++){
513  num1=BitCheckGeneral (list_1[j], isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
514  if(num1 != 0){
515  list_Diagonal[j] += dtmp_V;
516  }
517  }
518  }
519 
520  break;
521  default:
522  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
523  return -1;
524  }
525 
526  return 0;
527 }
528 
543 int SetDiagonalCoulombInter
544 (
545  long unsigned int isite1,
546  long unsigned int isite2,
547  double dtmp_V,
548  struct BindStruct *X
549  )
550 {
551 
552  long unsigned int is1_up, is1_down;
553  long unsigned int ibit1_up, ibit1_down;
554  long unsigned int num1;
555  long unsigned int is2_up, is2_down;
556  long unsigned int ibit2_up, ibit2_down;
557  long unsigned int num2;
558 
559  long unsigned int j;
560  long unsigned int i_max=X->Check.idim_max;
561 
562  /*
563  Force isite1 <= isite2
564  */
565  if (isite2 < isite1) {
566  j = isite2;
567  isite2 = isite1;
568  isite1 = j;
569  }/*if (isite2 < isite1)*/
570  /*
571  When isite1 & site2 are in the inter process region
572  */
573  if (/*isite2 => */ isite1 > X->Def.Nsite) {
574 
575  switch (X->Def.iCalcModel) {
576 
577  case HubbardGC:
578  case KondoGC:
579  case Hubbard:
580  case Kondo:
581 
582  is1_up = X->Def.Tpow[2 * isite1 - 2];
583  is1_down = X->Def.Tpow[2 * isite1 - 1];
584  is2_up = X->Def.Tpow[2 * isite2 - 2];
585  is2_down = X->Def.Tpow[2 * isite2 - 1];
586 
587  num1 = 0;
588  num2 = 0;
589 
590  ibit1_up = (unsigned long int)myrank&is1_up;
591  num1 += ibit1_up / is1_up;
592  ibit1_down = (unsigned long int)myrank&is1_down;
593  num1 += ibit1_down / is1_down;
594 
595  ibit2_up = (unsigned long int)myrank&is2_up;
596  num2 += ibit2_up / is2_up;
597  ibit2_down = (unsigned long int)myrank&is2_down;
598  num2 += ibit2_down / is2_down;
599 
600 #pragma omp parallel for default(none) shared(list_Diagonal) \
601  firstprivate(i_max, dtmp_V, num1, num2) private(j)
602  for (j = 1; j <= i_max; j++) list_Diagonal[j] += num1*num2*dtmp_V;
603 
604  break;/*case HubbardGC, KondoGC, Hubbard, Kondo:*/
605 
606  case Spin:
607  case SpinGC:
608 #pragma omp parallel for default(none) shared(list_Diagonal) firstprivate(i_max, dtmp_V)
609  for (j = 1; j <= i_max; j++) {
610  list_Diagonal[j] += dtmp_V;
611  }
612  break;/*case Spin, SpinGC*/
613 
614  default:
615  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
616  return -1;
617 
618  }/*switch (X->Def.iCalcModel)*/
619 
620  return 0;
621 
622  }/*if (isite1 > X->Def.Nsite)*/
623  else if (isite2 > X->Def.Nsite /* => isite1 */) {
624 
625  switch(X->Def.iCalcModel){
626  case HubbardGC:
627  case KondoGC:
628  case Hubbard:
629  case Kondo:
630  is1_up = X->Def.Tpow[2 * isite1 - 2];
631  is1_down = X->Def.Tpow[2 * isite1 - 1];
632  is2_up = X->Def.Tpow[2 * isite2 - 2];
633  is2_down = X->Def.Tpow[2 * isite2 - 1];
634  num2 = 0;
635  ibit2_up = (unsigned long int)myrank&is2_up;
636  num2 += ibit2_up / is2_up;
637  ibit2_down = (unsigned long int)myrank&is2_down;
638  num2 += ibit2_down / is2_down;
639  break;
640 
641  case Spin:
642  case SpinGC:
643  break;
644 
645  default:
646  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
647  return -1;
648  }
649 
650  switch (X->Def.iCalcModel) {
651 
652  case HubbardGC:
653 
654 #pragma omp parallel for default(none) shared(list_Diagonal) \
655 firstprivate(i_max, dtmp_V, num2, is1_up, is1_down) \
656 private(num1, ibit1_up, ibit1_down, j)
657  for (j = 1; j <= i_max; j++) {
658  num1 = 0;
659  ibit1_up = (j - 1)&is1_up;
660  num1 += ibit1_up / is1_up;
661  ibit1_down = (j - 1)&is1_down;
662  num1 += ibit1_down / is1_down;
663 
664  list_Diagonal[j] += num1*num2*dtmp_V;
665  }
666 
667  break;/*case HubbardGC*/
668 
669  case KondoGC:
670  case Hubbard:
671  case Kondo:
672 
673 #pragma omp parallel for default(none) shared(list_1, list_Diagonal) \
674 firstprivate(i_max, dtmp_V, is1_up, is1_down, num2) \
675 private(num1, ibit1_up, ibit1_down, j)
676  for (j = 1; j <= i_max; j++) {
677  num1 = 0;
678  ibit1_up = list_1[j] & is1_up;
679  num1 += ibit1_up / is1_up;
680  ibit1_down = list_1[j] & is1_down;
681  num1 += ibit1_down / is1_down;
682 
683  list_Diagonal[j] += num1*num2*dtmp_V;
684  }
685  break;/*case KondoGC, Hubbard, Kondo:*/
686 
687  case Spin:
688  case SpinGC:
689 #pragma omp parallel for default(none) shared(list_Diagonal) firstprivate(i_max, dtmp_V)
690  for (j = 1; j <= i_max; j++) {
691  list_Diagonal[j] += dtmp_V;
692  }
693  break;/* case Spin, SpinGC:*/
694 
695  default:
696  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
697  return -1;
698 
699  }/*switch (X->Def.iCalcModel)*/
700 
701  return 0;
702 
703  }/*else if (isite2 > X->Def.Nsite)*/
704  else{
705  switch (X->Def.iCalcModel){
706  case HubbardGC: //list_1[j] -> j-1
707  is1_up = X->Def.Tpow[2*isite1-2];
708  is1_down = X->Def.Tpow[2*isite1-1];
709  is2_up = X->Def.Tpow[2*isite2-2];
710  is2_down = X->Def.Tpow[2*isite2-1];
711 #pragma omp parallel for default(none) shared( list_Diagonal) firstprivate(i_max, dtmp_V, is1_up, is1_down, is2_up, is2_down) private(num1, ibit1_up, ibit1_down, num2, ibit2_up, ibit2_down)
712  for(j = 1;j <= i_max;j++){
713  num1=0;
714  num2=0;
715  ibit1_up=(j-1)&is1_up;
716  num1+=ibit1_up/is1_up;
717  ibit1_down=(j-1)&is1_down;
718  num1+=ibit1_down/is1_down;
719 
720  ibit2_up=(j-1)&is2_up;
721  num2+=ibit2_up/is2_up;
722  ibit2_down=(j-1)&is2_down;
723  num2+=ibit2_down/is2_down;
724 
725  list_Diagonal[j]+=num1*num2*dtmp_V;
726  }
727  break;
728  case KondoGC:
729  case Hubbard:
730  case Kondo:
731  is1_up = X->Def.Tpow[2*isite1-2];
732  is1_down = X->Def.Tpow[2*isite1-1];
733  is2_up = X->Def.Tpow[2*isite2-2];
734  is2_down = X->Def.Tpow[2*isite2-1];
735 
736 #pragma omp parallel for default(none) shared(list_1, list_Diagonal) firstprivate(i_max, dtmp_V, is1_up, is1_down, is2_up, is2_down) private(num1, ibit1_up, ibit1_down, num2, ibit2_up, ibit2_down)
737  for(j = 1;j <= i_max;j++){
738  num1=0;
739  num2=0;
740  ibit1_up=list_1[j]&is1_up;
741  num1+=ibit1_up/is1_up;
742  ibit1_down=list_1[j]&is1_down;
743  num1+=ibit1_down/is1_down;
744 
745  ibit2_up=list_1[j]&is2_up;
746  num2+=ibit2_up/is2_up;
747  ibit2_down=list_1[j]&is2_down;
748  num2+=ibit2_down/is2_down;
749 
750  list_Diagonal[j]+=num1*num2*dtmp_V;
751  }
752  break;
753 
754  case Spin:
755  case SpinGC:
756 #pragma omp parallel for default(none) shared(list_Diagonal) firstprivate(i_max, dtmp_V)
757  for(j = 1;j <= i_max; j++){
758  list_Diagonal[j] += dtmp_V;
759  }
760  break;
761  default:
762  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
763  return -1;
764  }
765  }
766 
767  return 0;
768 }
769 
784 int SetDiagonalHund
785 (
786  long unsigned int isite1,
787  long unsigned int isite2,
788  double dtmp_V,
789  struct BindStruct *X
790  ){
791 
792  long unsigned int is1_up, is1_down;
793  long unsigned int ibit1_up, ibit1_down;
794  long unsigned int num1_up, num1_down;
795  long unsigned int is2_up, is2_down;
796  long unsigned int ibit2_up, ibit2_down;
797  long unsigned int num2_up, num2_down;
798 
799  long unsigned int is_up;
800  long unsigned int ibit;
801  long unsigned int j;
802  long unsigned int i_max=X->Check.idim_max;
803  /*
804  Force isite1 <= isite2
805  */
806  if (isite2 < isite1) {
807  j = isite2;
808  isite2 = isite1;
809  isite1 = j;
810  }
811  /*
812  When isite1 & site2 are in the inter process region
813  */
814  if (/*isite2 >= */ isite1 > X->Def.Nsite){
815 
816  switch (X->Def.iCalcModel) {
817 
818  case HubbardGC:
819  case KondoGC:
820  case Hubbard:
821  case Kondo:
822 
823  is1_up = X->Def.Tpow[2 * isite1 - 2];
824  is1_down = X->Def.Tpow[2 * isite1 - 1];
825  is2_up = X->Def.Tpow[2 * isite2 - 2];
826  is2_down = X->Def.Tpow[2 * isite2 - 1];
827 
828  num1_up = 0;
829  num1_down = 0;
830  num2_up = 0;
831  num2_down = 0;
832 
833  ibit1_up = (unsigned long int)myrank &is1_up;
834  num1_up = ibit1_up / is1_up;
835  ibit1_down = (unsigned long int)myrank &is1_down;
836  num1_down = ibit1_down / is1_down;
837 
838  ibit2_up = (unsigned long int)myrank &is2_up;
839  num2_up = ibit2_up / is2_up;
840  ibit2_down = (unsigned long int)myrank &is2_down;
841  num2_down = ibit2_down / is2_down;
842 
843 #pragma omp parallel for default(none) shared(list_Diagonal) \
844  firstprivate(i_max, dtmp_V, num1_up, num1_down, num2_up, num2_down) private(j)
845  for (j = 1; j <= i_max; j++)
846  list_Diagonal[j] += dtmp_V*(num1_up*num2_up + num1_down*num2_down);
847 
848  break;/*case HubbardGC, KondoGC, Hubbard, Kondo:*/
849 
850  case SpinGC:
851  case Spin:
852 
853  is1_up = X->Def.Tpow[isite1 - 1];
854  is2_up = X->Def.Tpow[isite2 - 1];
855  is_up = is1_up + is2_up;
856  ibit = (unsigned long int)myrank & is_up;
857  if (ibit == 0 || ibit == is_up) {
858 #pragma omp parallel for default(none) shared(list_Diagonal) \
859 firstprivate(i_max, dtmp_V) private(j)
860  for (j = 1; j <= i_max; j++) list_Diagonal[j] += dtmp_V;
861  }
862  break;/*case SpinGC, Spin:*/
863 
864  default:
865  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
866  return -1;
867  }
868 
869  return 0;
870 
871  }/*if (isite1 > X->Def.Nsite)*/
872  else if (isite2 > X->Def.Nsite /* >= isite1 */) {
873 
874  switch (X->Def.iCalcModel) {
875 
876  case HubbardGC:
877 
878  is1_up = X->Def.Tpow[2 * isite1 - 2];
879  is1_down = X->Def.Tpow[2 * isite1 - 1];
880  is2_up = X->Def.Tpow[2 * isite2 - 2];
881  is2_down = X->Def.Tpow[2 * isite2 - 1];
882 
883  num2_up = 0;
884  num2_down = 0;
885 
886  ibit2_up = (unsigned long int)myrank &is2_up;
887  num2_up = ibit2_up / is2_up;
888  ibit2_down = (unsigned long int)myrank &is2_down;
889  num2_down = ibit2_down / is2_down;
890 
891 #pragma omp parallel for default(none) shared( list_Diagonal) \
892 firstprivate(i_max, dtmp_V, num2_up, num2_down, is1_up, is1_down) \
893 private(num1_up, num1_down, ibit1_up, ibit1_down, j)
894  for (j = 1; j <= i_max; j++) {
895  num1_up = 0;
896  num1_down = 0;
897 
898  ibit1_up = (j - 1)&is1_up;
899  num1_up = ibit1_up / is1_up;
900  ibit1_down = (j - 1)&is1_down;
901  num1_down = ibit1_down / is1_down;
902 
903  list_Diagonal[j] += dtmp_V*(num1_up*num2_up + num1_down*num2_down);
904  }
905  break;/*case HubbardGC:*/
906 
907  case KondoGC:
908  case Hubbard:
909  case Kondo:
910 
911  is1_up = X->Def.Tpow[2 * isite1 - 2];
912  is1_down = X->Def.Tpow[2 * isite1 - 1];
913  is2_up = X->Def.Tpow[2 * isite2 - 2];
914  is2_down = X->Def.Tpow[2 * isite2 - 1];
915 
916  num2_up = 0;
917  num2_down = 0;
918 
919  ibit2_up = (unsigned long int)myrank&is2_up;
920  num2_up = ibit2_up / is2_up;
921  ibit2_down = (unsigned long int)myrank&is2_down;
922  num2_down = ibit2_down / is2_down;
923 
924 #pragma omp parallel for default(none) shared(list_1, list_Diagonal) \
925 firstprivate(i_max, dtmp_V, num2_up, num2_down, is1_up, is1_down) \
926 private(num1_up, num1_down, ibit1_up, ibit1_down, j)
927  for (j = 1; j <= i_max; j++) {
928  num1_up = 0;
929  num1_down = 0;
930 
931  ibit1_up = list_1[j] & is1_up;
932  num1_up = ibit1_up / is1_up;
933  ibit1_down = list_1[j] & is1_down;
934  num1_down = ibit1_down / is1_down;
935 
936  list_Diagonal[j] += dtmp_V*(num1_up*num2_up + num1_down*num2_down);
937  }
938  break;/*case KondoGC, Hubbard, Kondo:*/
939 
940  case SpinGC:
941  is1_up = X->Def.Tpow[isite1 - 1];
942  is2_up = X->Def.Tpow[isite2 - 1];
943  ibit2_up = (unsigned long int)myrank & is2_up;
944 
945  if (ibit2_up == is2_up) {
946 #pragma omp parallel for default(none) shared(list_Diagonal) \
947 firstprivate(i_max, dtmp_V, is1_up) private(j, ibit1_up)
948  for (j = 1; j <= i_max; j++) {
949  ibit1_up = (j - 1) & is1_up;
950  if (ibit1_up == is1_up) {
951  list_Diagonal[j] += dtmp_V;
952  }
953  }
954  }
955  else if(ibit2_up == 0){
956 #pragma omp parallel for default(none) shared(list_Diagonal) \
957 firstprivate(i_max, dtmp_V, is1_up) private(j, ibit1_up)
958  for (j = 1; j <= i_max; j++) {
959  ibit1_up = (j - 1) & is1_up;
960  if (ibit1_up == 0) {
961  list_Diagonal[j] += dtmp_V;
962  }
963  }
964  }
965  break;/*case SpinGC:*/
966 
967  case Spin:
968  is1_up = X->Def.Tpow[isite1 - 1];
969  is2_up = X->Def.Tpow[isite2 - 1];
970  ibit2_up = (unsigned long int)myrank & is2_up;
971 
972  if (ibit2_up == is2_up) {
973 #pragma omp parallel for default(none) shared(list_1, list_Diagonal) \
974 firstprivate(i_max, dtmp_V, is1_up) private(j, ibit1_up)
975  for (j = 1; j <= i_max; j++) {
976  ibit1_up = list_1[j] & is1_up;
977  if (ibit1_up == is1_up) {
978  list_Diagonal[j] += dtmp_V;
979  }
980  }
981  }
982  else if (ibit2_up == 0) {
983 #pragma omp parallel for default(none) shared(list_1, list_Diagonal) \
984 firstprivate(i_max, dtmp_V, is1_up) private(j, ibit1_up)
985  for (j = 1; j <= i_max; j++) {
986  ibit1_up = list_1[j] & is1_up;
987  if (ibit1_up == 0) {
988  list_Diagonal[j] += dtmp_V;
989  }
990  }
991  }
992  break;/*case Spin:*/
993 
994  default:
995  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
996  return -1;
997 
998  }/*switch (X->Def.iCalcModel)*/
999 
1000  return 0;
1001 
1002  }/*else if (isite2 > X->Def.Nsite)*/
1003  else{
1004  switch (X->Def.iCalcModel){
1005  case HubbardGC: // list_1[j] -> j-1
1006  is1_up = X->Def.Tpow[2*isite1-2];
1007  is1_down = X->Def.Tpow[2*isite1-1];
1008  is2_up = X->Def.Tpow[2*isite2-2];
1009  is2_down = X->Def.Tpow[2*isite2-1];
1010 
1011 #pragma omp parallel for default(none) shared( list_Diagonal) firstprivate(i_max, dtmp_V, is1_up, is1_down, is2_up, is2_down) private(num1_up, num1_down, num2_up, num2_down, ibit1_up, ibit1_down, ibit2_up, ibit2_down)
1012  for(j = 1; j <= i_max;j++){
1013  num1_up=0;
1014  num1_down=0;
1015  num2_up=0;
1016  num2_down=0;
1017 
1018  ibit1_up=(j-1)&is1_up;
1019  num1_up=ibit1_up/is1_up;
1020  ibit1_down=(j-1)&is1_down;
1021  num1_down=ibit1_down/is1_down;
1022 
1023  ibit2_up=(j-1)&is2_up;
1024  num2_up=ibit2_up/is2_up;
1025  ibit2_down=(j-1)&is2_down;
1026  num2_down=ibit2_down/is2_down;
1027 
1028  list_Diagonal[j]+=dtmp_V*(num1_up*num2_up+num1_down*num2_down);
1029  }
1030  break;
1031  case KondoGC:
1032  case Hubbard:
1033  case Kondo:
1034  is1_up = X->Def.Tpow[2*isite1-2];
1035  is1_down = X->Def.Tpow[2*isite1-1];
1036  is2_up = X->Def.Tpow[2*isite2-2];
1037  is2_down = X->Def.Tpow[2*isite2-1];
1038 
1039 #pragma omp parallel for default(none) shared(list_1, list_Diagonal) firstprivate(i_max, dtmp_V, is1_up, is1_down, is2_up, is2_down) private(num1_up, num1_down, num2_up, num2_down, ibit1_up, ibit1_down, ibit2_up, ibit2_down)
1040  for(j = 1; j <= i_max;j++){
1041  num1_up=0;
1042  num1_down=0;
1043  num2_up=0;
1044  num2_down=0;
1045 
1046  ibit1_up=list_1[j]&is1_up;
1047  num1_up=ibit1_up/is1_up;
1048  ibit1_down=list_1[j]&is1_down;
1049  num1_down=ibit1_down/is1_down;
1050 
1051  ibit2_up=list_1[j]&is2_up;
1052  num2_up=ibit2_up/is2_up;
1053  ibit2_down=list_1[j]&is2_down;
1054  num2_down=ibit2_down/is2_down;
1055 
1056  list_Diagonal[j]+=dtmp_V*(num1_up*num2_up+num1_down*num2_down);
1057  }
1058  break;
1059 
1060  case SpinGC:
1061  is1_up = X->Def.Tpow[isite1-1];
1062  is2_up = X->Def.Tpow[isite2-1];
1063  is_up = is1_up+is2_up;
1064 #pragma omp parallel for default(none) shared(list_1, list_Diagonal) firstprivate(i_max, dtmp_V, is1_up, is2_up, is_up) private(j, ibit)
1065  for(j = 1;j <= i_max;j++){
1066  ibit = (j-1) & is_up;
1067  if(ibit == 0 || ibit == is_up){
1068  list_Diagonal[j]+= dtmp_V;
1069  }
1070  }
1071  break;
1072 
1073  case Spin:
1074  is1_up = X->Def.Tpow[isite1-1];
1075  is2_up = X->Def.Tpow[isite2-1];
1076  is_up = is1_up+is2_up;
1077 #pragma omp parallel for default(none) shared(list_1, list_Diagonal) firstprivate(i_max, dtmp_V, is1_up, is2_up, is_up) private(j, ibit)
1078  for(j = 1;j <= i_max;j++){
1079  ibit = list_1[j] & is_up;
1080  if(ibit == 0 || ibit == is_up){
1081  list_Diagonal[j]+= dtmp_V;
1082  }
1083  }
1084  break;
1085  default:
1086  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
1087  return -1;
1088  }
1089  }
1090  return 0;
1091 }
1092 
1109 int SetDiagonalInterAll
1110 (
1111  long unsigned int isite1,
1112  long unsigned int isite2,
1113  long unsigned int isigma1,
1114  long unsigned int isigma2,
1115  double dtmp_V,
1116  struct BindStruct *X
1117  )
1118 {
1119  long unsigned int is1_spin;
1120  long unsigned int is2_spin;
1121  long unsigned int is1_up;
1122  long unsigned int is2_up;
1123 
1124  long unsigned int ibit1_spin;
1125  long unsigned int ibit2_spin;
1126 
1127  long unsigned int num1;
1128  long unsigned int num2;
1129 
1130  long unsigned int j;
1131  long unsigned int i_max=X->Check.idim_max;
1132 
1133  /*
1134  Forse isite1 <= isite2
1135  */
1136  if (isite2 < isite1) {
1137  j = isite2;
1138  isite2 = isite1;
1139  isite1 = j;
1140  j = isigma2;
1141  isigma2 = isigma1;
1142  isigma1 = j;
1143  }
1144  /*
1145  When isite1 & site2 are in the inter process regino
1146  */
1147  if (isite1 > X->Def.Nsite) {
1148 
1149  switch (X->Def.iCalcModel) {
1150 
1151  case HubbardGC:
1152  case KondoGC:
1153  case Hubbard:
1154  case Kondo:
1155 
1156  is1_spin = X->Def.Tpow[2 * isite1 - 2 + isigma1];
1157  is2_spin = X->Def.Tpow[2 * isite2 - 2 + isigma2];
1158 
1159  num1 = 0;
1160  ibit1_spin = (unsigned long int)myrank&is1_spin;
1161  num1 += ibit1_spin / is1_spin;
1162 
1163  num2 = 0;
1164  ibit2_spin = (unsigned long int)myrank&is2_spin;
1165  num2 += ibit2_spin / is2_spin;
1166 
1167 #pragma omp parallel for default(none) shared(list_Diagonal) \
1168 firstprivate(i_max, dtmp_V, num2, num1) private(ibit1_spin, j)
1169  for (j = 1; j <= i_max; j++) list_Diagonal[j] += num1*num2*dtmp_V;
1170 
1171  break;/*case HubbardGC, KondoGC, Hubbard, Kondo:*/
1172 
1173  case SpinGC:
1174  case Spin:
1175 
1176  if (X->Def.iFlgGeneralSpin == FALSE) {
1177  is1_up = X->Def.Tpow[isite1 - 1];
1178  is2_up = X->Def.Tpow[isite2 - 1];
1179  num1 = X_SpinGC_CisAis((unsigned long int)myrank + 1, X, is1_up, isigma1);
1180  num2 = X_SpinGC_CisAis((unsigned long int)myrank + 1, X, is2_up, isigma2);
1181 
1182 #pragma omp parallel for default(none) shared(list_Diagonal) \
1183 firstprivate(i_max, dtmp_V, is1_up, isigma1, X, num1, num2) private(j)
1184  for (j = 1; j <= i_max; j++) {
1185  list_Diagonal[j] += num1*num2*dtmp_V;
1186  }
1187  }/*if (X->Def.iFlgGeneralSpin == FALSE)*/
1188  else {//start:generalspin
1189  num1 = BitCheckGeneral((unsigned long int)myrank, isite1, isigma1,
1190  X->Def.SiteToBit, X->Def.Tpow);
1191  num2 = BitCheckGeneral((unsigned long int)myrank, isite2, isigma2,
1192  X->Def.SiteToBit, X->Def.Tpow);
1193  if (num1 !=0 && num2 != 0) {
1194 #pragma omp parallel for default(none) shared(list_Diagonal) \
1195 firstprivate(i_max, dtmp_V, num1, X) private(j)
1196  for (j = 1; j <= i_max; j++) list_Diagonal[j] += dtmp_V*num1;
1197  }
1198  }/*if (X->Def.iFlgGeneralSpin == TRUE)*/
1199 
1200  break;/*case SpinGC, Spin:*/
1201 
1202  default:
1203  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
1204  return -1;
1205 
1206  }/*if (isite1 > X->Def.Nsite)*/
1207 
1208  return 0;
1209 
1210  }/*if (isite1 > X->Def.Nsite)*/
1211  else if (isite2 > X->Def.Nsite) {
1212 
1213  switch (X->Def.iCalcModel) {
1214 
1215  case HubbardGC:
1216 
1217  is1_spin = X->Def.Tpow[2 * isite1 - 2 + isigma1];
1218  is2_spin = X->Def.Tpow[2 * isite2 - 2 + isigma2];
1219 
1220  num2 = 0;
1221  ibit2_spin = (unsigned long int)myrank&is2_spin;
1222  num2 += ibit2_spin / is2_spin;
1223 
1224 #pragma omp parallel for default(none) shared(list_Diagonal) \
1225 firstprivate(i_max, dtmp_V, is1_spin, num2) private(num1, ibit1_spin, j)
1226  for (j = 1; j <= i_max; j++) {
1227  num1 = 0;
1228  ibit1_spin = (j - 1)&is1_spin;
1229  num1 += ibit1_spin / is1_spin;
1230  list_Diagonal[j] += num1*num2*dtmp_V;
1231  }
1232  break;/*case HubbardGC:*/
1233 
1234  case KondoGC:
1235  case Hubbard:
1236  case Kondo:
1237 
1238  is1_spin = X->Def.Tpow[2 * isite1 - 2 + isigma1];
1239  is2_spin = X->Def.Tpow[2 * isite2 - 2 + isigma2];
1240 
1241  num2 = 0;
1242  ibit2_spin = (unsigned long int)myrank&is2_spin;
1243  num2 += ibit2_spin / is2_spin;
1244 
1245 #pragma omp parallel for default(none) shared(list_Diagonal, list_1) \
1246 firstprivate(i_max, dtmp_V, is1_spin, num2) private(num1, ibit1_spin, j)
1247  for (j = 1; j <= i_max; j++) {
1248  num1 = 0;
1249  ibit1_spin = list_1[j] & is1_spin;
1250  num1 += ibit1_spin / is1_spin;
1251  list_Diagonal[j] += num1*num2*dtmp_V;
1252  }
1253  break;/*case KondoGC, Hubbard, Kondo:*/
1254 
1255  case SpinGC:
1256 
1257  if (X->Def.iFlgGeneralSpin == FALSE) {
1258  is1_up = X->Def.Tpow[isite1 - 1];
1259  is2_up = X->Def.Tpow[isite2 - 1];
1260  num2 = X_SpinGC_CisAis((unsigned long int)myrank + 1, X, is2_up, isigma2);
1261 
1262 #pragma omp parallel for default(none) shared(list_Diagonal) \
1263 firstprivate(i_max, dtmp_V, is1_up, isigma1, X, num2) private(j, num1)
1264  for (j = 1; j <= i_max; j++) {
1265  num1 = X_SpinGC_CisAis(j, X, is1_up, isigma1);
1266  list_Diagonal[j] += num1*num2*dtmp_V;
1267  }
1268  }/* if (X->Def.iFlgGeneralSpin == FALSE)*/
1269  else {//start:generalspin
1270  num2 = BitCheckGeneral((unsigned long int)myrank, isite2, isigma2,
1271  X->Def.SiteToBit, X->Def.Tpow);
1272  if (num2 != 0) {
1273 #pragma omp parallel for default(none) shared(list_Diagonal) \
1274 firstprivate(i_max, dtmp_V, isite1, isigma1, X) private(j, num1)
1275  for (j = 1; j <= i_max; j++) {
1276  num1 = BitCheckGeneral(j - 1, isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
1277  list_Diagonal[j] += dtmp_V*num1;
1278  }
1279  }
1280  }/* if (X->Def.iFlgGeneralSpin == TRUE)*/
1281 
1282  break;/*case SpinGC:*/
1283 
1284  case Spin:
1285 
1286  if (X->Def.iFlgGeneralSpin == FALSE) {
1287  is1_up = X->Def.Tpow[isite1 - 1];
1288  is2_up = X->Def.Tpow[isite2 - 1];
1289  num2 = X_SpinGC_CisAis((unsigned long int)myrank + 1, X, is2_up, isigma2);
1290 
1291 #pragma omp parallel for default(none) shared(list_Diagonal) \
1292 firstprivate(i_max, dtmp_V, is1_up, isigma1, X, num2) private(j, num1)
1293  for (j = 1; j <= i_max; j++) {
1294  num1 = X_Spin_CisAis(j, X, is1_up, isigma1);
1295  list_Diagonal[j] += num1*num2*dtmp_V;
1296  }
1297  }/* if (X->Def.iFlgGeneralSpin == FALSE)*/
1298  else /* if (X->Def.iFlgGeneralSpin == TRUE)*/{
1299  num2 = BitCheckGeneral((unsigned long int)myrank, isite2, isigma2, \
1300  X->Def.SiteToBit, X->Def.Tpow);
1301  if (num2 != 0) {
1302 #pragma omp parallel for default(none) shared(list_Diagonal, list_1) \
1303 firstprivate(i_max, dtmp_V, isite1, isigma1, X) private(j, num1)
1304  for (j = 1; j <= i_max; j++) {
1305  num1 = BitCheckGeneral(list_1[j], isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
1306  list_Diagonal[j] += dtmp_V*num1;
1307  }
1308  }
1309  } /* if (X->Def.iFlgGeneralSpin == TRUE)*/
1310 
1311  break;/*case Spin:*/
1312 
1313  default:
1314  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
1315  return -1;
1316 
1317  }/*switch (X->Def.iCalcModel)*/
1318 
1319  return 0;
1320 
1321  }/*else if (isite2 > X->Def.Nsite)*/
1322 
1323  switch (X->Def.iCalcModel){
1324  case HubbardGC: //list_1[j] -> j-1
1325  is1_spin = X->Def.Tpow[2*isite1-2+isigma1];
1326  is2_spin = X->Def.Tpow[2*isite2-2+isigma2];
1327 #pragma omp parallel for default(none) shared(list_Diagonal) firstprivate(i_max, dtmp_V, is1_spin, is2_spin) private(num1, ibit1_spin, num2, ibit2_spin)
1328  for(j = 1;j <= i_max;j++){
1329  num1=0;
1330  num2=0;
1331  ibit1_spin=(j-1)&is1_spin;
1332  num1+=ibit1_spin/is1_spin;
1333  ibit2_spin=(j-1)&is2_spin;
1334  num2+=ibit2_spin/is2_spin;
1335  list_Diagonal[j]+=num1*num2*dtmp_V;
1336  }
1337  break;
1338  case KondoGC:
1339  case Hubbard:
1340  case Kondo:
1341  is1_spin = X->Def.Tpow[2*isite1-2+isigma1];
1342  is2_spin = X->Def.Tpow[2*isite2-2+isigma2];
1343 
1344 #pragma omp parallel for default(none) shared(list_Diagonal, list_1) firstprivate(i_max, dtmp_V, is1_spin, is2_spin) private(num1, ibit1_spin, num2, ibit2_spin)
1345  for(j = 1;j <= i_max;j++){
1346  num1=0;
1347  num2=0;
1348  ibit1_spin=list_1[j]&is1_spin;
1349  num1+=ibit1_spin/is1_spin;
1350 
1351  ibit2_spin=list_1[j]&is2_spin;
1352  num2+=ibit2_spin/is2_spin;
1353  list_Diagonal[j]+=num1*num2*dtmp_V;
1354  }
1355  break;
1356 
1357  case Spin:
1358  if(X->Def.iFlgGeneralSpin==FALSE){
1359  is1_up = X->Def.Tpow[isite1-1];
1360  is2_up = X->Def.Tpow[isite2-1];
1361 #pragma omp parallel for default(none) shared(list_Diagonal) firstprivate(i_max, dtmp_V, is1_up, is2_up, isigma1, isigma2, X) private(j, num1, num2)
1362  for(j = 1;j <= i_max; j++){
1363  num1=X_Spin_CisAis(j, X, is1_up, isigma1);
1364  num2=X_Spin_CisAis(j, X, is2_up, isigma2);
1365  list_Diagonal[j] += num1*num2*dtmp_V;
1366  }
1367  }
1368  else{
1369 #pragma omp parallel for default(none) shared(list_Diagonal, list_1) firstprivate(i_max, dtmp_V, isite1, isite2, isigma1, isigma2, X) private(j, num1)
1370  for(j = 1;j <= i_max; j++){
1371  num1=BitCheckGeneral (list_1[j], isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
1372  if(num1 != 0){
1373  num1=BitCheckGeneral (list_1[j], isite2, isigma2, X->Def.SiteToBit, X->Def.Tpow);
1374  list_Diagonal[j] += dtmp_V*num1;
1375  }
1376  }
1377 
1378  }
1379  break;
1380 
1381  case SpinGC:
1382  if(X->Def.iFlgGeneralSpin==FALSE){
1383  is1_up = X->Def.Tpow[isite1-1];
1384  is2_up = X->Def.Tpow[isite2-1];
1385 #pragma omp parallel for default(none) shared(list_Diagonal) firstprivate(i_max, dtmp_V, is1_up, is2_up, isigma1, isigma2, X) private(j, num1, num2)
1386  for(j = 1;j <= i_max; j++){
1387  num1=X_SpinGC_CisAis(j, X, is1_up, isigma1);
1388  num2=X_SpinGC_CisAis(j, X, is2_up, isigma2);
1389  list_Diagonal[j] += num1*num2*dtmp_V;
1390  }
1391  }
1392  else{//start:generalspin
1393 #pragma omp parallel for default(none) shared(list_Diagonal) firstprivate(i_max, dtmp_V, isite1, isite2, isigma1, isigma2, X) private(j, num1)
1394  for(j = 1;j <= i_max; j++){
1395  num1=BitCheckGeneral (j-1, isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
1396  if(num1 != 0){
1397  num1=BitCheckGeneral (j-1, isite2, isigma2, X->Def.SiteToBit, X->Def.Tpow);
1398  list_Diagonal[j] += dtmp_V*num1;
1399  }
1400  }
1401  }
1402  break;
1403 
1404  default:
1405  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
1406  return -1;
1407  }
1408 
1409  return 0;
1410 
1411 }
1412 
1431 int SetDiagonalTEInterAll(
1432  long unsigned int isite1,
1433  long unsigned int isite2,
1434  long unsigned int isigma1,
1435  long unsigned int isigma2,
1436  double dtmp_V,
1437  struct BindStruct *X,
1438  double complex *tmp_v0,
1439  double complex *tmp_v1
1440 ){
1441  long unsigned int is1_spin;
1442  long unsigned int is2_spin;
1443  long unsigned int is1_up;
1444  long unsigned int is2_up;
1445 
1446  long unsigned int ibit1_spin;
1447  long unsigned int ibit2_spin;
1448 
1449  long unsigned int num1;
1450  long unsigned int num2;
1451 
1452  long unsigned int j;
1453  long unsigned int i_max=X->Check.idim_max;
1454  double complex dam_pr=0.0;
1455 
1456 
1457  /*
1458  Forse isite1 <= isite2
1459  */
1460  if (isite2 < isite1) {
1461  j = isite2;
1462  isite2 = isite1;
1463  isite1 = j;
1464  j = isigma2;
1465  isigma2 = isigma1;
1466  isigma1 = j;
1467  }
1468  /*
1469  When isite1 & site2 are in the inter process regino
1470  */
1471  if (isite1 > X->Def.Nsite) {
1472 
1473  switch (X->Def.iCalcModel) {
1474 
1475  case HubbardGC:
1476  case KondoGC:
1477  case Hubbard:
1478  case Kondo:
1479  is1_spin = X->Def.Tpow[2 * isite1 - 2 + isigma1];
1480  is2_spin = X->Def.Tpow[2 * isite2 - 2 + isigma2];
1481  num1 = 0;
1482  ibit1_spin = (unsigned long int)myrank&is1_spin;
1483  num1 += ibit1_spin / is1_spin;
1484  num2 = 0;
1485  ibit2_spin = (unsigned long int)myrank&is2_spin;
1486  num2 += ibit2_spin / is2_spin;
1487  break;/*case HubbardGC, KondoGC, Hubbard, Kondo:*/
1488 
1489  case SpinGC:
1490  case Spin:
1491  if (X->Def.iFlgGeneralSpin == FALSE) {
1492  is1_up = X->Def.Tpow[isite1 - 1];
1493  is2_up = X->Def.Tpow[isite2 - 1];
1494  num1 = X_SpinGC_CisAis((unsigned long int) myrank + 1, X, is1_up, isigma1);
1495  num2 = X_SpinGC_CisAis((unsigned long int) myrank + 1, X, is2_up, isigma2);
1496  }/*if (X->Def.iFlgGeneralSpin == FALSE)*/
1497  else {//start:generalspin
1498  num1 = BitCheckGeneral((unsigned long int) myrank, isite1, isigma1,
1499  X->Def.SiteToBit, X->Def.Tpow);
1500  num2 = BitCheckGeneral((unsigned long int) myrank, isite2, isigma2,
1501  X->Def.SiteToBit, X->Def.Tpow);
1502  }
1503  break;/*case SpinGC, Spin:*/
1504 
1505  default:
1506  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
1507  return -1;
1508  }/*if (isite1 > X->Def.Nsite)*/
1509 
1510  if (num1 * num2 != 0) {
1511 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1) \
1512 firstprivate(i_max, dtmp_V) private(j)
1513  for (j = 1; j <= i_max; j++) {
1514  tmp_v0[j] += dtmp_V * tmp_v1[j];
1515  dam_pr += dtmp_V * conj(tmp_v1[j]) * tmp_v1[j];
1516  }
1517  }
1518  dam_pr=SumMPI_dc(dam_pr);
1519  X->Large.prdct += dam_pr;
1520  return 0;
1521 
1522  }/*if (isite1 > X->Def.Nsite)*/
1523  else if (isite2 > X->Def.Nsite) {
1524 
1525  switch (X->Def.iCalcModel) {
1526 
1527  case HubbardGC:
1528 
1529  is1_spin = X->Def.Tpow[2 * isite1 - 2 + isigma1];
1530  is2_spin = X->Def.Tpow[2 * isite2 - 2 + isigma2];
1531 
1532  num2 = 0;
1533  ibit2_spin = (unsigned long int)myrank&is2_spin;
1534  num2 += ibit2_spin / is2_spin;
1535  if(num2 !=0) {
1536 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1)\
1537  firstprivate(i_max, dtmp_V, is1_spin) private(num1, ibit1_spin, j)
1538  for (j = 1; j <= i_max; j++) {
1539  num1 = 0;
1540  ibit1_spin = (j - 1) & is1_spin;
1541  num1 += ibit1_spin / is1_spin;
1542  tmp_v0[j] += dtmp_V * num1 * tmp_v1[j];
1543  dam_pr += dtmp_V * num1 * conj(tmp_v1[j]) * tmp_v1[j];
1544  }
1545  }
1546  break;/*case HubbardGC:*/
1547 
1548  case KondoGC:
1549  case Hubbard:
1550  case Kondo:
1551 
1552  is1_spin = X->Def.Tpow[2 * isite1 - 2 + isigma1];
1553  is2_spin = X->Def.Tpow[2 * isite2 - 2 + isigma2];
1554 
1555  num2 = 0;
1556  ibit2_spin = (unsigned long int)myrank&is2_spin;
1557  num2 += ibit2_spin / is2_spin;
1558  if(num2 !=0) {
1559 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1, list_1)\
1560  firstprivate(i_max, dtmp_V, is1_spin) private(num1, ibit1_spin, j)
1561  for (j = 1; j <= i_max; j++) {
1562  num1 = 0;
1563  ibit1_spin = list_1[j] & is1_spin;
1564  num1 += ibit1_spin / is1_spin;
1565  tmp_v0[j] += dtmp_V *num1*tmp_v1[j];
1566  dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j];
1567  }
1568  }
1569  break;/*case KondoGC, Hubbard, Kondo:*/
1570 
1571  case SpinGC:
1572 
1573  if (X->Def.iFlgGeneralSpin == FALSE) {
1574  is1_up = X->Def.Tpow[isite1 - 1];
1575  is2_up = X->Def.Tpow[isite2 - 1];
1576  num2 = X_SpinGC_CisAis((unsigned long int)myrank + 1, X, is2_up, isigma2);
1577 
1578  if(num2 !=0) {
1579 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1)\
1580  firstprivate(i_max, dtmp_V, is1_up, isigma1, X) private(num1, j)
1581  for (j = 1; j <= i_max; j++) {
1582  num1 = X_SpinGC_CisAis(j, X, is1_up, isigma1);
1583  tmp_v0[j] += dtmp_V * num1 * tmp_v1[j];
1584  dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j];
1585  }
1586  }
1587  }/* if (X->Def.iFlgGeneralSpin == FALSE)*/
1588  else {//start:generalspin
1589  num2 = BitCheckGeneral((unsigned long int)myrank, isite2, isigma2,
1590  X->Def.SiteToBit, X->Def.Tpow);
1591  if (num2 != 0) {
1592 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1) \
1593 firstprivate(i_max, dtmp_V, isite1, isigma1, X) private(j, num1)
1594  for (j = 1; j <= i_max; j++) {
1595  num1 = BitCheckGeneral(j - 1, isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
1596  tmp_v0[j] += dtmp_V * num1 * tmp_v1[j];
1597  dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j];
1598  }
1599  }
1600  }/* if (X->Def.iFlgGeneralSpin == TRUE)*/
1601 
1602  break;/*case SpinGC:*/
1603 
1604  case Spin:
1605 
1606  if (X->Def.iFlgGeneralSpin == FALSE) {
1607  is1_up = X->Def.Tpow[isite1 - 1];
1608  is2_up = X->Def.Tpow[isite2 - 1];
1609  num2 = X_SpinGC_CisAis((unsigned long int)myrank + 1, X, is2_up, isigma2);
1610 
1611  if(num2 !=0) {
1612 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1) \
1613 firstprivate(i_max, dtmp_V, is1_up, isigma1, X, num2) private(j, num1)
1614  for (j = 1; j <= i_max; j++) {
1615  num1 = X_Spin_CisAis(j, X, is1_up, isigma1);
1616  tmp_v0[j] += dtmp_V * num1 * tmp_v1[j];
1617  dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j];
1618  }
1619  }
1620  }/* if (X->Def.iFlgGeneralSpin == FALSE)*/
1621  else /* if (X->Def.iFlgGeneralSpin == TRUE)*/{
1622  num2 = BitCheckGeneral((unsigned long int)myrank, isite2, isigma2, \
1623  X->Def.SiteToBit, X->Def.Tpow);
1624  if (num2 != 0) {
1625 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1, list_1)\
1626 firstprivate(i_max, dtmp_V, isite1, isigma1, X) private(j, num1)
1627  for (j = 1; j <= i_max; j++) {
1628  num1 = BitCheckGeneral(list_1[j], isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
1629  tmp_v0[j] += dtmp_V * num1 * tmp_v1[j];
1630  dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j];
1631  }
1632  }
1633  } /* if (X->Def.iFlgGeneralSpin == TRUE)*/
1634 
1635  break;/*case Spin:*/
1636 
1637  default:
1638  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
1639  return -1;
1640 
1641  }/*switch (X->Def.iCalcModel)*/
1642  dam_pr=SumMPI_dc(dam_pr);
1643  X->Large.prdct += dam_pr;
1644  return 0;
1645  }/*else if (isite2 > X->Def.Nsite)*/
1646 
1647  switch (X->Def.iCalcModel){
1648  case HubbardGC: //list_1[j] -> j-1
1649  is1_spin = X->Def.Tpow[2*isite1-2+isigma1];
1650  is2_spin = X->Def.Tpow[2*isite2-2+isigma2];
1651 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1) firstprivate(i_max, dtmp_V, is1_spin, is2_spin) private(num1, ibit1_spin, num2, ibit2_spin)
1652  for(j = 1;j <= i_max;j++){
1653  num1=0;
1654  num2=0;
1655  ibit1_spin=(j-1)&is1_spin;
1656  num1+=ibit1_spin/is1_spin;
1657  ibit2_spin=(j-1)&is2_spin;
1658  num2+=ibit2_spin/is2_spin;
1659  tmp_v0[j] += dtmp_V * num1*num2*tmp_v1[j];
1660  dam_pr += dtmp_V * num1*num2*conj(tmp_v1[j]) * tmp_v1[j];
1661  }
1662  break;
1663  case KondoGC:
1664  case Hubbard:
1665  case Kondo:
1666  is1_spin = X->Def.Tpow[2*isite1-2+isigma1];
1667  is2_spin = X->Def.Tpow[2*isite2-2+isigma2];
1668 
1669 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1, list_1) firstprivate(i_max, dtmp_V, is1_spin, is2_spin) private(num1, ibit1_spin, num2, ibit2_spin)
1670  for(j = 1;j <= i_max;j++){
1671  num1=0;
1672  num2=0;
1673  ibit1_spin=list_1[j]&is1_spin;
1674  num1+=ibit1_spin/is1_spin;
1675 
1676  ibit2_spin=list_1[j]&is2_spin;
1677  num2+=ibit2_spin/is2_spin;
1678  tmp_v0[j] += dtmp_V * num1*num2*tmp_v1[j];
1679  dam_pr += dtmp_V * num1*num2*conj(tmp_v1[j]) * tmp_v1[j];
1680  }
1681  break;
1682 
1683  case Spin:
1684  if(X->Def.iFlgGeneralSpin==FALSE){
1685  is1_up = X->Def.Tpow[isite1-1];
1686  is2_up = X->Def.Tpow[isite2-1];
1687 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1) firstprivate(i_max, dtmp_V, is1_up, is2_up, isigma1, isigma2, X) private(j, num1, num2)
1688  for(j = 1;j <= i_max; j++){
1689  num1=X_Spin_CisAis(j, X, is1_up, isigma1);
1690  num2=X_Spin_CisAis(j, X, is2_up, isigma2);
1691  tmp_v0[j] += dtmp_V * num1*num2*tmp_v1[j];
1692  dam_pr += dtmp_V * num1*num2*conj(tmp_v1[j]) * tmp_v1[j];
1693  }
1694  }
1695  else{
1696 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1, list_1) firstprivate(i_max, dtmp_V, isite1, isite2, isigma1, isigma2, X) private(j, num1)
1697  for(j = 1;j <= i_max; j++){
1698  num1=BitCheckGeneral (list_1[j], isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
1699  if(num1 != 0){
1700  num1=BitCheckGeneral (list_1[j], isite2, isigma2, X->Def.SiteToBit, X->Def.Tpow);
1701  tmp_v0[j] += dtmp_V *num1*tmp_v1[j];
1702  dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j];
1703  }
1704  }
1705 
1706  }
1707  break;
1708 
1709  case SpinGC:
1710  if(X->Def.iFlgGeneralSpin==FALSE){
1711  is1_up = X->Def.Tpow[isite1-1];
1712  is2_up = X->Def.Tpow[isite2-1];
1713 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1) firstprivate(i_max, dtmp_V, is1_up, is2_up, isigma1, isigma2, X) private(j, num1, num2)
1714  for(j = 1;j <= i_max; j++){
1715  num1=X_SpinGC_CisAis(j, X, is1_up, isigma1);
1716  num2=X_SpinGC_CisAis(j, X, is2_up, isigma2);
1717  tmp_v0[j] += dtmp_V * num1*num2*tmp_v1[j];
1718  dam_pr += dtmp_V * num1*num2*conj(tmp_v1[j]) * tmp_v1[j];
1719  }
1720  }
1721  else{//start:generalspin
1722 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1) firstprivate(i_max, dtmp_V, isite1, isite2, isigma1, isigma2, X) private(j, num1)
1723  for(j = 1;j <= i_max; j++){
1724  num1=BitCheckGeneral (j-1, isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
1725  if(num1 != 0){
1726  num1=BitCheckGeneral (j-1, isite2, isigma2, X->Def.SiteToBit, X->Def.Tpow);
1727  tmp_v0[j] += dtmp_V *num1*tmp_v1[j];
1728  dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j];
1729  }
1730  }
1731  }
1732  break;
1733 
1734  default:
1735  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
1736  return -1;
1737  }
1738  dam_pr=SumMPI_dc(dam_pr);
1739  X->Large.prdct += dam_pr;
1740  return 0;
1741 }
1742 
1762 int SetDiagonalTEChemi(
1763  long unsigned int isite1,
1764  long unsigned int spin,
1765  double dtmp_V,
1766  struct BindStruct *X,
1767  double complex *tmp_v0,
1768  double complex *tmp_v1
1769 ){
1770  long unsigned int is1_up;
1771  long unsigned int num1;
1772  long unsigned int isigma1 =spin;
1773  long unsigned int is1,ibit1;
1774 
1775  long unsigned int j;
1776  long unsigned int i_max=X->Check.idim_max;
1777  double complex dam_pr=0;
1778 
1779  /*
1780  When isite1 is in the inter process region
1781  */
1782  if (isite1 > X->Def.Nsite){
1783 
1784  switch (X->Def.iCalcModel) {
1785 
1786  case HubbardGC:
1787  case KondoGC:
1788  case Hubbard:
1789  case Kondo:
1790 
1791  if (spin == 0) {
1792  is1 = X->Def.Tpow[2 * isite1 - 2];
1793  }
1794  else {
1795  is1 = X->Def.Tpow[2 * isite1 - 1];
1796  }
1797  ibit1 = (unsigned long int)myrank & is1;
1798  num1 = ibit1 / is1;
1799  break;/*case HubbardGC, case KondoGC, Hubbard, Kondo:*/
1800 
1801  case SpinGC:
1802  case Spin:
1803 
1804  if (X->Def.iFlgGeneralSpin == FALSE) {
1805  is1_up = X->Def.Tpow[isite1 - 1];
1806  num1 = (((unsigned long int)myrank& is1_up) / is1_up) ^ (1 - spin);
1807  } /*if (X->Def.iFlgGeneralSpin == FALSE)*/
1808  else /*if (X->Def.iFlgGeneralSpin == TRUE)*/ {
1809  num1 = BitCheckGeneral((unsigned long int)myrank,
1810  isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
1811  }/*if (X->Def.iFlgGeneralSpin == TRUE)*/
1812  break;/*case SpinGC, Spin:*/
1813 
1814  default:
1815  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
1816  return -1;
1817 
1818  } /*switch (X->Def.iCalcModel)*/
1819  if (num1 != 0) {
1820 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1) \
1821 firstprivate(i_max, dtmp_V) private(j)
1822  for (j = 1; j <= i_max; j++){
1823  tmp_v0[j] += dtmp_V * tmp_v1[j];
1824  dam_pr += dtmp_V * conj(tmp_v1[j]) * tmp_v1[j];
1825  }
1826  }/*if (num1 != 0)*/
1827  dam_pr=SumMPI_dc(dam_pr);
1828  X->Large.prdct += dam_pr;
1829  return 0;
1830 
1831  }/*if (isite1 >= X->Def.Nsite*/
1832 
1833  switch (X->Def.iCalcModel){
1834  case HubbardGC:
1835  if(spin==0){
1836  is1 = X->Def.Tpow[2*isite1-2];
1837  }else{
1838  is1 = X->Def.Tpow[2*isite1-1];
1839  }
1840 
1841  #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1) firstprivate(i_max, dtmp_V, is1) private(num1, ibit1)
1842  for(j = 1;j <= i_max;j++){
1843  ibit1 = (j-1)&is1;
1844  num1 = ibit1/is1;
1845  tmp_v0[j] += dtmp_V * num1*tmp_v1[j];
1846  dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j];
1847  }
1848  break;
1849  case KondoGC:
1850  case Hubbard:
1851  case Kondo:
1852  if(spin==0){
1853  is1 = X->Def.Tpow[2*isite1-2];
1854  }else{
1855  is1 = X->Def.Tpow[2*isite1-1];
1856  }
1857 
1858 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(list_1, tmp_v0, tmp_v1) firstprivate(i_max, dtmp_V, is1) private(num1, ibit1)
1859  for(j = 1;j <= i_max;j++){
1860 
1861  ibit1 = list_1[j]&is1;
1862  num1 = ibit1/is1;
1863  tmp_v0[j] += dtmp_V * num1*tmp_v1[j];
1864  dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j];
1865  }
1866  break;
1867 
1868  case SpinGC:
1869  if(X->Def.iFlgGeneralSpin==FALSE){
1870  is1_up = X->Def.Tpow[isite1-1];
1871 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(list_1, tmp_v0, tmp_v1) firstprivate(i_max, dtmp_V, is1_up, spin) private(num1)
1872  for(j = 1;j <= i_max;j++){
1873  num1=(((j-1)& is1_up)/is1_up)^(1-spin);
1874  tmp_v0[j] += dtmp_V * num1*tmp_v1[j];
1875  dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j];
1876  }
1877  }
1878  else{
1879 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1) firstprivate(i_max, dtmp_V, isite1, isigma1, X) private(j, num1)
1880  for(j = 1;j <= i_max; j++){
1881  num1=BitCheckGeneral (j-1, isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
1882  if(num1 != 0){
1883  tmp_v0[j] += dtmp_V * tmp_v1[j];
1884  dam_pr += dtmp_V * conj(tmp_v1[j]) * tmp_v1[j];
1885  }
1886  }
1887  }
1888  break;
1889 
1890  case Spin:
1891  if(X->Def.iFlgGeneralSpin==FALSE){
1892  is1_up = X->Def.Tpow[isite1-1];
1893 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(list_1, tmp_v0, tmp_v1) firstprivate(i_max, dtmp_V, is1_up, spin) private(num1)
1894  for(j = 1;j <= i_max;j++){
1895  num1=((list_1[j]& is1_up)/is1_up)^(1-spin);
1896  tmp_v0[j] += dtmp_V * num1*tmp_v1[j];
1897  dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j];
1898  }
1899  }
1900  else{
1901 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1, list_1) firstprivate(i_max, dtmp_V, isite1, isigma1, X) private(j, num1)
1902  for(j = 1;j <= i_max; j++){
1903  num1=BitCheckGeneral (list_1[j], isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
1904  if(num1 != 0){
1905  tmp_v0[j] += dtmp_V * tmp_v1[j];
1906  dam_pr += dtmp_V * conj(tmp_v1[j]) * tmp_v1[j];
1907 
1908  }
1909  }
1910  }
1911 
1912  break;
1913  default:
1914  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
1915  return -1;
1916  }
1917  dam_pr=SumMPI_dc(dam_pr);
1918  X->Large.prdct += dam_pr;
1919  return 0;
1920 }
1921 
1939 int SetDiagonalTETransfer
1940  (
1941  long unsigned int isite1,
1942  double dtmp_V,
1943  long unsigned int spin,
1944  struct BindStruct *X,
1945  double complex *tmp_v0,
1946  double complex *tmp_v1
1947  ){
1948  long unsigned int is1_up;
1949  long unsigned int ibit1_up;
1950  long unsigned int num1;
1951  long unsigned int isigma1 =spin;
1952  long unsigned int is1,ibit1;
1953  double dam_pr=0.0;
1954 
1955  long unsigned int j;
1956  long unsigned int i_max=X->Check.idim_max;
1957 
1958  /*
1959  When isite1 is in the inter process region
1960  */
1961  if (isite1 > X->Def.Nsite){
1962 
1963  switch (X->Def.iCalcModel) {
1964 
1965  case HubbardGC:
1966  case KondoGC:
1967  case Hubbard:
1968  case Kondo:
1969  if (spin == 0) {
1970  is1 = X->Def.Tpow[2 * isite1 - 2];
1971  }
1972  else {
1973  is1 = X->Def.Tpow[2 * isite1 - 1];
1974  }
1975  ibit1 = (unsigned long int)myrank & is1;
1976  num1 = ibit1 / is1;
1977  break;/*case HubbardGC, case KondoGC, Hubbard, Kondo:*/
1978 
1979  case SpinGC:
1980  case Spin:
1981  if (X->Def.iFlgGeneralSpin == FALSE) {
1982  is1_up = X->Def.Tpow[isite1 - 1];
1983  num1 = (((unsigned long int)myrank& is1_up) / is1_up) ^ (1 - spin);
1984  } /*if (X->Def.iFlgGeneralSpin == FALSE)*/
1985  else /*if (X->Def.iFlgGeneralSpin == TRUE)*/ {
1986  num1 = BitCheckGeneral((unsigned long int)myrank,
1987  isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
1988  }/*if (X->Def.iFlgGeneralSpin == TRUE)*/
1989  break;/*case SpinGC, Spin:*/
1990 
1991  default:
1992  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
1993  return -1;
1994 
1995  } /*switch (X->Def.iCalcModel)*/
1996 
1997  if(num1 !=0) {
1998 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1)\
1999  firstprivate(i_max, dtmp_V) private(j)
2000  for (j = 1; j <= i_max; j++) {
2001  tmp_v0[j] += dtmp_V * tmp_v1[j];
2002  dam_pr += dtmp_V * conj(tmp_v1[j]) * tmp_v1[j];
2003  }
2004  }
2005  }/*if (isite1 >= X->Def.Nsite*/
2006  else {//(isite1 < X->Def.Nsite)
2007  switch (X->Def.iCalcModel) {
2008  case HubbardGC:
2009  if (spin == 0) {
2010  is1 = X->Def.Tpow[2 * isite1 - 2];
2011  } else {
2012  is1 = X->Def.Tpow[2 * isite1 - 1];
2013  }
2014 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(list_1, tmp_v0, tmp_v1) \
2015  firstprivate(i_max, dtmp_V, is1) private(num1, ibit1)
2016  for (j = 1; j <= i_max; j++) {
2017  ibit1 = (j - 1) & is1;
2018  num1 = ibit1 / is1;
2019  tmp_v0[j] += dtmp_V * num1*tmp_v1[j];
2020  dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j];
2021  }
2022  break;
2023 
2024  case KondoGC:
2025  case Hubbard:
2026  case Kondo:
2027  if (spin == 0) {
2028  is1 = X->Def.Tpow[2 * isite1 - 2];
2029  } else {
2030  is1 = X->Def.Tpow[2 * isite1 - 1];
2031  }
2032 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(list_1, tmp_v0, tmp_v1) \
2033  firstprivate(i_max, dtmp_V, is1) private(num1, ibit1)
2034  for (j = 1; j <= i_max; j++) {
2035  ibit1 = list_1[j] & is1;
2036  num1 = ibit1 / is1;
2037  tmp_v0[j] += dtmp_V * num1*tmp_v1[j];
2038  dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j];
2039  }
2040  break;
2041 
2042  case SpinGC:
2043  if (X->Def.iFlgGeneralSpin == FALSE) {
2044  is1_up = X->Def.Tpow[isite1 - 1];
2045 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(list_1, tmp_v0, tmp_v1) \
2046  firstprivate(i_max, dtmp_V, is1_up, spin) private(num1, ibit1_up)
2047  for (j = 1; j <= i_max; j++) {
2048  ibit1_up = (((j - 1) & is1_up) / is1_up) ^ (1 - spin);
2049  tmp_v0[j] += dtmp_V * ibit1_up*tmp_v1[j];
2050  dam_pr += dtmp_V * ibit1_up*conj(tmp_v1[j]) * tmp_v1[j];
2051  }
2052  } else {
2053 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(tmp_v0, tmp_v1) \
2054  firstprivate(i_max, dtmp_V, isite1, isigma1, X) private(j, num1)
2055  for (j = 1; j <= i_max; j++) {
2056  num1 = BitCheckGeneral(j - 1, isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
2057  if (num1 != 0) {
2058  tmp_v0[j] += dtmp_V *tmp_v1[j];
2059  dam_pr += dtmp_V *conj(tmp_v1[j]) * tmp_v1[j];
2060  }
2061  }
2062  }
2063  break;
2064 
2065  case Spin:
2066  if (X->Def.iFlgGeneralSpin == FALSE) {
2067  is1_up = X->Def.Tpow[isite1 - 1];
2068 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(list_1, tmp_v0, tmp_v1)\
2069  firstprivate(i_max, dtmp_V, is1_up, spin) private(num1, ibit1_up)
2070  for (j = 1; j <= i_max; j++) {
2071  ibit1_up = ((list_1[j] & is1_up) / is1_up) ^ (1 - spin);
2072  tmp_v0[j] += dtmp_V * ibit1_up * tmp_v1[j];
2073  dam_pr += dtmp_V * ibit1_up * conj(tmp_v1[j]) * tmp_v1[j];
2074  }
2075  } else {
2076 #pragma omp parallel for default(none) reduction(+:dam_pr) shared(list_1, tmp_v0, tmp_v1)\
2077  firstprivate(i_max, dtmp_V, isite1, isigma1, X) private(j, num1)
2078  for (j = 1; j <= i_max; j++) {
2079  num1 = BitCheckGeneral(list_1[j], isite1, isigma1, X->Def.SiteToBit, X->Def.Tpow);
2080  tmp_v0[j] += dtmp_V * num1 * tmp_v1[j];
2081  dam_pr += dtmp_V * num1 * conj(tmp_v1[j]) * tmp_v1[j];
2082  }
2083  }
2084  break;
2085 
2086  default:
2087  fprintf(stdoutMPI, cErrNoModel, X->Def.iCalcModel);
2088  return -1;
2089  }
2090  }
2091  dam_pr=SumMPI_dc(dam_pr);
2092  X->Large.prdct += dam_pr;
2093  return 0;
2094 }
cFileNameCheckChemi
const char * cFileNameCheckChemi
Definition: global.c:28
BindStruct::Def
struct DefineList Def
Definision of system (Hamiltonian) etc.
Definition: struct.h:410
DefineList::NHundCoupling
unsigned int NHundCoupling
Number of Hund coupling.
Definition: struct.h:133
DefineList::EDSpinChemi
int * EDSpinChemi
[DefineList::Nsite]
Definition: struct.h:99
DefineList::SpinTEChemi
int ** SpinTEChemi
Definition: struct.h:295
SetDiagonalTEChemi
int SetDiagonalTEChemi(long unsigned int isite1, long unsigned int spin, double dtmp_V, struct BindStruct *X, double complex *tmp_v0, double complex *tmp_v1)
Update the vector by the chemical potential generated by the commutation relation in terms of the g...
Definition: diagonalcalc.c:1762
SumMPI_dc
double complex SumMPI_dc(double complex norm)
MPI wrapper function to obtain sum of Double complex across processes.
Definition: wrapperMPI.c:205
CheckList::idim_max
unsigned long int idim_max
The dimension of the Hilbert space of this process.
Definition: struct.h:303
cDiagonalCalcStart
const char * cDiagonalCalcStart
Definition: LogMessage.c:56
DefineList::ParaInterAll_Diagonal
double * ParaInterAll_Diagonal
[DefineList::NInterAll_Diagonal] Coupling constant of diagonal inter-all term. malloc in setmem_def()...
Definition: struct.h:168
cFileNameCheckHund
const char * cFileNameCheckHund
Definition: global.c:30
DefineList::TETransferDiagonal
int *** TETransferDiagonal
Definition: struct.h:262
LargeList::prdct
double complex prdct
The expectation value of the energy.
Definition: struct.h:314
X_SpinGC_CisAis
int X_SpinGC_CisAis(long unsigned int j, struct BindStruct *X, long unsigned int is1_spin, long unsigned int sigma1)
Compute the grandcanonical spin state with bit mask is1_spin.
Definition: mltplySpinCore.c:183
SetDiagonalTEInterAll
int SetDiagonalTEInterAll(long unsigned int isite1, long unsigned int isite2, long unsigned int isigma1, long unsigned int isigma2, double dtmp_V, struct BindStruct *X, double complex *tmp_v0, double complex *tmp_v1)
Update the vector by the general two-body diagonal interaction, . (Using in Time Evolution mode)...
Definition: diagonalcalc.c:1431
DefineList::TEChemi
int ** TEChemi
Definition: struct.h:293
BindStruct::Large
struct LargeList Large
Variables for Matrix-Vector product.
Definition: struct.h:412
childfopenMPI
int childfopenMPI(const char *_cPathChild, const char *_cmode, FILE **_fp)
Only the root process open file in output/ directory.
Definition: FileIO.c:27
SetDiagonalInterAll
int SetDiagonalInterAll(long unsigned int isite1, long unsigned int isite2, long unsigned int isigma1, long unsigned int isigma2, double dtmp_V, struct BindStruct *X)
Calculate the components for general two-body diagonal interaction, .
Definition: diagonalcalc.c:1110
X_Spin_CisAis
int X_Spin_CisAis(long unsigned int j, struct BindStruct *X, long unsigned int is1_spin, long unsigned int sigma1)
Compute the spin state with bit mask is1_spin.
Definition: mltplySpinCore.c:166
cProStartCalcDiag
const char * cProStartCalcDiag
Definition: ProgressMessage.c:33
cFileNameCheckCoulombIntra
const char * cFileNameCheckCoulombIntra
Definition: global.c:27
DefineList::ParaTEChemi
double ** ParaTEChemi
Definition: struct.h:296
DefineList::NTEChemi
unsigned int * NTEChemi
Definition: struct.h:294
DefineList::InterAll_Diagonal
int ** InterAll_Diagonal
[DefineList::NinterAll_Diagonal][4] Interacted quartet
Definition: struct.h:162
DefineList::ParaTETransferDiagonal
double ** ParaTETransferDiagonal
Definition: struct.h:266
DefineList::HundCoupling
int ** HundCoupling
[DefineList::NHundCoupling][2] Index of Hund coupling. malloc in setmem_def().
Definition: struct.h:134
DefineList::NInterAll_Diagonal
unsigned int NInterAll_Diagonal
Number of interall term (diagonal)
Definition: struct.h:164
diagonalcalc
int diagonalcalc(struct BindStruct *X)
Calculate diagonal components and obtain the list, list_diagonal.
Definition: diagonalcalc.c:85
BindStruct
Bind.
Definition: struct.h:409
cProEndCalcDiag
const char * cProEndCalcDiag
Definition: ProgressMessage.c:34
DefineList::Nsite
unsigned int Nsite
Number of sites in the INTRA process region.
Definition: struct.h:56
DefineList::ParaCoulombInter
double * ParaCoulombInter
[DefineList::NCoulombInter]Coupling constant of off-site Coulomb interaction. malloc in setmem_def()...
Definition: struct.h:130
DefineList::Tpow
long unsigned int * Tpow
[2 * DefineList::NsiteMPI] malloc in setmem_def().
Definition: struct.h:90
DefineList::CoulombInter
int ** CoulombInter
Definition: struct.h:128
DefineList::NCoulombIntra
unsigned int NCoulombIntra
Definition: struct.h:121
SetDiagonalHund
int SetDiagonalHund(long unsigned int isite1, long unsigned int isite2, double dtmp_V, struct BindStruct *X)
Calculate the components for Hund interaction, .
Definition: diagonalcalc.c:785
DefineList::EDParaChemi
double * EDParaChemi
[DefineList::Nsite] On-site potential parameter. malloc in setmem_def().
Definition: struct.h:100
cFileNameCheckInterAll
const char * cFileNameCheckInterAll
Definition: global.c:31
SetDiagonalCoulombIntra
int SetDiagonalCoulombIntra(long unsigned int isite1, double dtmp_V, struct BindStruct *X)
Calculate the components for Coulombintra interaction, .
Definition: diagonalcalc.c:261
DefineList::iFlgGeneralSpin
int iFlgGeneralSpin
Flag for the general (Sz/=1/2) spin.
Definition: struct.h:86
FALSE
#define FALSE
Definition: global.h:25
list_1
long unsigned int * list_1
Definition: global.h:47
SetDiagonalTETransfer
int SetDiagonalTETransfer(long unsigned int isite1, double dtmp_V, long unsigned int spin, struct BindStruct *X, double complex *tmp_v0, double complex *tmp_v1)
Update the vector by the general one-body diagonal interaction, . (Using in Time Evolution mode)...
Definition: diagonalcalc.c:1940
cFileNameCheckInterU
const char * cFileNameCheckInterU
Definition: global.c:29
cErrNoModel
char * cErrNoModel
Error Message in diagonal calc.c.
Definition: ErrorMessage.c:99
DefineList::SiteToBit
long int * SiteToBit
[DefineList::NsiteMPI] Similar to DefineList::Tpow. For general spin.
Definition: struct.h:94
cFileNameTimeKeep
const char * cFileNameTimeKeep
Definition: global.c:23
DefineList::NTEInterAllDiagonal
unsigned int * NTEInterAllDiagonal
Definition: struct.h:276
DefineList::ParaTEInterAllDiagonal
double ** ParaTEInterAllDiagonal
Definition: struct.h:291
BitCheckGeneral
int BitCheckGeneral(const long unsigned int org_bit, const unsigned int org_isite, const unsigned int target_ispin, const long int *SiteToBit, const long unsigned int *Tpow)
bit check function for general spin
Definition: bitcalc.c:393
SetDiagonalCoulombInter
int SetDiagonalCoulombInter(long unsigned int isite1, long unsigned int isite2, double dtmp_V, struct BindStruct *X)
Calculate the components for Coulombinter interaction, .
Definition: diagonalcalc.c:544
DefineList::CoulombIntra
int ** CoulombIntra
Definition: struct.h:122
cDiagonalCalcFinish
const char * cDiagonalCalcFinish
Definition: LogMessage.c:57
DefineList::ParaHundCoupling
double * ParaHundCoupling
[DefineList::NHundCoupling] Hund coupling constant. malloc in setmem_def().
Definition: struct.h:136
X
struct EDMainCalStruct X
Definition: struct.h:432
list_Diagonal
double * list_Diagonal
Definition: global.h:46
DefineList::iCalcModel
int iCalcModel
Switch for model. 0:Hubbard, 1:Spin, 2:Kondo, 3:HubbardGC, 4:SpinGC, 5:KondoGC, 6:HubbardNConserved.
Definition: struct.h:198
diagonalcalcForTE
int diagonalcalcForTE(const int _istep, struct BindStruct *X, double complex *tmp_v0, double complex *tmp_v1)
Definition: diagonalcalc.c:197
SetDiagonalChemi
int SetDiagonalChemi(long unsigned int isite1, double dtmp_V, long unsigned int spin, struct BindStruct *X)
Calculate the components for the chemical potential .
Definition: diagonalcalc.c:374
DefineList::NCoulombInter
unsigned int NCoulombInter
Number of off-site Coulomb interaction.
Definition: struct.h:127
DefineList::ParaCoulombIntra
double * ParaCoulombIntra
Definition: struct.h:124
DefineList::EDNChemi
unsigned int EDNChemi
Number of on-site term.
Definition: struct.h:97
DefineList::TEInterAllDiagonal
int *** TEInterAllDiagonal
Definition: struct.h:284
myrank
int myrank
Process ID, defined in InitializeMPI()
Definition: global.h:163
DefineList::EDChemi
int * EDChemi
[DefineList::Nsite] Chemical potential. malloc in setmem_def().
Definition: struct.h:98
BindStruct::Check
struct CheckList Check
Size of the Hilbert space.
Definition: struct.h:411
DefineList::NTETransferDiagonal
unsigned int * NTETransferDiagonal
Definition: struct.h:258
TimeKeeper
int TimeKeeper(struct BindStruct *X, const char *cFileName, const char *cTimeKeeper_Message, const char *cWriteType)
Functions for writing a time log.
Definition: log.c:42
stdoutMPI
FILE * stdoutMPI
File pointer to the standard output defined in InitializeMPI()
Definition: global.h:165