HPhi++  3.1.0
StdFace_ModelUtil.cpp
Go to the documentation of this file.
1 /*
2 HPhi-mVMC-StdFace - Common input generator
3 Copyright (C) 2015 The University of Tokyo
4 
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14 
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
21 #include <cstdio>
22 #include <cstdlib>
23 #include <cmath>
24 #include <complex>
25 #include "StdFace_vals.hpp"
26 #include <cstring>
27 #ifdef __MPI
28 #include <mpi.h>
29 #endif
30 #include <iostream>
31 
36 void StdFace_exit(int errorcode
37 )
38 {
39  int ierr = 0;
40  fflush(stdout);
41  fflush(stderr);
42 #ifdef __MPI
43  fprintf(stdout, "\n\n ####### You DO NOT have to WORRY about the following MPI-ERROR MESSAGE. #######\n\n");
44  ierr = MPI_Abort(MPI_COMM_WORLD, errorcode);
45  ierr = MPI_Finalize();
46 #endif
47  if (ierr != 0) fprintf(stderr, "\n MPI_Finalize() = %d\n\n", ierr);
48  exit(errorcode);
49 }
56 void StdFace_trans(
57  struct StdIntList *StdI,
58  std::complex<double> trans0,
59  int isite,
60  int ispin,
61  int jsite,
62  int jspin
63 )
64 {
65  StdI->trans[StdI->ntrans] = trans0;
66  StdI->transindx[StdI->ntrans][0] = isite;
67  StdI->transindx[StdI->ntrans][1] = ispin;
68  StdI->transindx[StdI->ntrans][2] = jsite;
69  StdI->transindx[StdI->ntrans][3] = jspin;
70  StdI->ntrans = StdI->ntrans + 1;
71 }/*void StdFace_trans*/
76 void StdFace_Hopping(
77 struct StdIntList *StdI,
78  std::complex<double> trans0,
79  int isite,
80  int jsite,
81  double *dR
82 )
83 {
84  int ispin;
90 #if defined(_HPhi)
91  int it, ii;
92  double Cphase;
93  std::complex<double> coef;
94 
95  if (strcmp(StdI->method, "timeevolution") == 0 && StdI->PumpBody == 1) {
96  for (it = 0; it < StdI->Lanczos_max; it++) {
97  Cphase = 0.0f;
98  for (ii = 0; ii < 3; ii++) Cphase += /*2.0*StdI->pi */ StdI->At[it][ii] * dR[ii];
99  coef = std::complex<double>(cos(Cphase), sin(-Cphase));
100  for (ispin = 0; ispin < 2; ispin++) {
101  StdI->pump[it][StdI->npump[it]] = coef * trans0;
102  StdI->pumpindx[it][StdI->npump[it]][0] = isite;
103  StdI->pumpindx[it][StdI->npump[it]][1] = ispin;
104  StdI->pumpindx[it][StdI->npump[it]][2] = jsite;
105  StdI->pumpindx[it][StdI->npump[it]][3] = ispin;
106  StdI->npump[it] = StdI->npump[it] + 1;
107 
108  StdI->pump[it][StdI->npump[it]] = conj(coef * trans0);
109  StdI->pumpindx[it][StdI->npump[it]][0] = jsite;
110  StdI->pumpindx[it][StdI->npump[it]][1] = ispin;
111  StdI->pumpindx[it][StdI->npump[it]][2] = isite;
112  StdI->pumpindx[it][StdI->npump[it]][3] = ispin;
113  StdI->npump[it] = StdI->npump[it] + 1;
114  }/*for (ispin = 0; ispin < 2; ispin++)*/
115  }/*for (it = 0; it < StdI->Lanczos_max; it++)*/
116  }/*if (strcmp(StdI->method, "timeevolution") == 0)*/
117  else
118 #endif
119  for (ispin = 0; ispin < 2; ispin++) {
120  StdFace_trans(StdI, trans0, jsite, ispin, isite, ispin);
121  StdFace_trans(StdI, conj(trans0), isite, ispin, jsite, ispin);
122  }/*for (ispin = 0; ispin < 2; ispin++)*/
123 }/*void StdFace_Hopping*/
129 void StdFace_HubbardLocal(
130  struct StdIntList *StdI,
131  double mu0,
132  double h0,
133  double Gamma0,
134  double U0,
135  int isite
136 )
137 {
138  StdFace_trans(StdI, mu0 + 0.5 * h0, isite, 0, isite, 0);
139  StdFace_trans(StdI, mu0 - 0.5 * h0, isite, 1, isite, 1);
140  StdFace_trans(StdI, -0.5 * Gamma0, isite, 1, isite, 0);
141  StdFace_trans(StdI, -0.5 * Gamma0, isite, 0, isite, 1);
147  StdI->Cintra[StdI->NCintra] = U0;
148  StdI->CintraIndx[StdI->NCintra][0] = isite;
149  StdI->NCintra += 1;
150 }/*void StdFace_HubbardLocal*/
155 void StdFace_MagField(
156  struct StdIntList *StdI,
157  int S2,
158  double h,
159  double Gamma,
160  int isite
161 )
162 {
163  int ispin;
164  double S, Sz;
165 
166  S = (double)S2 * 0.5;
170  for (ispin = 0; ispin <= S2; ispin++){
177  Sz = (double)ispin - S;
178  StdFace_trans(StdI, -h * Sz, isite, ispin, isite, ispin);
189  if (ispin < S2) {
190  StdFace_trans(StdI, -0.5 * Gamma * sqrt(S*(S + 1.0) - Sz*(Sz + 1.0)),
191  isite, ispin + 1, isite, ispin);
192  StdFace_trans(StdI, -0.5 * Gamma * sqrt(S*(S + 1.0) - Sz*(Sz + 1.0)),
193  isite, ispin, isite, ispin + 1);
194  }/*if (ispin < S2)*/
195  }/*for (ispin = 0; ispin <= S2; ispin++)*/
196 }/*void StdFace_MagField*/
203 void StdFace_intr(
204  struct StdIntList *StdI,
205  std::complex<double> intr0,
206  int site1,
207  int spin1,
208  int site2,
209  int spin2,
210  int site3,
211  int spin3,
212  int site4,
213  int spin4
214 )
215 {
216  StdI->intr[StdI->nintr] = intr0;
217  StdI->intrindx[StdI->nintr][0] = site1; StdI->intrindx[StdI->nintr][1] = spin1;
218  StdI->intrindx[StdI->nintr][2] = site2; StdI->intrindx[StdI->nintr][3] = spin2;
219  StdI->intrindx[StdI->nintr][4] = site3; StdI->intrindx[StdI->nintr][5] = spin3;
220  StdI->intrindx[StdI->nintr][6] = site4; StdI->intrindx[StdI->nintr][7] = spin4;
221  StdI->nintr = StdI->nintr + 1;
222 }/*void StdFace_intr*/
227 void StdFace_GeneralJ(
228 struct StdIntList *StdI,
229  double J[3][3],
230  int Si2,
231  int Sj2,
232  int isite,
233  int jsite
234 )
235 {
236  int ispin, jspin, ZGeneral, ExGeneral;
237  double Si, Sj, Siz, Sjz;
238  std::complex<double> intr0;
242  ZGeneral = 1;
243  ExGeneral = 1;
244  if (Si2 == 1 || Sj2 == 1) {
251  ZGeneral = 0;
252 
253  StdI->Hund[StdI->NHund] = -0.5 * J[2][2];
254  StdI->HundIndx[StdI->NHund][0] = isite;
255  StdI->HundIndx[StdI->NHund][1] = jsite;
256  StdI->NHund += 1;
257 
258  StdI->Cinter[StdI->NCinter] = -0.25 * J[2][2];
259  StdI->CinterIndx[StdI->NCinter][0] = isite;
260  StdI->CinterIndx[StdI->NCinter][1] = jsite;
261  StdI->NCinter += 1;
262 
263  if (fabs(J[0][1]) < 0.000001 && fabs(J[1][0]) < 0.000001
264 #if defined(_mVMC)
265  && abs(J[0][0] - J[1][1]) < 0.000001
266 #endif
267  ) {
274  ExGeneral = 0;
275 
276 #if defined(_mVMC)
277  StdI->Ex[StdI->NEx] = - 0.25 * (J[0][0] + J[1][1]);
278 #else
279  if (strcmp(StdI->model, "kondo") == 0)
280  StdI->Ex[StdI->NEx] = -0.25 * (J[0][0] + J[1][1]);
281  else
282  StdI->Ex[StdI->NEx] = 0.25 * (J[0][0] + J[1][1]);
283 #endif
284  StdI->ExIndx[StdI->NEx][0] = isite;
285  StdI->ExIndx[StdI->NEx][1] = jsite;
286  StdI->NEx += 1;
287 
288  StdI->PairLift[StdI->NPairLift] = 0.25 * (J[0][0] - J[1][1]);
289  StdI->PLIndx[StdI->NPairLift][0] = isite;
290  StdI->PLIndx[StdI->NPairLift][1] = jsite;
291  StdI->NPairLift += 1;
292  }/*if (fabs(J[0][1]) < 0.000001 && fabs(J[1][0]) < 0.000001)*/
293  }
298  Si = 0.5 * (double)Si2;
299  Sj = 0.5 * (double)Sj2;
300 
301  for (ispin = 0; ispin <= Si2; ispin++) {
302  Siz = (double)ispin - Si;
303  for (jspin = 0; jspin <= Sj2; jspin++) {
304  Sjz = (double)jspin - Sj;
312  if (ZGeneral == 1) {
313  intr0 = J[2][2] * Siz * Sjz;
314  StdFace_intr(StdI, intr0,
315  isite, ispin, isite, ispin, jsite, jspin, jsite, jspin);
316  }
328  if ((ispin < Si2 && jspin < Sj2) && ExGeneral == 1) {
329  intr0 = 0.25 * std::complex<double>(J[0][0] + J[1][1], J[0][1] - J[1][0])
330  * sqrt(Si * (Si + 1.0) - Siz * (Siz + 1.0))
331  * sqrt(Sj * (Sj + 1.0) - Sjz * (Sjz + 1.0));
332  StdFace_intr(StdI, intr0,
333  isite, ispin + 1, isite, ispin, jsite, jspin, jsite, jspin + 1);
334  StdFace_intr(StdI, conj(intr0),
335  isite, ispin, isite, ispin + 1, jsite, jspin + 1, jsite, jspin);
336  }
348  if ((ispin < Si2 && jspin < Sj2) && ExGeneral == 1) {
349  intr0 = 0.5 * 0.5 * std::complex<double>(J[0][0] - J[1][1], - (J[0][1] + J[1][0]))
350  * sqrt(Si * (Si + 1.0) - Siz * (Siz + 1.0))
351  * sqrt(Sj * (Sj + 1.0) - Sjz * (Sjz + 1.0));
352  StdFace_intr(StdI, intr0,
353  isite, ispin + 1, isite, ispin, jsite, jspin + 1, jsite, jspin);
354  StdFace_intr(StdI, conj(intr0),
355  isite, ispin, isite, ispin + 1, jsite, jspin, jsite, jspin + 1);
356  }
367  if (ispin < Si2) {
368  intr0 = 0.5 * std::complex<double>(J[0][2], - J[1][2]) * sqrt(Si * (Si + 1.0) - Siz * (Siz + 1.0)) * Sjz;
369  StdFace_intr(StdI, intr0,
370  isite, ispin + 1, isite, ispin, jsite, jspin, jsite, jspin);
371  StdFace_intr(StdI, conj(intr0),
372  jsite, jspin, jsite, jspin, isite, ispin, isite, ispin + 1);
373  }/*if (ispin < Si2)*/
384  if (jspin < Sj2) {
385  intr0 = 0.5 * std::complex<double>(J[2][0], - J[2][1]) * Siz * sqrt(Sj * (Sj + 1.0) - Sjz * (Sjz + 1.0));
386  StdFace_intr(StdI, intr0,
387  isite, ispin, isite, ispin, jsite, jspin + 1, jsite, jspin);
388  StdFace_intr(StdI, conj(intr0),
389  jsite, jspin, jsite, jspin + 1, isite, ispin, isite, ispin);
390  }/*if (jspin < Sj2)*/
391  }/*for (jspin = 0; jspin <= Sj2; jspin++)*/
392  }/*for (ispin = 0; ispin <= Si2; ispin++)*/
393 }/*StdFace_GeneralJ*/
400 void StdFace_Coulomb(
401 struct StdIntList *StdI,
402  double V,
403  int isite,
404  int jsite
405 )
406 {
407  StdI->Cinter[StdI->NCinter] = V;
408  StdI->CinterIndx[StdI->NCinter][0] = isite;
409  StdI->CinterIndx[StdI->NCinter][1] = jsite;
410  StdI->NCinter += 1;
411 }/*void StdFace_Coulomb*/
418 void StdFace_PrintVal_d(
419  const char* valname,
420  double *val,
421  double val0
422 )
423 {
424  if (std::isnan(*val) == 1) {
425  *val = val0;
426  fprintf(stdout, " %15s = %-10.5f ###### DEFAULT VALUE IS USED ######\n", valname, *val);
427  }
428  else fprintf(stdout, " %15s = %-10.5f\n", valname, *val);
429 }/*void StdFace_PrintVal_d*/
436 void StdFace_PrintVal_dd(
437  char* valname,
438  double *val,
439  double val0,
440  double val1
441 )
442 {
443  if (std::isnan(*val) == 1) {
447  if (std::isnan(val0) == 1) *val = val1;
448  else *val = val0;
449  fprintf(stdout, " %15s = %-10.5f ###### DEFAULT VALUE IS USED ######\n", valname, *val);
450  }
451  else fprintf(stdout, " %15s = %-10.5f\n", valname, *val);
452 }/*void StdFace_PrintVal_dd*/
459 void StdFace_PrintVal_c(
460  const char* valname,
461  std::complex<double> *val,
462  std::complex<double> val0
463 )
464 {
465  if (std::isnan(real(*val)) == 1) {
466  *val = val0;
467  fprintf(stdout, " %15s = %-10.5f %-10.5f ###### DEFAULT VALUE IS USED ######\n", valname, real(*val), imag(*val));
468  }
469  else fprintf(stdout, " %15s = %-10.5f %-10.5f\n", valname, real(*val), imag(*val));
470 }/*void StdFace_PrintVal_c*/
477 void StdFace_PrintVal_i(
478  const char* valname,
479  int *val,
480  int val0
481 )
482 {
483  int NaN_i = 2147483647;/*The upper limt of Int*/
484 
485  if (*val == NaN_i) {
486  *val = val0;
487  fprintf(stdout, " %15s = %-10d ###### DEFAULT VALUE IS USED ######\n", valname, *val);
488  }
489  else fprintf(stdout, " %15s = %-10d\n", valname, *val);
490 }/*void StdFace_PrintVal_i*/
496 void StdFace_NotUsed_d(
497  const char* valname,
498  double val
499 )
500 {
501  if (std::isnan(val) == 0) {
502  fprintf(stdout, "\n Check ! %s is SPECIFIED but will NOT be USED. \n", valname);
503  fprintf(stdout, " Please COMMENT-OUT this line \n");
504  fprintf(stdout, " or check this input is REALLY APPROPRIATE for your purpose ! \n\n");
505  StdFace_exit(-1);
506  }
507 }/*void StdFace_NotUsed_d*/
513 void StdFace_NotUsed_c(
514  const char* valname,
515  std::complex<double> val
516 )
517 {
518  if (std::isnan(real(val)) == 0) {
519  fprintf(stdout, "\n Check ! %s is SPECIFIED but will NOT be USED. \n", valname);
520  fprintf(stdout, " Please COMMENT-OUT this line \n");
521  fprintf(stdout, " or check this input is REALLY APPROPRIATE for your purpose ! \n\n");
522  StdFace_exit(-1);
523  }
524 }/*void StdFace_NotUsed_c*/
530 void StdFace_NotUsed_J(
531  const char* valname,
532  double JAll,
533  double J[3][3]
534 )
535 {
536  int i1, i2;
537  char Jname[3][3][10];
538 
539  sprintf(Jname[0][0], "%sx", valname);
540  sprintf(Jname[0][1], "%sxy", valname);
541  sprintf(Jname[0][2], "%sxz", valname);
542  sprintf(Jname[1][0], "%syx", valname);
543  sprintf(Jname[1][1], "%sy", valname);
544  sprintf(Jname[1][2], "%syz", valname);
545  sprintf(Jname[2][0], "%szx", valname);
546  sprintf(Jname[2][1], "%szy", valname);
547  sprintf(Jname[2][2], "%sz", valname);
548 
549  StdFace_NotUsed_d(valname, JAll);
550 
551  for (i1 = 0; i1 < 3; i1++) {
552  for (i2 = 0; i2 < 3; i2++) {
553  StdFace_NotUsed_d(Jname[i1][i2], J[i1][i2]);
554  }/*for (j = 0; j < 3; j++)*/
555  }/*for (i = 0; i < 3; i++)*/
556 }/*void StdFace_NotUsed_J*/
562 void StdFace_NotUsed_i(
563  const char* valname,
564  int val
565 )
566 {
567  int NaN_i = 2147483647;
568 
569  if (val != NaN_i) {
570  fprintf(stdout, "\n Check ! %s is SPECIFIED but will NOT be USED. \n", valname);
571  fprintf(stdout, " Please COMMENT-OUT this line \n");
572  fprintf(stdout, " or check this input is REALLY APPROPRIATE for your purpose ! \n\n");
573  StdFace_exit(-1);
574  }
575 }/*void StdFace_NotUsed_i*/
581 void StdFace_RequiredVal_i(
582  const char* valname,
583  int val
584 )
585 {
586  int NaN_i = 2147483647;
587 
588  if (val == NaN_i){
589  fprintf(stdout, "ERROR ! %s is NOT specified !\n", valname);
590  StdFace_exit(-1);
591  }
592  else fprintf(stdout, " %15s = %-3d\n", valname, val);
593 }/*void StdFace_RequiredVal_i*/
599 static void StdFace_FoldSite(
600  struct StdIntList *StdI,
601  int iCellV[3],
602  int nBox[3],
603  int iCellV_fold[3]
605 )
606 {
607  int ii, jj, iCellV_frac[3];
611  for (ii = 0; ii < 3; ii++) {
612  iCellV_frac[ii] = 0;
613  for (jj = 0; jj < 3; jj++)iCellV_frac[ii] += StdI->rbox[ii][jj] * iCellV[jj];
614  }
618  for (ii = 0; ii < 3; ii++)
619  nBox[ii] = (iCellV_frac[ii] + StdI->NCell * 1000) / StdI->NCell - 1000;
623  for (ii = 0; ii < 3; ii++)
624  iCellV_frac[ii] -= StdI->NCell*(nBox[ii]);
625 
626  for (ii = 0; ii < 3; ii++) {
627  iCellV_fold[ii] = 0;
628  for (jj = 0; jj < 3; jj++) iCellV_fold[ii] += StdI->box[jj][ii] * iCellV_frac[jj];
629  iCellV_fold[ii] = (iCellV_fold[ii] + StdI->NCell * 1000) / StdI->NCell - 1000;
630  }/*for (ii = 0; ii < 3; ii++)*/
631 }/*static void StdFace_FoldSite*/
636 void StdFace_InitSite(
637  struct StdIntList *StdI,
638  FILE *fp,
639  int dim
640 )
641 {
642  int bound[3][2], edge, ii, jj;
643  int ipos;
644  int nBox[3], iCellV_fold[3], iCellV[3];
645  double pos[4][2], xmin, xmax/*, offset[2], scale*/;
646 
647  fprintf(stdout, "\n @ Super-Lattice setting\n\n");
651  if (
652  (StdI->L != StdI->NaN_i || StdI->W != StdI->NaN_i || StdI->Height != StdI->NaN_i)
653  &&
654  (StdI->box[0][0] != StdI->NaN_i || StdI->box[0][1] != StdI->NaN_i || StdI->box[0][2] != StdI->NaN_i ||
655  StdI->box[1][0] != StdI->NaN_i || StdI->box[1][1] != StdI->NaN_i || StdI->box[1][2] != StdI->NaN_i ||
656  StdI->box[2][0] != StdI->NaN_i || StdI->box[2][1] != StdI->NaN_i || StdI->box[2][2] != StdI->NaN_i)
657  )
658  {
659  fprintf(stdout, "\nERROR ! (L, W, Height) and (a0W, ..., a2H) conflict !\n\n");
660  StdFace_exit(-1);
661  }
662  else if (StdI->L != StdI->NaN_i || StdI->W != StdI->NaN_i || StdI->Height != StdI->NaN_i)
663  {
664  StdFace_PrintVal_i("L", &StdI->L, 1);
665  StdFace_PrintVal_i("W", &StdI->W, 1);
666  StdFace_PrintVal_i("Height", &StdI->Height, 1);
667  for (ii = 0; ii < 3; ii++) for (jj = 0; jj < 3; jj++)
668  StdI->box[ii][jj] = 0;
669  StdI->box[0][0] = StdI->W;
670  StdI->box[1][1] = StdI->L;
671  StdI->box[2][2] = StdI->Height;
672  }
673  else
674  {
675  StdFace_PrintVal_i("a0W", &StdI->box[0][0], 1);
676  StdFace_PrintVal_i("a0L", &StdI->box[0][1], 0);
677  StdFace_PrintVal_i("a0H", &StdI->box[0][2], 0);
678  StdFace_PrintVal_i("a1W", &StdI->box[1][0], 0);
679  StdFace_PrintVal_i("a1L", &StdI->box[1][1], 1);
680  StdFace_PrintVal_i("a1H", &StdI->box[1][2], 0);
681  StdFace_PrintVal_i("a2W", &StdI->box[2][0], 0);
682  StdFace_PrintVal_i("a2L", &StdI->box[2][1], 0);
683  StdFace_PrintVal_i("a2H", &StdI->box[2][2], 1);
684  }
685  /*
686  Parameters for the 3D system will not used.
687  */
688  if (dim == 2) {
689  StdI->direct[0][2] = 0.0;
690  StdI->direct[1][2] = 0.0;
691  StdI->direct[2][0] = 0.0;
692  StdI->direct[2][1] = 0.0;
693  StdI->direct[2][2] = 1.0;
694  }
699  if (dim == 2) StdI->phase[2] = 0.0;
700  for (ii = 0; ii < 3; ii++) {
701  StdI->ExpPhase[ii] = std::complex<double>(cos(StdI->pi180 * StdI->phase[ii]), sin(StdI->pi180 * StdI->phase[ii]));
702  if (abs(StdI->ExpPhase[ii] + 1.0) < 0.000001) StdI->AntiPeriod[ii] = 1;
703  else StdI->AntiPeriod[ii] = 0;
704  }
708  StdI->tau = (double **)malloc(sizeof(double*) * StdI->NsiteUC);
709  for (ii = 0; ii < StdI->NsiteUC; ii++) {
710  StdI->tau[ii] = (double *)malloc(sizeof(double) * 3);
711  }
716  StdI->NCell = 0;
717  for (ii = 0; ii < 3; ii++) {
718  StdI->NCell += StdI->box[0][ii]
719  * StdI->box[1][(ii + 1) % 3]
720  * StdI->box[2][(ii + 2) % 3]
721  - StdI->box[0][ii]
722  * StdI->box[1][(ii + 2) % 3]
723  * StdI->box[2][(ii + 1) % 3];
724  }
725  printf(" Number of Cell = %d\n", abs(StdI->NCell));
726  if (StdI->NCell == 0) {
727  StdFace_exit(-1);
728  }
729 
730  for (ii = 0; ii < 3; ii++) {
731  for (jj = 0; jj < 3; jj++) {
732  StdI->rbox[ii][jj] = StdI->box[(ii + 1) % 3][(jj + 1) % 3] * StdI->box[(ii + 2) % 3][(jj + 2) % 3]
733  - StdI->box[(ii + 1) % 3][(jj + 2) % 3] * StdI->box[(ii + 2) % 3][(jj + 1) % 3];
734  }
735  }
736  if (StdI->NCell < 0) {
737  for (ii = 0; ii < 3; ii++)
738  for (jj = 0; jj < 3; jj++)
739  StdI->rbox[ii][jj] *= -1;
740  StdI->NCell *= -1;
741  }/*if (StdI->NCell < 0)*/
746  for (ii = 0; ii < 3; ii++) {
747  bound[ii][0] = 0;
748  bound[ii][1] = 0;
749  for (nBox[2] = 0; nBox[2] < 2; nBox[2]++) {
750  for (nBox[1] = 0; nBox[1] < 2; nBox[1]++) {
751  for (nBox[0] = 0; nBox[0] < 2; nBox[0]++) {
752  edge = 0;
753  for (jj = 0; jj < 3; jj++) edge += nBox[jj] * StdI->box[jj][ii];
754  if (edge < bound[ii][0]) bound[ii][0] = edge;
755  if (edge > bound[ii][1]) bound[ii][1] = edge;
756  }
757  }
758  }
759  }
763  StdI->Cell = (int **)malloc(sizeof(int*) * StdI->NCell);
764  for (ii = 0; ii < StdI->NCell; ii++) {
765  StdI->Cell[ii] = (int *)malloc(sizeof(int) * 3);
766  }/*for (ii = 0; ii < StdI->NCell; ii++)*/
767  jj = 0;
768  for (iCellV[2] = bound[2][0]; iCellV[2] <= bound[2][1]; iCellV[2]++) {
769  for (iCellV[1] = bound[1][0]; iCellV[1] <= bound[1][1]; iCellV[1]++) {
770  for (iCellV[0] = bound[0][0]; iCellV[0] <= bound[0][1]; iCellV[0]++) {
771  StdFace_FoldSite(StdI, iCellV, nBox, iCellV_fold);
772  if (nBox[0] == 0 && nBox[1] == 0 && nBox[2] == 0) {
773  for (ii = 0; ii < 3; ii++)
774  StdI->Cell[jj][ii] = iCellV[ii];
775  jj += 1;
776  }/*if (lUC == 1)*/
777  }/*for (iCellV[0] = bound[0][0]; iCellV[0] <= bound[0][1]; iCellV[0]++*/
778  }/*for (iCellV[1] = bound[1][0]; iCellV[1] <= bound[1][1]; iCellV[1]++)*/
779  }/*for (iCellV[2] = bound[2][0]; iCellV[2] <= bound[2][1]; iCellV[2]++)*/
783  if (dim == 2) {
784  pos[0][0] = 0.0;
785  pos[0][1] = 0.0;
786  pos[1][0] = StdI->direct[0][0] * (double)StdI->box[0][0] + StdI->direct[1][0] * (double)StdI->box[0][1];
787  pos[1][1] = StdI->direct[0][1] * (double)StdI->box[0][0] + StdI->direct[1][1] * (double)StdI->box[0][1];
788  pos[2][0] = StdI->direct[0][0] * (double)StdI->box[1][0] + StdI->direct[1][0] * (double)StdI->box[1][1];
789  pos[2][1] = StdI->direct[0][1] * (double)StdI->box[1][0] + StdI->direct[1][1] * (double)StdI->box[1][1];
790  pos[3][0] = pos[1][0] + pos[2][0];
791  pos[3][1] = pos[1][1] + pos[2][1];
792 
793  xmin = 0.0;
794  xmax = 0.0;
795  for (ipos = 0; ipos < 4; ipos++) {
796  if (pos[ipos][0] < xmin) xmin = pos[ipos][0];
797  if (pos[ipos][0] > xmax) xmax = pos[ipos][0];
798  if (pos[ipos][1] < xmin) xmin = pos[ipos][1];
799  if (pos[ipos][1] > xmax) xmax = pos[ipos][1];
800  }
801  xmin -= 2.0;
802  xmax += 2.0;
803 
804  fprintf(fp, "#set terminal pdf color enhanced \\\n");
805  fprintf(fp, "#dashed dl 1.0 size 20.0cm, 20.0cm \n");
806  fprintf(fp, "#set output \"lattice.pdf\"\n");
807  fprintf(fp, "set xrange [%f: %f]\n", xmin, xmax);
808  fprintf(fp, "set yrange [%f: %f]\n", xmin, xmax);
809  fprintf(fp, "set size square\n");
810  fprintf(fp, "unset key\n");
811  fprintf(fp, "unset tics\n");
812  fprintf(fp, "unset border\n");
813 
814  fprintf(fp, "set style line 1 lc 1 lt 1\n");
815  fprintf(fp, "set style line 2 lc 5 lt 1\n");
816  fprintf(fp, "set style line 3 lc 0 lt 1\n");
817 
818  fprintf(fp, "set arrow from %f, %f to %f, %f nohead front ls 3\n", pos[0][0], pos[0][1], pos[1][0], pos[1][1]);
819  fprintf(fp, "set arrow from %f, %f to %f, %f nohead front ls 3\n", pos[1][0], pos[1][1], pos[3][0], pos[3][1]);
820  fprintf(fp, "set arrow from %f, %f to %f, %f nohead front ls 3\n", pos[3][0], pos[3][1], pos[2][0], pos[2][1]);
821  fprintf(fp, "set arrow from %f, %f to %f, %f nohead front ls 3\n", pos[2][0], pos[2][1], pos[0][0], pos[0][1]);
822  }/*if (dim == 2)*/
823 }/*void StdFace_InitSite2D*/
827 void StdFace_FindSite(
828  struct StdIntList *StdI,
829  int iW,
830  int iL,
831  int iH,
832  int diW,
833  int diL,
834  int diH,
835  int isiteUC,
836  int jsiteUC,
837  int *isite,
838  int *jsite,
839  std::complex<double> *Cphase,
840  double *dR
841 )
842 {
843  int iCell, jCell, kCell, ii;
844  int nBox[3], jCellV[3];
845 
846  dR[0] = - (double)diW + StdI->tau[isiteUC][0] - StdI->tau[jsiteUC][0];
847  dR[1] = - (double)diL + StdI->tau[isiteUC][1] - StdI->tau[jsiteUC][1];
848  dR[2] = - (double)diH + StdI->tau[isiteUC][2] - StdI->tau[jsiteUC][2];
849 
850  jCellV[0] = iW + diW;
851  jCellV[1] = iL + diL;
852  jCellV[2] = iH + diH;
853  StdFace_FoldSite(StdI, jCellV, nBox, jCellV);
854  *Cphase = 1.0;
855  for (ii = 0; ii < 3; ii++) *Cphase *= pow(StdI->ExpPhase[ii], (double)nBox[ii]);
856 
857  for (kCell = 0; kCell < StdI->NCell; kCell++) {
858  if (jCellV[0] == StdI->Cell[kCell][0] &&
859  jCellV[1] == StdI->Cell[kCell][1] &&
860  jCellV[2] == StdI->Cell[kCell][2])
861  {
862  jCell = kCell;
863  }
864  if (iW == StdI->Cell[kCell][0] &&
865  iL == StdI->Cell[kCell][1] &&
866  iH == StdI->Cell[kCell][2])
867  {
868  iCell = kCell;
869  }
870  }/*for (iCell = 0; iCell < StdI->NCell; iCell++)*/
871  *isite = iCell * StdI->NsiteUC + isiteUC;
872  *jsite = jCell * StdI->NsiteUC + jsiteUC;
873  if (strcmp(StdI->model, "kondo") == 0) {
874  *isite += StdI->NCell * StdI->NsiteUC;
875  *jsite += StdI->NCell * StdI->NsiteUC;
876  }
877 }/*void StdFace_FindSite*/
881 void StdFace_SetLabel(
882  struct StdIntList *StdI,
883  FILE *fp,
884  int iW,
885  int iL,
886  int diW,
887  int diL,
888  int isiteUC,
889  int jsiteUC,
890  int *isite,
891  int *jsite,
892  int connect,
893  std::complex<double> *Cphase,
894  double *dR
895 )
896 {
897  double xi, yi, xj, yj;
901  StdFace_FindSite(StdI, iW, iL, 0, -diW, -diL, 0, jsiteUC, isiteUC, isite, jsite, Cphase, dR);
902 
903  xi = StdI->direct[0][0] * ((double)iW + StdI->tau[jsiteUC][0])
904  + StdI->direct[1][0] * ((double)iL + StdI->tau[jsiteUC][1]);
905  yi = StdI->direct[0][1] * ((double)iW + StdI->tau[jsiteUC][0])
906  + StdI->direct[1][1] * ((double)iL + StdI->tau[jsiteUC][1]);
907 
908  xj = StdI->direct[0][0] * ((double)(iW - diW) + StdI->tau[isiteUC][0])
909  + StdI->direct[1][0] * ((double)(iL - diL) + StdI->tau[isiteUC][1]);
910  yj = StdI->direct[0][1] * ((double)(iW - diW) + StdI->tau[isiteUC][0])
911  + StdI->direct[1][1] * ((double)(iL - diL) + StdI->tau[isiteUC][1]);
912 
913  if (*isite < 10)fprintf(fp, "set label \"%1d\" at %f, %f center front\n", *isite, xi, yi);
914  else fprintf(fp, "set label \"%2d\" at %f, %f center front\n", *isite, xi, yi);
915  if (*jsite < 10)fprintf(fp, "set label \"%1d\" at %f, %f center front\n", *jsite, xj, yj);
916  else fprintf(fp, "set label \"%2d\" at %f, %f center front\n", *jsite, xj, yj);
917  if (connect < 3)
918  fprintf(fp, "set arrow from %f, %f to %f, %f nohead ls %d\n", xi, yi, xj, yj, connect);
922  StdFace_FindSite(StdI, iW, iL, 0, diW, diL, 0, isiteUC, jsiteUC, isite, jsite, Cphase, dR);
923 
924  xi = StdI->direct[1][0] * ((double)iL + StdI->tau[isiteUC][1])
925  + StdI->direct[0][0] * ((double)iW + StdI->tau[isiteUC][0]);
926  yi = StdI->direct[1][1] * ((double)iL + StdI->tau[isiteUC][1])
927  + StdI->direct[0][1] * ((double)iW + StdI->tau[isiteUC][0]);
928 
929  xj = StdI->direct[0][0] * ((double)(iW + diW) + StdI->tau[jsiteUC][0])
930  + StdI->direct[1][0] * ((double)(iL + diL) + StdI->tau[jsiteUC][1]);
931  yj = StdI->direct[0][1] * ((double)(iW + diW) + StdI->tau[jsiteUC][0])
932  + StdI->direct[1][1] * ((double)(iL + diL) + StdI->tau[jsiteUC][1]);
933 
934  if (*isite < 10)fprintf(fp, "set label \"%1d\" at %f, %f center front\n", *isite, xi, yi);
935  else fprintf(fp, "set label \"%2d\" at %f, %f center front\n", *isite, xi, yi);
936  if (*jsite < 10)fprintf(fp, "set label \"%1d\" at %f, %f center front\n", *jsite, xj, yj);
937  else fprintf(fp, "set label \"%2d\" at %f, %f center front\n", *jsite, xj, yj);
938  if (connect < 3)
939  fprintf(fp, "set arrow from %f, %f to %f, %f nohead ls %d\n", xi, yi, xj, yj, connect);
940 }/*void StdFace_SetLabel*/
944 void StdFace_PrintXSF(struct StdIntList *StdI) {
945  FILE *fp;
946  int ii, jj, kk, isite, iCell;
947  double vec[3];
948 
949  fp = fopen("lattice.xsf", "w");
950  fprintf(fp, "CRYSTAL\n");
951  fprintf(fp, "PRIMVEC\n");
952  for (ii = 0; ii < 3; ii++) {
953  for (jj = 0; jj < 3; jj++) {
954  vec[jj] = 0.0;
955  for (kk = 0; kk < 3; kk++)
956  vec[jj] += (double)StdI->box[ii][kk] * StdI->direct[kk][jj];
957  }
958  fprintf(fp, "%15.5f %15.5f %15.5f\n", vec[0], vec[1], vec[2]);
959  }
960  fprintf(fp, "PRIMCOORD\n");
961  fprintf(fp, "%d 1\n", StdI->NCell * StdI->NsiteUC);
962  for (iCell = 0; iCell < StdI->NCell; iCell++) {
963  for (isite = 0; isite < StdI->NsiteUC; isite++) {
964  for (jj = 0; jj < 3; jj++) {
965  vec[jj] = 0.0;
966  for (kk = 0; kk < 3; kk++)
967  vec[jj] += ((double)StdI->Cell[iCell][kk] + StdI->tau[isite][kk])
968  * StdI->direct[kk][jj];
969  }
970  fprintf(fp, "H %15.5f %15.5f %15.5f\n", vec[0], vec[1], vec[2]);
971  }
972  }
973  fflush(fp);
974  fclose(fp);
975 }/*void StdFace_PrintXSF*/
979 void StdFace_InputSpinNN(
980  double J[3][3],
981  double JAll,
982  double J0[3][3],
983  double J0All,
984  const char *J0name
985 )
986 {
987  int i1, i2, i3, i4;
988  char Jname[3][3][10];
989 
990  strcpy(Jname[0][0], "x\0");
991  strcpy(Jname[0][1], "xy\0");
992  strcpy(Jname[0][2], "xz\0");
993  strcpy(Jname[1][0], "yx\0");
994  strcpy(Jname[1][1], "y\0");
995  strcpy(Jname[1][2], "yz\0");
996  strcpy(Jname[2][0], "zx\0");
997  strcpy(Jname[2][1], "zy\0");
998  strcpy(Jname[2][2], "z\0");
999 
1000  if (std::isnan(JAll) == 0 && std::isnan(J0All) == 0) {
1001  fprintf(stdout, "\n ERROR! %s conflict !\n\n", J0name);
1002  StdFace_exit(-1);
1003  }
1004  for (i1 = 0; i1 < 3; i1++) {
1005  for (i2 = 0; i2 < 3; i2++) {
1006  if (std::isnan(JAll) == 0 && std::isnan(J[i1][i2]) == 0) {
1007  fprintf(stdout, "\n ERROR! J%s conflict !\n\n", Jname[i1][i2]);
1008  StdFace_exit(-1);
1009  }
1010  else if (std::isnan(J0All) == 0 && std::isnan(J[i1][i2]) == 0) {
1011  fprintf(stdout, "\n ERROR! %s and J%s conflict !\n\n",
1012  J0name, Jname[i1][i2]);
1013  StdFace_exit(-1);
1014  }
1015  else if (std::isnan(J0All) == 0 && std::isnan(J0[i1][i2]) == 0) {
1016  fprintf(stdout, "\n ERROR! %s and %s%s conflict !\n\n", J0name,
1017  J0name, Jname[i1][i2]);
1018  StdFace_exit(-1);
1019  }
1020  else if (std::isnan(J0[i1][i2]) == 0 && std::isnan(JAll) == 0) {
1021  fprintf(stdout, "\n ERROR! %s%s conflict !\n\n",
1022  J0name, Jname[i1][i2]);
1023  StdFace_exit(-1);
1024  }
1025  }/*for (j = 0; j < 3; j++)*/
1026  }/*for (i = 0; i < 3; i++)*/
1027 
1028  for (i1 = 0; i1 < 3; i1++) {
1029  for (i2 = 0; i2 < 3; i2++) {
1030  for (i3 = 0; i3 < 3; i3++) {
1031  for (i4 = 0; i4 < 3; i4++) {
1032  if (std::isnan(J0[i1][i2]) == 0 && std::isnan(J[i3][i4]) == 0) {
1033  fprintf(stdout, "\n ERROR! %s%s and J%s conflict !\n\n",
1034  J0name, Jname[i1][i2], Jname[i3][i4]);
1035  StdFace_exit(-1);
1036  }
1037  }/*for (i4 = 0; i4 < 3; i4++)*/
1038  }/*for (i3 = 0; i3 < 3; i3++)*/
1039  }/*for (j = 0; j < 3; j++)*/
1040  }/*for (i = 0; i < 3; i++)*/
1041 
1042  for (i1 = 0; i1 < 3; i1++) {
1043  for (i2 = 0; i2 < 3; i2++) {
1044  if (std::isnan(J0[i1][i2]) == 0)
1045  fprintf(stdout, " %14s%s = %-10.5f\n", J0name, Jname[i1][i2], J0[i1][i2]);
1046  else if (std::isnan(J[i1][i2]) == 0) {
1047  J0[i1][i2] = J[i1][i2];
1048  fprintf(stdout, " %14s%s = %-10.5f\n", J0name, Jname[i1][i2], J0[i1][i2]);
1049  }
1050  else if (i1 == i2 && std::isnan(J0All) == 0) {
1051  J0[i1][i2] = J0All;
1052  fprintf(stdout, " %14s%s = %-10.5f\n", J0name, Jname[i1][i2], J0[i1][i2]);
1053  }
1054  else if (i1 == i2 && std::isnan(JAll) == 0) {
1055  J0[i1][i2] = JAll;
1056  fprintf(stdout, " %14s%s = %-10.5f\n", J0name, Jname[i1][i2], J0[i1][i2]);
1057  }
1058  else {
1059  J0[i1][i2] = 0.0;
1060  }
1061  }/*for (i2 = 0; i2 < 3; i2++)*/
1062  }/*for (i = 0; i < 3; i++)*/
1063 }/*void StdFace_InputSpinNN*/
1067 void StdFace_InputSpin(
1068  double Jp[3][3],
1069  double JpAll,
1070  const char *Jpname
1071 )
1072 {
1073  int i1, i2;
1074  char Jname[3][3][10];
1075 
1076  strcpy(Jname[0][0], "x\0");
1077  strcpy(Jname[0][1], "xy\0");
1078  strcpy(Jname[0][2], "xz\0");
1079  strcpy(Jname[1][0], "yx\0");
1080  strcpy(Jname[1][1], "y\0");
1081  strcpy(Jname[1][2], "yz\0");
1082  strcpy(Jname[2][0], "zx\0");
1083  strcpy(Jname[2][1], "zy\0");
1084  strcpy(Jname[2][2], "z\0");
1085 
1086  for (i1 = 0; i1 < 3; i1++) {
1087  for (i2 = 0; i2 < 3; i2++) {
1088  if (std::isnan(JpAll) == 0 && std::isnan(Jp[i1][i2]) == 0) {
1089  fprintf(stdout, "\n ERROR! %s and %s%s conflict !\n\n", Jpname,
1090  Jpname, Jname[i1][i2]);
1091  StdFace_exit(-1);
1092  }
1093  }/*for (j = 0; j < 3; j++)*/
1094  }/*for (i = 0; i < 3; i++)*/
1095 
1096  for (i1 = 0; i1 < 3; i1++) {
1097  for (i2 = 0; i2 < 3; i2++) {
1098  if (std::isnan(Jp[i1][i2]) == 0)
1099  fprintf(stdout, " %14s%s = %-10.5f\n", Jpname, Jname[i1][i2], Jp[i1][i2]);
1100  else if (i1 == i2 && std::isnan(JpAll) == 0) {
1101  Jp[i1][i2] = JpAll;
1102  fprintf(stdout, " %14s%s = %-10.5f\n", Jpname, Jname[i1][i2], Jp[i1][i2]);
1103  }
1104  else {
1105  Jp[i1][i2] = 0.0;
1106  }
1107  }/*for (i2 = 0; i2 < 3; i2++)*/
1108  }/*for (i = 0; i < 3; i++)*/
1109 }/*void StdFace_InputSpin*/
1115 void StdFace_InputCoulombV(
1116  double V,
1117  double *V0,
1118  const char *V0name
1119 )
1120 {
1121  if (std::isnan(V) == 0 && std::isnan(*V0) == 0) {
1122  fprintf(stdout, "\n ERROR! %s conflicts !\n\n", V0name);
1123  StdFace_exit(-1);
1124  }
1125  else if (std::isnan(*V0) == 0)
1126  fprintf(stdout, " %15s = %-10.5f\n", V0name, *V0);
1127  else if (std::isnan(V) == 0) {
1128  *V0 = V;
1129  fprintf(stdout, " %15s = %-10.5f\n", V0name, *V0);
1130  }
1131  else {
1132  *V0 = 0.0;
1133  }
1134 }/*void StdFace_InputCoulombV*/
1140 void StdFace_InputHopp(
1141  std::complex<double> t,
1142  std::complex<double> *t0,
1143  const char *t0name
1144 )
1145 {
1146  if (std::isnan(real(t)) == 0 && std::isnan(real(*t0)) == 0) {
1147  fprintf(stdout, "\n ERROR! %s conflicts !\n\n", t0name);
1148  StdFace_exit(-1);
1149  }
1150  else if (std::isnan(real(*t0)) == 0)
1151  fprintf(stdout, " %15s = %-10.5f %-10.5f\n", t0name, real(*t0), imag(*t0));
1152  else if (std::isnan(real(t)) == 0) {
1153  *t0 = t;
1154  fprintf(stdout, " %15s = %-10.5f %-10.5f\n", t0name, real(*t0), imag(*t0));
1155  }
1156  else {
1157  *t0 = 0.0;
1158  }
1159 }/*void StdFace_InputHopp*/
1163 void StdFace_PrintGeometry(struct StdIntList *StdI) {
1164  FILE *fp;
1165  int isite, iCell, ii;
1166 
1167  fp = fopen("geometry.dat", "w");
1168 
1169  for (ii = 0; ii < 3; ii++)
1170  fprintf(fp, "%25.15e %25.15e %25.15e\n",
1171  StdI->direct[ii][0], StdI->direct[ii][1], StdI->direct[ii][2]);
1172  fprintf(fp, "%25.15e %25.15e %25.15e\n",
1173  StdI->phase[0], StdI->phase[1], StdI->phase[2]);
1174  for (ii = 0; ii < 3; ii++)
1175  fprintf(fp, "%d %d %d\n",
1176  StdI->box[ii][0], StdI->box[ii][1], StdI->box[ii][2]);
1177 
1178  for (iCell = 0; iCell < StdI->NCell; iCell++) {
1179  for (isite = 0; isite < StdI->NsiteUC; isite++) {
1180  fprintf(fp, "%d %d %d %d\n",
1181  StdI->Cell[iCell][0] - StdI->Cell[0][0],
1182  StdI->Cell[iCell][1] - StdI->Cell[0][1],
1183  StdI->Cell[iCell][2] - StdI->Cell[0][2],
1184  isite);
1185  }/*for (isite = 0; isite < StdI->NsiteUC; isite++)*/
1186  }/* for (iCell = 0; iCell < StdI->NCell; iCell++)*/
1187  if (strcmp(StdI->model, "kondo") == 0) {
1188  for (iCell = 0; iCell < StdI->NCell; iCell++) {
1189  for (isite = 0; isite < StdI->NsiteUC; isite++) {
1190  fprintf(fp, "%d %d %d %d\n",
1191  StdI->Cell[iCell][0] - StdI->Cell[0][0],
1192  StdI->Cell[iCell][1] - StdI->Cell[0][1],
1193  StdI->Cell[iCell][2] - StdI->Cell[0][2],
1194  isite + StdI->NsiteUC);
1195  }/*for (isite = 0; isite < StdI->NsiteUC; isite++)*/
1196  }/* for (iCell = 0; iCell < StdI->NCell; iCell++)*/
1197  }
1198  fflush(fp);
1199  fclose(fp);
1200 }/*void StdFace_PrintGeometry()*/
1204 void StdFace_MallocInteractions(
1205  struct StdIntList *StdI,
1206  int ntransMax,
1207  int nintrMax
1208 ) {
1209  int ii;
1210 #if defined(_HPhi)
1211  int it;
1212 #endif
1213 
1216  StdI->transindx = (int **)malloc(sizeof(int*) * ntransMax);
1217  StdI->trans = (std::complex<double> *)malloc(sizeof(std::complex<double>) * ntransMax);
1218  for (ii = 0; ii < ntransMax; ii++) {
1219  StdI->transindx[ii] = (int *)malloc(sizeof(int) * 4);
1220  }
1221  StdI->ntrans = 0;
1222 #if defined(_HPhi)
1223  if (strcmp(StdI->method, "timeevolution") == 0 && StdI->PumpBody == 1) {
1224  StdI->npump = (int *)malloc(sizeof(int) * StdI->Lanczos_max);
1225  StdI->pumpindx = (int ***)malloc(sizeof(int**) * StdI->Lanczos_max);
1226  StdI->pump = (std::complex<double> **)malloc(sizeof(std::complex<double>*) * StdI->Lanczos_max);
1227  for (it = 0; it < StdI->Lanczos_max; it++) {
1228  StdI->npump[it] = 0;
1229  StdI->pumpindx[it] = (int **)malloc(sizeof(int*) * ntransMax);
1230  StdI->pump[it] = (std::complex<double> *)malloc(sizeof(std::complex<double>) * ntransMax);
1231  for (ii = 0; ii < ntransMax; ii++) {
1232  StdI->pumpindx[it][ii] = (int *)malloc(sizeof(int) * 4);
1233  }
1234  }/*for (it = 0; it < StdI->Lanczos_max;)*/
1235  }/*if (strcmp(StdI->method, "timeevolution") == 0*/
1236 #endif
1237 
1240  StdI->intrindx = (int **)malloc(sizeof(int*) * nintrMax);
1241  StdI->intr = (std::complex<double> *)malloc(sizeof(std::complex<double>) * nintrMax);
1242  for (ii = 0; ii < nintrMax; ii++) {
1243  StdI->intrindx[ii] = (int *)malloc(sizeof(int) * 8);
1244  }
1245  StdI->nintr = 0;
1249  StdI->CintraIndx = (int **)malloc(sizeof(int*) * nintrMax);
1250  StdI->Cintra = (double *)malloc(sizeof(double) * nintrMax);
1251  for (ii = 0; ii < nintrMax; ii++) {
1252  StdI->CintraIndx[ii] = (int *)malloc(sizeof(int) * 1);
1253  }
1254  StdI->NCintra = 0;
1258  StdI->CinterIndx = (int **)malloc(sizeof(int*) * nintrMax);
1259  StdI->Cinter = (double *)malloc(sizeof(double) * nintrMax);
1260  for (ii = 0; ii < nintrMax; ii++) {
1261  StdI->CinterIndx[ii] = (int *)malloc(sizeof(int) * 2);
1262  }
1263  StdI->NCinter = 0;
1267  StdI->HundIndx = (int **)malloc(sizeof(int*) * nintrMax);
1268  StdI->Hund = (double *)malloc(sizeof(double) * nintrMax);
1269  for (ii = 0; ii < nintrMax; ii++) {
1270  StdI->HundIndx[ii] = (int *)malloc(sizeof(int) * 2);
1271  }
1272  StdI->NHund = 0;
1276  StdI->ExIndx = (int **)malloc(sizeof(int*) * nintrMax);
1277  StdI->Ex = (double *)malloc(sizeof(double) * nintrMax);
1278  for (ii = 0; ii < nintrMax; ii++) {
1279  StdI->ExIndx[ii] = (int *)malloc(sizeof(int) * 2);
1280  }
1281  StdI->NEx = 0;
1285  StdI->PLIndx = (int **)malloc(sizeof(int*) * nintrMax);
1286  StdI->PairLift = (double *)malloc(sizeof(double) * nintrMax);
1287  for (ii = 0; ii < nintrMax; ii++) {
1288  StdI->PLIndx[ii] = (int *)malloc(sizeof(int) * 2);
1289  }
1290  StdI->NPairLift = 0;
1294  StdI->PHIndx = (int **)malloc(sizeof(int*) * nintrMax);
1295  StdI->PairHopp = (double *)malloc(sizeof(double) * nintrMax);
1296  for (ii = 0; ii < nintrMax; ii++) {
1297  StdI->PHIndx[ii] = (int *)malloc(sizeof(int) * 2);
1298  }
1299  StdI->NPairHopp = 0;
1300 }/*void StdFace_MallocInteractions*/
1301 #if defined(_mVMC)
1302 
1306 static void StdFace_FoldSiteSub(
1307  struct StdIntList *StdI,
1308  int iCellV[3],
1309  int nBox[3],
1310  int iCellV_fold[3]
1311 )
1312 {
1313  int ii, jj, iCellV_frac[3];
1317  for (ii = 0; ii < 3; ii++) {
1318  iCellV_frac[ii] = 0;
1319  for (jj = 0; jj < 3; jj++)iCellV_frac[ii] += StdI->rboxsub[ii][jj] * iCellV[jj];
1320  }
1324  for (ii = 0; ii < 3; ii++)
1325  nBox[ii] = (iCellV_frac[ii] + StdI->NCellsub * 1000) / StdI->NCellsub - 1000;
1329  for (ii = 0; ii < 3; ii++)
1330  iCellV_frac[ii] -= StdI->NCellsub*(nBox[ii]);
1331 
1332  for (ii = 0; ii < 3; ii++) {
1333  iCellV_fold[ii] = 0;
1334  for (jj = 0; jj < 3; jj++) iCellV_fold[ii] += StdI->boxsub[jj][ii] * iCellV_frac[jj];
1335  iCellV_fold[ii] = (iCellV_fold[ii] + StdI->NCellsub * 1000) / StdI->NCellsub - 1000;
1336  }
1337 }/*static void StdFace_FoldSiteSub*/
1342 void StdFace_Proj(struct StdIntList *StdI)
1343 {
1344  FILE *fp;
1345  int jsite, iCell, jCell, kCell;
1346  int nBox[3], iCellV[3], jCellV[3], ii;
1347  int iSym;
1348  int **Sym, **Anti;
1349 
1350  Sym = (int **)malloc(sizeof(int*) * StdI->nsite);
1351  Anti = (int **)malloc(sizeof(int*) * StdI->nsite);
1352  for (jsite = 0; jsite < StdI->nsite; jsite++) {
1353  Sym[jsite] = (int *)malloc(sizeof(int) * StdI->nsite);
1354  Anti[jsite] = (int *)malloc(sizeof(int) * StdI->nsite);
1355  }
1359  StdI->NSym = 0;
1360  for (iCell = 0; iCell < StdI->NCell; iCell++) {
1361 
1362  StdFace_FoldSiteSub(StdI, StdI->Cell[iCell], nBox, iCellV);
1363 
1364  StdFace_FoldSite(StdI, iCellV, nBox, iCellV);
1365 
1366  if (iCellV[0] == StdI->Cell[iCell][0] &&
1367  iCellV[1] == StdI->Cell[iCell][1] &&
1368  iCellV[2] == StdI->Cell[iCell][2]) {
1372  for (jCell = 0; jCell < StdI->NCell; jCell++) {
1373 
1374  for (ii = 0; ii < 3; ii++)jCellV[ii] = StdI->Cell[jCell][ii] + iCellV[ii];
1375  StdFace_FoldSite(StdI, jCellV, nBox, jCellV);
1376 
1377  for (kCell = 0; kCell < StdI->NCell; kCell++) {
1378  if (jCellV[0] == StdI->Cell[kCell][0] &&
1379  jCellV[1] == StdI->Cell[kCell][1] &&
1380  jCellV[2] == StdI->Cell[kCell][2])
1381  {
1382 
1383  for (jsite = 0; jsite < StdI->NsiteUC; jsite++) {
1384 
1385  Sym[StdI->NSym][jCell*StdI->NsiteUC + jsite] = kCell*StdI->NsiteUC + jsite;
1386  Anti[StdI->NSym][jCell*StdI->NsiteUC + jsite]
1387  = StdI->AntiPeriod[0] * nBox[0]
1388  + StdI->AntiPeriod[1] * nBox[1]
1389  + StdI->AntiPeriod[2] * nBox[2];
1390 
1391  if (strcmp(StdI->model, "kondo") == 0) {
1392  Sym[StdI->NSym][StdI->nsite / 2 + jCell*StdI->NsiteUC + jsite] = StdI->nsite / 2 + kCell*StdI->NsiteUC + jsite;
1393  Anti[StdI->NSym][StdI->nsite / 2 + jCell*StdI->NsiteUC + jsite]
1394  = StdI->AntiPeriod[0] * nBox[0]
1395  + StdI->AntiPeriod[1] * nBox[1]
1396  + StdI->AntiPeriod[2] * nBox[2];
1397  }/*if (strcmp(StdI->model, "kondo") == 0)*/
1398 
1399  }/*for (jsite = 0; jsite < StdI->NsiteUC; jsite++)*/
1400 
1401  }/*if (jWfold == StdI->Cell[kCell][0] && jLfold == StdI->Cell[kCell][1])*/
1402  }/*for (kCell = 0; kCell < StdI->NCell; kCell++)*/
1403  }/*for (jCell = 0; jCell < StdI->NCell; jCell++)*/
1404  StdI->NSym += 1;
1405  }/*if (iWfold == iW && iLfold == iL)*/
1406  }/*for (iCell = 0; iCell < StdI->NCell; iCell++)*/
1407 
1408  fp = fopen("qptransidx.def", "w");
1409  fprintf(fp, "=============================================\n");
1410  fprintf(fp, "NQPTrans %10d\n", StdI->NSym);
1411  fprintf(fp, "=============================================\n");
1412  fprintf(fp, "======== TrIdx_TrWeight_and_TrIdx_i_xi ======\n");
1413  fprintf(fp, "=============================================\n");
1414  for (iSym = 0; iSym < StdI->NSym; iSym++) {
1415  fprintf(fp, "%d %10.5f\n", iSym, 1.0);
1416  }
1417  for (iSym = 0; iSym < StdI->NSym; iSym++) {
1418  for (jsite = 0; jsite < StdI->nsite; jsite++) {
1419  if (Anti[iSym][jsite] % 2 == 0) Anti[iSym][jsite] = 1;
1420  else Anti[iSym][jsite] = -1;
1421  if (StdI->AntiPeriod[0] == 1 || StdI->AntiPeriod[1] == 1 || StdI->AntiPeriod[2] == 1) {
1422  fprintf(fp, "%5d %5d %5d %5d\n", iSym, jsite, Sym[iSym][jsite], Anti[iSym][jsite]);
1423  }
1424  else {
1425  fprintf(fp, "%5d %5d %5d\n", iSym, jsite, Sym[iSym][jsite]);
1426  }
1427  }
1428  }
1429  fflush(fp);
1430  fclose(fp);
1431  fprintf(stdout, " qptransidx.def is written.\n");
1432 
1433  for (jsite = 0; jsite < StdI->nsite; jsite++) {
1434  free(Anti[jsite]);
1435  free(Sym[jsite]);
1436  }
1437  free(Sym);
1438  free(Anti);
1439 }/*void StdFace_Proj(struct StdIntList *StdI)*/
1444 static void StdFace_InitSiteSub(struct StdIntList *StdI)
1445 {
1446  int ii, jj, kk, prod;
1450  if ((StdI->Lsub != StdI->NaN_i || StdI->Wsub != StdI->NaN_i || StdI->Hsub != StdI->NaN_i)
1451  && (StdI->boxsub[0][0] != StdI->NaN_i || StdI->boxsub[0][1] != StdI->NaN_i || StdI->boxsub[0][2] != StdI->NaN_i ||
1452  StdI->boxsub[1][0] != StdI->NaN_i || StdI->boxsub[1][1] != StdI->NaN_i || StdI->boxsub[1][2] != StdI->NaN_i ||
1453  StdI->boxsub[2][0] != StdI->NaN_i || StdI->boxsub[2][1] != StdI->NaN_i || StdI->boxsub[2][2] != StdI->NaN_i))
1454  {
1455  fprintf(stdout, "\nERROR ! (Lsub, Wsub, Hsub) and (a0Wsub, ..., a2Hsub) conflict !\n\n");
1456  StdFace_exit(-1);
1457  }
1458  else if (StdI->Wsub != StdI->NaN_i || StdI->Lsub != StdI->NaN_i || StdI->Hsub != StdI->NaN_i) {
1459  StdFace_PrintVal_i("Lsub", &StdI->Lsub, 1);
1460  StdFace_PrintVal_i("Wsub", &StdI->Wsub, 1);
1461  StdFace_PrintVal_i("Hsub", &StdI->Hsub, 1);
1462  for (ii = 0; ii < 3; ii++) for (jj = 0; jj < 3; jj++)
1463  StdI->boxsub[ii][jj] = 0;
1464  StdI->boxsub[0][0] = StdI->Wsub;
1465  StdI->boxsub[1][1] = StdI->Lsub;
1466  StdI->boxsub[2][2] = StdI->Hsub;
1467  }
1468  else {
1469  StdFace_PrintVal_i("a0Wsub", &StdI->boxsub[0][0], StdI->box[0][0]);
1470  StdFace_PrintVal_i("a0Lsub", &StdI->boxsub[0][1], StdI->box[0][1]);
1471  StdFace_PrintVal_i("a0Hsub", &StdI->boxsub[0][2], StdI->box[0][2]);
1472  StdFace_PrintVal_i("a1Wsub", &StdI->boxsub[1][0], StdI->box[1][0]);
1473  StdFace_PrintVal_i("a1Lsub", &StdI->boxsub[1][1], StdI->box[1][1]);
1474  StdFace_PrintVal_i("a1Hsub", &StdI->boxsub[1][2], StdI->box[1][2]);
1475  StdFace_PrintVal_i("a2Wsub", &StdI->boxsub[2][0], StdI->box[2][0]);
1476  StdFace_PrintVal_i("a2Lsub", &StdI->boxsub[2][1], StdI->box[2][1]);
1477  StdFace_PrintVal_i("a2Hsub", &StdI->boxsub[2][2], StdI->box[2][2]);
1478  }
1482  StdI->NCellsub = 0;
1483  for (ii = 0; ii < 3; ii++) {
1484  StdI->NCellsub += StdI->boxsub[0][ii]
1485  * StdI->boxsub[1][(ii + 1) % 3]
1486  * StdI->boxsub[2][(ii + 2) % 3]
1487  - StdI->boxsub[0][ii]
1488  * StdI->boxsub[1][(ii + 2) % 3]
1489  * StdI->boxsub[2][(ii + 1) % 3];
1490  }
1491  printf(" Number of Cell in the sublattice: %d\n", abs(StdI->NCellsub));
1492  if (StdI->NCellsub == 0) {
1493  StdFace_exit(-1);
1494  }
1495 
1496  for (ii = 0; ii < 3; ii++) {
1497  for (jj = 0; jj < 3; jj++) {
1498  StdI->rboxsub[ii][jj] = StdI->boxsub[(ii + 1) % 3][(jj + 1) % 3] * StdI->boxsub[(ii + 2) % 3][(jj + 2) % 3]
1499  - StdI->boxsub[(ii + 1) % 3][(jj + 2) % 3] * StdI->boxsub[(ii + 2) % 3][(jj + 1) % 3];
1500  }
1501  }
1502  if (StdI->NCellsub < 0) {
1503  for (ii = 0; ii < 3; ii++)
1504  for (jj = 0; jj < 3; jj++)
1505  StdI->rboxsub[ii][jj] *= -1;
1506  StdI->NCellsub *= -1;
1507  }/*if (StdI->NCell < 0)*/
1511  for (ii = 0; ii < 3; ii++) {
1512  for (jj = 0; jj < 3; jj++) {
1513  prod = 0.0;
1514  for (kk = 0; kk < 3; kk++) prod += StdI->rboxsub[ii][kk] * (double)StdI->box[jj][kk];
1515  if (prod % StdI->NCellsub != 0) {
1516  printf("\n ERROR ! Sublattice is INCOMMENSURATE !\n\n");
1517  StdFace_exit(-1);
1518  }/*if (prod % StdI->NCellsub != 0)*/
1519  }
1520  }
1521 }/*void StdFace_InitSiteSub*/
1525 void StdFace_generate_orb(struct StdIntList *StdI) {
1526  int iCell, jCell, kCell, iCell2, jCell2, iOrb, isite, jsite, Anti;
1527  int nBox[3], iCellV[3], jCellV[3], dCellV[3], ii;
1528  int **CellDone;
1529 
1530  StdFace_InitSiteSub(StdI);
1531 
1532  StdI->Orb = (int **)malloc(sizeof(int*) * StdI->nsite);
1533  StdI->AntiOrb = (int **)malloc(sizeof(int*) * StdI->nsite);
1534  for (isite = 0; isite < StdI->nsite; isite++) {
1535  StdI->Orb[isite] = (int *)malloc(sizeof(int) * StdI->nsite);
1536  StdI->AntiOrb[isite] = (int *)malloc(sizeof(int) * StdI->nsite);
1537  }
1538  CellDone = (int **)malloc(sizeof(int*) * StdI->NCell);
1539  for (iCell = 0; iCell < StdI->NCell; iCell++) {
1540  CellDone[iCell] = (int *)malloc(sizeof(int) * StdI->NCell);
1541  for (jCell = 0; jCell < StdI->NCell; jCell++) {
1542  CellDone[iCell][jCell] = 0;
1543  }
1544  }
1545 
1546  iOrb = 0;
1547  for (iCell = 0; iCell < StdI->NCell; iCell++) {
1548 
1549  StdFace_FoldSiteSub(StdI, StdI->Cell[iCell], nBox, iCellV);
1550 
1551  StdFace_FoldSite(StdI, iCellV, nBox, iCellV);
1552 
1553  for (kCell = 0; kCell < StdI->NCell; kCell++) {
1554  if (iCellV[0] == StdI->Cell[kCell][0] &&
1555  iCellV[1] == StdI->Cell[kCell][1] &&
1556  iCellV[2] == StdI->Cell[kCell][2])
1557  {
1558  iCell2 = kCell;
1559  }
1560  }/*for (kCell = 0; kCell < StdI->NCell; kCell++)*/
1561 
1562  for (jCell = 0; jCell < StdI->NCell; jCell++) {
1563 
1564  for (ii = 0; ii < 3; ii++)
1565  jCellV[ii] = StdI->Cell[jCell][ii] + iCellV[ii] - StdI->Cell[iCell][ii];
1566 
1567  StdFace_FoldSite(StdI, jCellV, nBox, jCellV);
1568 
1569  for (kCell = 0; kCell < StdI->NCell; kCell++) {
1570  if (jCellV[0] == StdI->Cell[kCell][0] &&
1571  jCellV[1] == StdI->Cell[kCell][1] &&
1572  jCellV[2] == StdI->Cell[kCell][2])
1573  {
1574  jCell2 = kCell;
1575  }
1576  }/*for (kCell = 0; kCell < StdI->NCell; kCell++)*/
1577  /*
1578  AntiPeriodic factor
1579  */
1580  for (ii = 0; ii < 3; ii++)
1581  dCellV[ii] = StdI->Cell[jCell][ii] - StdI->Cell[iCell][ii];
1582  StdFace_FoldSite(StdI, dCellV, nBox, dCellV);
1583  Anti = 0;
1584  for (ii = 0; ii < 3; ii++)Anti += StdI->AntiPeriod[ii] * nBox[ii];
1585  if (Anti % 2 == 0) Anti = 1;
1586  else Anti = -1;
1587 
1588  for (isite = 0; isite < StdI->NsiteUC; isite++) {
1589  for (jsite = 0; jsite < StdI->NsiteUC; jsite++) {
1590 
1591  if (CellDone[iCell2][jCell2] == 0) {
1592  StdI->Orb[iCell2*StdI->NsiteUC + isite][jCell2*StdI->NsiteUC + jsite] = iOrb;
1593  StdI->AntiOrb[iCell2*StdI->NsiteUC + isite][jCell2*StdI->NsiteUC + jsite] = Anti;
1594  iOrb += 1;
1595  }
1596  StdI->Orb[iCell*StdI->NsiteUC + isite][jCell*StdI->NsiteUC + jsite]
1597  = StdI->Orb[iCell2*StdI->NsiteUC + isite][jCell2*StdI->NsiteUC + jsite];
1598  StdI->AntiOrb[iCell*StdI->NsiteUC + isite][jCell*StdI->NsiteUC + jsite] = Anti;
1599 
1600  if (strcmp(StdI->model, "kondo") == 0) {
1601  if (CellDone[iCell2][jCell2] == 0) {
1602  StdI->Orb[StdI->nsite / 2 + iCell2*StdI->NsiteUC + isite]
1603  [ jCell2*StdI->NsiteUC + jsite] = iOrb;
1604  StdI->AntiOrb[StdI->nsite / 2 + iCell2*StdI->NsiteUC + isite]
1605  [ jCell2*StdI->NsiteUC + jsite] = Anti;
1606  iOrb += 1;
1607  StdI->Orb[ iCell2*StdI->NsiteUC + isite]
1608  [StdI->nsite / 2 + jCell2*StdI->NsiteUC + jsite] = iOrb;
1609  StdI->AntiOrb[ iCell2*StdI->NsiteUC + isite]
1610  [StdI->nsite / 2 + jCell2*StdI->NsiteUC + jsite] = Anti;
1611  iOrb += 1;
1612  StdI->Orb[StdI->nsite / 2 + iCell2*StdI->NsiteUC + isite]
1613  [StdI->nsite / 2 + jCell2*StdI->NsiteUC + jsite] = iOrb;
1614  StdI->AntiOrb[StdI->nsite / 2 + iCell2*StdI->NsiteUC + isite]
1615  [StdI->nsite / 2 + jCell2*StdI->NsiteUC + jsite] = Anti;
1616  iOrb += 1;
1617  }
1618  StdI->Orb[StdI->nsite / 2 + iCell*StdI->NsiteUC + isite]
1619  [ jCell*StdI->NsiteUC + jsite]
1620  = StdI->Orb[StdI->nsite / 2 + iCell2*StdI->NsiteUC + isite]
1621  [ jCell2*StdI->NsiteUC + jsite];
1622  StdI->AntiOrb[StdI->nsite / 2 + iCell*StdI->NsiteUC + isite]
1623  [ jCell*StdI->NsiteUC + jsite] = Anti;
1624  StdI->Orb[ iCell*StdI->NsiteUC + isite]
1625  [StdI->nsite / 2 + jCell*StdI->NsiteUC + jsite]
1626  = StdI->Orb[ iCell2*StdI->NsiteUC + isite]
1627  [StdI->nsite / 2 + jCell2*StdI->NsiteUC + jsite];
1628  StdI->AntiOrb[iCell*StdI->NsiteUC + isite]
1629  [StdI->nsite / 2 + jCell*StdI->NsiteUC + jsite] = Anti;
1630  StdI->Orb[StdI->nsite / 2 + iCell*StdI->NsiteUC + isite]
1631  [StdI->nsite / 2 + jCell*StdI->NsiteUC + jsite]
1632  = StdI->Orb[StdI->nsite / 2 + iCell2*StdI->NsiteUC + isite]
1633  [StdI->nsite / 2 + jCell2*StdI->NsiteUC + jsite];
1634  StdI->AntiOrb[StdI->nsite / 2 + iCell*StdI->NsiteUC + isite]
1635  [StdI->nsite / 2 + jCell*StdI->NsiteUC + jsite] = Anti;
1636  }/*if (strcmp(StdI->model, "kondo") == 0)*/
1637 
1638  }/*for (jsite = 0; jsite < StdI->NsiteUC; jsite++)*/
1639  }/*for (isite = 0; isite < StdI->NsiteUC; isite++)*/
1640  CellDone[iCell2][jCell2] = 1;
1641  }/*for (jCell = 0; jCell < StdI->NCell; jCell++)*/
1642 
1643  }/*for (iCell = 0; iCell < StdI->NCell; iCell++)*/
1644  StdI->NOrb = iOrb;
1645 
1646  for (iCell = 0; iCell < StdI->NCell; iCell++) free(CellDone[iCell]);
1647  free(CellDone);
1648 }/*void StdFace_generate_orb*/
1653 void PrintJastrow(struct StdIntList *StdI) {
1654  FILE *fp;
1655  int isite, jsite, isiteUC, jsiteUC, revarsal, isite1, jsite1, iorb;
1656  int NJastrow, iJastrow;
1657  int dCell, iCell;//, jCell, dCellv[3];
1658  int **Jastrow;
1659  std::complex<double> Cphase;
1660  double dR[3];
1661 
1662  Jastrow = (int **)malloc(sizeof(int*) * StdI->nsite);
1663  for (isite = 0; isite < StdI->nsite; isite++)
1664  Jastrow[isite] = (int *)malloc(sizeof(int) * StdI->nsite);
1665 
1666  if (abs(StdI->NMPTrans) == 1 || StdI->NMPTrans == StdI->NaN_i) {
1670  for (isite = 0; isite < StdI->nsite; isite++) {
1671  for (jsite = 0; jsite < StdI->nsite; jsite++) {
1672  Jastrow[isite][jsite] = StdI->Orb[isite][jsite];
1673  }/*for (jsite = 0; jsite < isite; jsite++)*/
1674  }/*for (isite = 0; isite < StdI->nsite; isite++)*/
1678  for (iorb = 0; iorb < StdI->NOrb; iorb++) {
1679  for (isite = 0; isite < StdI->nsite; isite++) {
1680  for (jsite = 0; jsite < StdI->nsite; jsite++) {
1681  if (Jastrow[isite][jsite] == iorb) {
1682  Jastrow[jsite][isite] = Jastrow[isite][jsite];
1683  }
1684  }/*for (jsite = 0; jsite < isite; jsite++)*/
1685  }/*for (isite = 0; isite < StdI->nsite; isite++)*/
1686  }/*for (iorb = 0; iorb < StdI->NOrb; iorb++)*/
1687 
1688  if (strcmp(StdI->model, "hubbard") == 0) NJastrow = 0;
1689  else NJastrow = -1;
1690  for (isite = 0; isite < StdI->nsite; isite++) {
1691  /*
1692  For Local spin
1693  */
1694  if (StdI->locspinflag[isite] != 0) {
1695  for (jsite = 0; jsite < StdI->nsite; jsite++) {
1696  Jastrow[isite][jsite] = -1;
1697  Jastrow[jsite][isite] = -1;
1698  }
1699  continue;
1700  }
1701 
1702  for (jsite = 0; jsite < isite; jsite++) {
1703  if (Jastrow[isite][jsite] >= 0) {
1704  iJastrow = Jastrow[isite][jsite];
1705  NJastrow -= 1;
1706  for (isite1 = 0; isite1 < StdI->nsite; isite1++) {
1707  for (jsite1 = 0; jsite1 < StdI->nsite; jsite1++) {
1708  if (Jastrow[isite1][jsite1] == iJastrow)
1709  Jastrow[isite1][jsite1] = NJastrow;
1710  }/*for (jsite1 = 0; jsite1 < StdI->nsite; jsite1++)*/
1711  }/*for (isite1 = 0; isite1 < StdI->nsite; isite1++)*/
1712  }/*if (Jastrow[isite][jsite] >= 0)*/
1713  }/*for (jsite = 0; jsite < isite; jsite++)*/
1714  }/*for (isite = 0; isite < StdI->nsite; isite++)*/
1715 
1716  NJastrow = -NJastrow;
1717  for (isite = 0; isite < StdI->nsite; isite++) {
1718  for (jsite = 0; jsite < StdI->nsite; jsite++) {
1719  Jastrow[isite][jsite] = -1 - Jastrow[isite][jsite];
1720  }/*for (jsite = 0; jsite < isite; jsite++)*/
1721  }/*for (isite = 0; isite < StdI->nsite; isite++)*/
1722  }/*if (abs(StdI->NMPTrans) == 1)*/
1723  else {
1724 
1725  if (strcmp(StdI->model, "spin") == 0) {
1726  NJastrow = 1;
1727 
1728  for (isite = 0; isite < StdI->nsite; isite++) {
1729  for (jsite = 0; jsite < StdI->nsite; jsite++) {
1730  Jastrow[isite][jsite] = 0;
1731  }/*for (jsite = 0; jsite < StdI->nsite; jsite++)*/
1732  }/*for (isite = 0; isite < StdI->nsite; isite++)*/
1733  }/*if (strcmp(StdI->model, "spin") == 0)*/
1734  else {
1735 
1736  NJastrow = 0;
1737 
1738  if (strcmp(StdI->model, "kondo") == 0) {
1739  /*
1740  Local spin - itererant electron part
1741  */
1742  for (isite = 0; isite < StdI->nsite; isite++) {
1743  for (jsite = 0; jsite < StdI->nsite / 2; jsite++) {
1744  Jastrow[isite][jsite] = 0;
1745  Jastrow[jsite][isite] = 0;
1746  }/*for (jsite = 0; jsite < StdI->nsite; jsite++)*/
1747  }/*for (isite = 0; isite < StdI->nsite; isite++)*/
1748 
1749  NJastrow += 1;
1750  }/*if (strcmp(StdI->model, "kondo") == 0)*/
1751 
1752  for (dCell = 0; dCell < StdI->NCell; dCell++) {
1753  StdFace_FindSite(StdI,
1754  0, 0, 0,
1755  -StdI->Cell[dCell][0], -StdI->Cell[dCell][1], -StdI->Cell[dCell][2],
1756  0, 0, &isite, &jsite, &Cphase, dR);
1757  if (strcmp(StdI->model, "kondo") == 0) jsite += -StdI->NCell * StdI->NsiteUC;
1758  iCell = jsite / StdI->NsiteUC;
1759  if (iCell < dCell) {
1760  /*
1761  If -R has been already done, skip.
1762  */
1763  continue;
1764  }
1765  else if (iCell == dCell) {
1766  /*
1767  If revarsal symmetry [Fold(-R) = R], J(R,i,j) = J(R,j,i)
1768  */
1769  revarsal = 1;
1770  }
1771  else revarsal = 0;
1772 
1773  for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++) {
1774  for (jsiteUC = 0; jsiteUC < StdI->NsiteUC; jsiteUC++) {
1775  if (revarsal == 1 && jsiteUC > isiteUC) continue;/*If [Fold(-R) = R]*/
1776  if (isiteUC == jsiteUC &&
1777  StdI->Cell[dCell][0] == 0 &&
1778  StdI->Cell[dCell][1] == 0 &&
1779  StdI->Cell[dCell][2] == 0) continue;/*Diagonal*/
1780 
1781  for (iCell = 0; iCell < StdI->NCell; iCell++) {
1782  StdFace_FindSite(StdI,
1783  StdI->Cell[iCell][0], StdI->Cell[iCell][1], StdI->Cell[iCell][2],
1784  StdI->Cell[dCell][0], StdI->Cell[dCell][1], StdI->Cell[dCell][2],
1785  isiteUC, jsiteUC, &isite, &jsite, &Cphase, dR);
1786 
1787  Jastrow[isite][jsite] = NJastrow;
1788  Jastrow[jsite][isite] = NJastrow;
1789 
1790  }/*for (iCell = 0; iCell < StdI->NCell; iCell++)*/
1791 
1792  NJastrow += 1;
1793 
1794  }/*for (jsiteUC = 0; jsiteUC < StdI->NsiteUC; jsiteUC++)*/
1795  }/*for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++)*/
1796  }/*for (dCell = 0; dCell < StdI->NCell; dCell++)*/
1797  }/*if (strcmp(StdI->model, "spin") != 0)*/
1798  }/*if (abs(StdI->NMPTrans) != 1)*/
1799 
1800  fp = fopen("jastrowidx.def", "w");
1801  fprintf(fp, "=============================================\n");
1802  fprintf(fp, "NJastrowIdx %10d\n", NJastrow);
1803  fprintf(fp, "ComplexType %10d\n", 0);
1804  fprintf(fp, "=============================================\n");
1805  fprintf(fp, "=============================================\n");
1806 
1807  for (isite = 0; isite < StdI->nsite; isite++) {
1808  for (jsite = 0; jsite < StdI->nsite; jsite++) {
1809  if (isite == jsite) continue;
1810  fprintf(fp, "%5d %5d %5d\n", isite, jsite, Jastrow[isite][jsite]);
1811  }/*for (jsite = 0; jsite < isite; jsite++)*/
1812  }/*for (isite = 0; isite < StdI->nsite; isite++)*/
1813 
1814  for (iJastrow = 0; iJastrow < NJastrow; iJastrow++){
1815  if (strcmp(StdI->model, "hubbard") == 0 || iJastrow > 0)
1816  fprintf(fp, "%5d %5d\n", iJastrow, 1);
1817  else
1818  fprintf(fp, "%5d %5d\n", iJastrow, 0);
1819  }/*for (iJastrow = 0; iJastrow < NJastrow; iJastrow++)*/
1820  fflush(fp);
1821  fclose(fp);
1822  fprintf(stdout, " jastrowidx.def is written.\n");
1823 
1824  for (isite = 0; isite < StdI->nsite; isite++) free(Jastrow[isite]);
1825  free(Jastrow);
1826 }/*void PrintJastrow*/
1827 #endif
StdFace_InputHopp
void StdFace_InputHopp(std::complex< double > t, std::complex< double > *t0, const char *t0name)
Input hopping integral from the input file, if it is not specified, use the default value(0 or the is...
Definition: StdFace_ModelUtil.cpp:1140
StdFace_SetLabel
void StdFace_SetLabel(struct StdIntList *StdI, FILE *fp, int iW, int iL, int diW, int diL, int isiteUC, int jsiteUC, int *isite, int *jsite, int connect, std::complex< double > *Cphase, double *dR)
Set Label in the gnuplot display (Only used in 2D system)
Definition: StdFace_ModelUtil.cpp:881
StdFace_Coulomb
void StdFace_Coulomb(struct StdIntList *StdI, double V, int isite, int jsite)
Add onsite/offsite Coulomb term to the list StdIntList::Cinter and StdIntList::CinterIndx, and increase the number of them (StdIntList::NCinter).
Definition: StdFace_ModelUtil.cpp:400
StdFace_GeneralJ
void StdFace_GeneralJ(struct StdIntList *StdI, double J[3][3], int Si2, int Sj2, int isite, int jsite)
Treat J as a 3*3 matrix [(6S + 1)*(6S&#39; + 1) interactions].
Definition: StdFace_ModelUtil.cpp:227
StdFace_InputSpin
void StdFace_InputSpin(double Jp[3][3], double JpAll, const char *Jpname)
Input spin-spin interaction other than nearest-neighbor.
Definition: StdFace_ModelUtil.cpp:1067
StdFace_PrintVal_d
void StdFace_PrintVal_d(const char *valname, double *val, double val0)
Print a valiable (real) read from the input file if it is not specified in the input file (=NaN)...
Definition: StdFace_ModelUtil.cpp:418
StdFace_intr
void StdFace_intr(struct StdIntList *StdI, std::complex< double > intr0, int site1, int spin1, int site2, int spin2, int site3, int spin3, int site4, int spin4)
Add interaction (InterAll) to the list Set StdIntList::intr and StdIntList::intrindx and increase the...
Definition: StdFace_ModelUtil.cpp:203
StdFace_RequiredVal_i
void StdFace_RequiredVal_i(const char *valname, int val)
Stop HPhi if a variable (integer) which must be specified is absent in the input file (=2147483647...
Definition: StdFace_ModelUtil.cpp:581
StdFace_FindSite
void StdFace_FindSite(struct StdIntList *StdI, int iW, int iL, int iH, int diW, int diL, int diH, int isiteUC, int jsiteUC, int *isite, int *jsite, std::complex< double > *Cphase, double *dR)
Find the index of transfer and interaction.
Definition: StdFace_ModelUtil.cpp:827
StdFace_exit
void StdFace_exit(int errorcode)
MPI Abortation wrapper.
Definition: StdFace_ModelUtil.cpp:36
StdFace_Hopping
void StdFace_Hopping(struct StdIntList *StdI, std::complex< double > trans0, int isite, int jsite, double *dR)
Add Hopping for the both spin.
Definition: StdFace_ModelUtil.cpp:76
StdFace_NotUsed_c
void StdFace_NotUsed_c(const char *valname, std::complex< double > val)
Stop HPhi if a variable (complex) not used is specified in the input file (!=NaN).
Definition: StdFace_ModelUtil.cpp:513
StdFace_HubbardLocal
void StdFace_HubbardLocal(struct StdIntList *StdI, double mu0, double h0, double Gamma0, double U0, int isite)
Add intra-Coulomb, magnetic field, chemical potential for the itenerant electron. ...
Definition: StdFace_ModelUtil.cpp:129
StdFace_NotUsed_J
void StdFace_NotUsed_J(const char *valname, double JAll, double J[3][3])
Stop HPhi if variables (real) not used is specified in the input file (!=NaN).
Definition: StdFace_ModelUtil.cpp:530
StdFace_FoldSite
static void StdFace_FoldSite(struct StdIntList *StdI, int iCellV[3], int nBox[3], int iCellV_fold[3])
Move a site into the original supercell if it is outside the original supercell.
Definition: StdFace_ModelUtil.cpp:599
StdFace_MagField
void StdFace_MagField(struct StdIntList *StdI, int S2, double h, double Gamma, int isite)
Add longitudinal and transvars magnetic field to the list.
Definition: StdFace_ModelUtil.cpp:155
StdFace_PrintGeometry
void StdFace_PrintGeometry(struct StdIntList *StdI)
Print geometry of sites for the pos-process of correlation function.
Definition: StdFace_ModelUtil.cpp:1163
StdFace_PrintXSF
void StdFace_PrintXSF(struct StdIntList *StdI)
Print lattice.xsf (XCrysDen format)
Definition: StdFace_ModelUtil.cpp:944
StdFace_InputCoulombV
void StdFace_InputCoulombV(double V, double *V0, const char *V0name)
Input off-site Coulomb interaction from the input file, if it is not specified, use the default value...
Definition: StdFace_ModelUtil.cpp:1115
StdFace_trans
void StdFace_trans(struct StdIntList *StdI, std::complex< double > trans0, int isite, int ispin, int jsite, int jspin)
Add transfer to the list set StdIntList::trans and StdIntList::transindx and increment StdIntList::nt...
Definition: StdFace_ModelUtil.cpp:56
StdFace_MallocInteractions
void StdFace_MallocInteractions(struct StdIntList *StdI, int ntransMax, int nintrMax)
Malloc Arrays for interactions.
Definition: StdFace_ModelUtil.cpp:1204
StdFace_NotUsed_i
void StdFace_NotUsed_i(const char *valname, int val)
Stop HPhi if a variable (integer) not used is specified in the input file (!=2147483647, the upper limt of Int).
Definition: StdFace_ModelUtil.cpp:562
StdFace_PrintVal_i
void StdFace_PrintVal_i(const char *valname, int *val, int val0)
Print a valiable (integer) read from the input file if it is not specified in the input file (=214748...
Definition: StdFace_ModelUtil.cpp:477
StdFace_PrintVal_dd
void StdFace_PrintVal_dd(char *valname, double *val, double val0, double val1)
Print a valiable (real) read from the input file if it is not specified in the input file (=NaN)...
Definition: StdFace_ModelUtil.cpp:436
StdFace_InitSite
void StdFace_InitSite(struct StdIntList *StdI, FILE *fp, int dim)
Initialize the super-cell where simulation is performed.
Definition: StdFace_ModelUtil.cpp:636
StdFace_PrintVal_c
void StdFace_PrintVal_c(const char *valname, std::complex< double > *val, std::complex< double > val0)
Print a valiable (complex) read from the input file if it is not specified in the input file (=NaN)...
Definition: StdFace_ModelUtil.cpp:459
StdFace_InputSpinNN
void StdFace_InputSpinNN(double J[3][3], double JAll, double J0[3][3], double J0All, const char *J0name)
Input nearest-neighbor spin-spin interaction.
Definition: StdFace_ModelUtil.cpp:979
StdFace_NotUsed_d
void StdFace_NotUsed_d(const char *valname, double val)
Stop HPhi if a variable (real) not used is specified in the input file (!=NaN).
Definition: StdFace_ModelUtil.cpp:496