HPhi++  3.1.0
Functions
Wannier90.cpp File Reference

Standard mode for wannier90. More...

#include "StdFace_vals.hpp"
#include "StdFace_ModelUtil.hpp"
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <complex>
#include <cstring>

Go to the source code of this file.

Functions

static void geometry_W90 (struct StdIntList *StdI)
 Read Geometry file for wannier90. More...
 
static void read_W90 (struct StdIntList *StdI, char *filename, double cutoff, int *cutoff_R, double cutoff_length, int itUJ, int *NtUJ, int ***tUJindx, std::complex< double > **tUJ)
 Read Wannier90 hamiltonian file (*_hr) More...
 
static std::complex< double > ***** read_density_matrix (struct StdIntList *StdI, char *filename)
 Read RESPACK Density-matrix file (*_dr.dat) More...
 
static void PrintUHFinitial (struct StdIntList *StdI, int *NtUJ, std::complex< double > *****DenMat, int ***tUJindx)
 Print the initial guess of UHF. More...
 
void StdFace_Wannier90 (struct StdIntList *StdI)
 Setup a Hamiltonian for the Wannier90 *_hr.dat. More...
 

Detailed Description

Standard mode for wannier90.

Definition in file Wannier90.cpp.

Function Documentation

◆ geometry_W90()

static void geometry_W90 ( struct StdIntList *  StdI)
static

Read Geometry file for wannier90.

Author
Mitsuaki Kawamura (The University of Tokyo)

Direct lattice vector StdIntList::direct

Intrinsic site position StdIntList::tau and its number StdIntList::NsiteUC

Parameters
[in,out]StdI

Definition at line 32 of file Wannier90.cpp.

Referenced by StdFace_Wannier90().

35 {
36  int isite, ii, ierr;
37  char filename[265];
38  FILE *fp;
39 
40  sprintf(filename, "%s_geom.dat", StdI->CDataFileHead);
41  fprintf(stdout, " Wannier90 Geometry file = %s\n", filename);
42 
43  fp = fopen(filename, "r");
47  for (ii = 0; ii < 3; ii++)
48  ierr = fscanf(fp, "%lf%lf%lf", &StdI->direct[ii][0], &StdI->direct[ii][1], &StdI->direct[ii][2]);
49  if(ierr == EOF) printf("%d\n", ierr);
53  for (isite = 0; isite < StdI->NsiteUC; isite++) free(StdI->tau[isite]);
54  free(StdI->tau);
55  ierr = fscanf(fp, "%d", &StdI->NsiteUC);
56  fprintf(stdout, " Number of Correlated Sites = %d\n", StdI->NsiteUC);
57 
58  StdI->tau = (double **)malloc(sizeof(double*) * StdI->NsiteUC);
59  for (ii = 0; ii < StdI->NsiteUC; ii++) StdI->tau[ii] = (double *)malloc(sizeof(double) * 3);
60 
61  for (isite = 0; isite < StdI->NsiteUC; isite++)
62  ierr = fscanf(fp, "%lf%lf%lf", &StdI->tau[isite][0], &StdI->tau[isite][1], &StdI->tau[isite][2]);
63  fclose(fp);
64 
65  printf(" Direct lattice vectors:\n");
66  for (ii = 0; ii < 3; ii++) printf(" %10.5f %10.5f %10.5f\n",
67  StdI->direct[ii][0], StdI->direct[ii][1], StdI->direct[ii][2]);
68  printf(" Wannier centres:\n");
69  for (isite = 0; isite < StdI->NsiteUC; isite++) printf(" %10.5f %10.5f %10.5f\n",
70  StdI->tau[isite][0], StdI->tau[isite][1], StdI->tau[isite][2]);
71 
72 }/*static void geometry_W90(struct StdIntList *StdI) */

◆ PrintUHFinitial()

static void PrintUHFinitial ( struct StdIntList *  StdI,
int *  NtUJ,
std::complex< double > *****  DenMat,
int ***  tUJindx 
)
static

Print the initial guess of UHF.

Author
Mitsuaki Kawamura (The University of Tokyo)

Definition at line 333 of file Wannier90.cpp.

References StdFace_FindSite().

Referenced by StdFace_Wannier90().

338  {
339  FILE *fp;
340  int iW, iL, iH, kCell, it, isite, jsite ,NIniGuess, ispin;
341  double dR[3];
342  std::complex<double> Cphase, **IniGuess;
343 
344  IniGuess = (std::complex<double> **)malloc(sizeof(std::complex<double>*) * StdI->nsite);
345  for (isite = 0; isite < StdI->nsite; isite++)
346  IniGuess[isite] = (std::complex<double> *)malloc(sizeof(std::complex<double>) * StdI->nsite);
347 
348  for (kCell = 0; kCell < StdI->NCell; kCell++) {
349 
350  iW = StdI->Cell[kCell][0];
351  iL = StdI->Cell[kCell][1];
352  iH = StdI->Cell[kCell][2];
353  /*
354  Diagonal term
355  */
356  for (isite = 0; isite < StdI->NsiteUC; isite++) {
357  jsite = isite + StdI->NsiteUC*kCell;
358  IniGuess[jsite][jsite] = DenMat[0][0][0][isite][isite];
359  }
360  /*
361  Coulomb integral (U)
362  */
363  for (it = 0; it < NtUJ[1]; it++) {
364  StdFace_FindSite(StdI, iW, iL, iH,
365  tUJindx[1][it][0], tUJindx[1][it][1], tUJindx[1][it][2],
366  tUJindx[1][it][3], tUJindx[1][it][4], &isite, &jsite, &Cphase, dR);
367  IniGuess[isite][jsite] = DenMat[tUJindx[1][it][0]][tUJindx[1][it][1]][tUJindx[1][it][2]]
368  [tUJindx[1][it][3]][tUJindx[1][it][4]];
369  IniGuess[jsite][isite] = conj(DenMat[tUJindx[1][it][0]][tUJindx[1][it][1]][tUJindx[1][it][2]]
370  [tUJindx[1][it][3]][tUJindx[1][it][4]]);
371  }/*for (it = 0; it < NtUJ[1]; it++)*/
372  /*
373  Wxchange integral (J)
374  */
375  for (it = 0; it < NtUJ[2]; it++) {
376  StdFace_FindSite(StdI, iW, iL, iH,
377  tUJindx[2][it][0], tUJindx[2][it][1], tUJindx[2][it][2],
378  tUJindx[2][it][3], tUJindx[2][it][4], &isite, &jsite, &Cphase, dR);
379  IniGuess[isite][jsite] = DenMat[tUJindx[2][it][0]][tUJindx[2][it][1]][tUJindx[2][it][2]]
380  [tUJindx[2][it][3]][tUJindx[2][it][4]];
381  IniGuess[jsite][isite] = conj(DenMat[tUJindx[2][it][0]][tUJindx[2][it][1]][tUJindx[2][it][2]]
382  [tUJindx[2][it][3]][tUJindx[2][it][4]]);
383  }/*for (it = 0; it < NtUJ[1]; it++)*/
384  }/*for (kCell = 0; kCell < StdI->NCell; kCell++)*/
385 
386  NIniGuess = 0;
387  for (isite = 0; isite < StdI->nsite; isite++) {
388  for (jsite = 0; jsite < StdI->nsite; jsite++) {
389  if (abs(IniGuess[isite][jsite]) > 1.0e-6)NIniGuess += 1;
390  }/*for (jsite = 0; jsite < StdI->nsite; jsite++)*/
391  }/*for (isite = 0; isite < StdI->nsite; isite++)*/
392 
393  fp = fopen("initial.def", "w");
394  fprintf(fp, "======================== \n");
395  fprintf(fp, "NInitialGuess %7d \n", NIniGuess*2);
396  fprintf(fp, "======================== \n");
397  fprintf(fp, "========i_j_s_tijs====== \n");
398  fprintf(fp, "======================== \n");
399 
400  for (isite = 0; isite < StdI->nsite; isite++) {
401  for (jsite = 0; jsite < StdI->nsite; jsite++) {
402  if (abs(IniGuess[isite][jsite]) > 1.0e-6)
403  for (ispin = 0; ispin < 2; ispin++) {
404  fprintf(fp, "%5d %5d %5d %5d %25.15f %25.15f\n",
405  jsite, ispin, isite, ispin,
406  0.5*real(IniGuess[isite][jsite]),
407  0.5*imag(IniGuess[isite][jsite]));
408  }
409  }/*for (jsite = 0; jsite < StdI->nsite; jsite++)*/
410  }/*for (isite = 0; isite < StdI->nsite; isite++)*/
411 
412  fflush(fp);
413  fclose(fp);
414  fprintf(stdout, " initial.def is written.\n");
415 
416  for (isite = 0; isite < StdI->nsite; isite++)
417  free(IniGuess[isite]);
418  free(IniGuess);
419 }/*static void PrintTrans*/
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

◆ read_density_matrix()

static std::complex<double>***** read_density_matrix ( struct StdIntList *  StdI,
char *  filename 
)
static

Read RESPACK Density-matrix file (*_dr.dat)

Author
Mitsuaki Kawamura (The University of Tokyo)
Parameters
[in,out]StdI
[in]filenameInput file name

Definition at line 231 of file Wannier90.cpp.

Referenced by StdFace_Wannier90().

235 {
236  FILE *fp;
237  int ierr, nWan, nWSC, iWSC, iWan, jWan, iWan0, jWan0, ii, Rmax[3], Rmin[3], NR[3], i0, i1, i2;
238  double dtmp[2];
239  char ctmp[256], *ctmp2;
240  std::complex<double> ***Mat_tot, *****DenMat;
241  int **indx_tot;
242  /*
243  Header part
244  */
245  fp = fopen(filename, "r");
246  ctmp2 = fgets(ctmp, 256, fp);
247  ierr = fscanf(fp, "%d", &nWan);
248  if (ierr == EOF) printf("%d %s\n", ierr, ctmp2);
249  ierr = fscanf(fp, "%d", &nWSC);
250  for (iWSC = 0; iWSC < nWSC; iWSC++) {
251  ierr = fscanf(fp, "%d", &ii);
252  }
253  /*
254  Malloc Matrix elements and their indices
255  */
256  Mat_tot = (std::complex<double> ***)malloc(sizeof(std::complex<double> **) * nWSC);
257  indx_tot = (int **)malloc(sizeof(int*) * nWSC);
258  for (iWSC = 0; iWSC < nWSC; iWSC++) {
259  Mat_tot[iWSC] = (std::complex<double> **)malloc(sizeof(std::complex<double> *) * nWan);
260  indx_tot[iWSC] = (int *)malloc(sizeof(int) * 3);
261  for (iWan = 0; iWan < nWan; iWan++) {
262  Mat_tot[iWSC][iWan] = (std::complex<double> *)malloc(sizeof(std::complex<double>) * nWan);
263  }
264  }
265  /*
266  Read body
267  */
268  for (ii = 0; ii < 3; ii++) {
269  Rmin[ii] = 0;
270  Rmax[ii] = 0;
271  }
272  for (iWSC = 0; iWSC < nWSC; iWSC++) {
273  for (iWan = 0; iWan < nWan; iWan++) {
274  for (jWan = 0; jWan < nWan; jWan++) {
275  ierr = fscanf(fp, "%d%d%d%d%d%lf%lf",
276  &indx_tot[iWSC][0], &indx_tot[iWSC][1], &indx_tot[iWSC][2],
277  &iWan0, &jWan0,
278  &dtmp[0], &dtmp[1]);
279  if (iWan0 <= StdI->NsiteUC && jWan0 <= StdI->NsiteUC)
280  Mat_tot[iWSC][iWan0 - 1][jWan0 - 1] = std::complex<double>(dtmp[0], dtmp[1]);
281  for (ii = 0; ii < 3; ii++) {
282  if (indx_tot[iWSC][ii] < Rmin[ii]) Rmin[ii] = indx_tot[iWSC][ii];
283  if (indx_tot[iWSC][ii] > Rmax[ii]) Rmax[ii] = indx_tot[iWSC][ii];
284  }
285  }
286  }
287  }/*for (iWSC = 0; iWSC < nWSC; iWSC++)*/
288  fclose(fp);
289  for (ii = 0; ii < 3; ii++) NR[ii] = Rmax[ii] - Rmin[ii] + 1;
290  fprintf(stdout, " Minimum R : %d %d %d\n", Rmin[0], Rmin[1], Rmin[2]);
291  fprintf(stdout, " Maximum R : %d %d %d\n", Rmax[0], Rmax[1], Rmax[2]);
292  fprintf(stdout, " Numver of R : %d %d %d\n", NR[0], NR[1], NR[2]);
293 
294  DenMat = (std::complex<double> *****)malloc(sizeof(std::complex<double> ****)*NR[0]);
295  DenMat = DenMat - Rmin[0]; // base shift
296  for (i0 = Rmin[0]; i0 <= Rmax[0]; i0++) {
297  DenMat[i0] = (std::complex<double> ****)malloc(sizeof(std::complex<double> ***)*NR[1]);
298  DenMat[i0] = DenMat[i0] - Rmin[1]; // base shift
299  for (i1 = Rmin[1]; i1 <= Rmax[1]; i1++) {
300  DenMat[i0][i1] = (std::complex<double> ***)malloc(sizeof(std::complex<double> **)*NR[2]);
301  DenMat[i0][i1] = DenMat[i0][i1] - Rmin[2]; // base shift
302  for (i2 = Rmin[2]; i2 <= Rmax[2]; i2++) {
303  DenMat[i0][i1][i2] = (std::complex<double> **)malloc(sizeof(std::complex<double> *) * StdI->NsiteUC);
304  for (iWan = 0; iWan < StdI->NsiteUC; iWan++) {
305  DenMat[i0][i1][i2][iWan] = (std::complex<double> *)malloc(sizeof(std::complex<double>) * StdI->NsiteUC);
306  for (jWan = 0; jWan < StdI->NsiteUC; jWan++)DenMat[i0][i1][i2][iWan][jWan] = 0.0;
307  }
308  }
309  }
310  }
311  for (iWSC = 0; iWSC < nWSC; iWSC++)
312  for (iWan = 0; iWan < nWan; iWan++)
313  for (jWan = 0; jWan < nWan; jWan++)
314  DenMat[indx_tot[iWSC][0]][indx_tot[iWSC][1]][indx_tot[iWSC][2]][iWan][jWan]
315  = Mat_tot[iWSC][iWan][jWan];
316 
317  for (iWSC = 0; iWSC < nWSC; iWSC++) {
318  for (iWan = 0; iWan < nWan; iWan++) {
319  free(Mat_tot[iWSC][iWan]);
320  }
321  free(Mat_tot[iWSC]);
322  free(indx_tot[iWSC]);
323  }
324  free(Mat_tot);
325  free(indx_tot);
326 
327  return DenMat;
328 }/*static int read_W90(struct StdIntList *StdI, char *model)*/

◆ read_W90()

static void read_W90 ( struct StdIntList *  StdI,
char *  filename,
double  cutoff,
int *  cutoff_R,
double  cutoff_length,
int  itUJ,
int *  NtUJ,
int ***  tUJindx,
std::complex< double > **  tUJ 
)
static

Read Wannier90 hamiltonian file (*_hr)

Author
Mitsuaki Kawamura (The University of Tokyo)

(1) Apply inversion symmetry and delete duplication

(3-1) Compute the number of terms lerger than cut-off.

Then Store to the hopping Integeral and its site index.

Parameters
[in,out]StdI
[in]filenameInput file name
[in]cutoffThreshold for the Hamiltonian
[out]tUJindxR, band index of matrix element
[out]tUJMatrix element

Definition at line 77 of file Wannier90.cpp.

Referenced by StdFace_Wannier90().

88 {
89  FILE *fp;
90  int ierr, nWan, nWSC, iWSC, jWSC, iWan, jWan, iWan0, jWan0, ii, jj;
91  double dtmp[2], dR[3], length;
92  char ctmp[256], *ctmp2;
93  std::complex<double> ***Mat_tot;
94  int **indx_tot;
95  /*
96  Header part
97  */
98  fp = fopen(filename, "r");
99  ctmp2 = fgets(ctmp, 256, fp);
100  ierr = fscanf(fp, "%d", &nWan);
101  if (ierr == EOF) printf("%d %s\n", ierr, ctmp2);
102  ierr = fscanf(fp, "%d", &nWSC);
103  for (iWSC = 0; iWSC < nWSC; iWSC++) {
104  ierr = fscanf(fp, "%d", &ii);
105  }
106  /*
107  Malloc Matrix elements and their indices
108  */
109  Mat_tot = (std::complex<double> ***)malloc(sizeof(std::complex<double> **) * nWSC);
110  indx_tot = (int **)malloc(sizeof(int*) * nWSC);
111  for (iWSC = 0; iWSC < nWSC; iWSC++) {
112  Mat_tot[iWSC] = (std::complex<double> **)malloc(sizeof(std::complex<double> *) * nWan);
113  indx_tot[iWSC] = (int *)malloc(sizeof(int) * 3);
114  for (iWan = 0; iWan < nWan; iWan++) {
115  Mat_tot[iWSC][iWan] = (std::complex<double> *)malloc(sizeof(std::complex<double>) * nWan);
116  }
117  }
118  /*
119  Read body
120  */
121  for (iWSC = 0; iWSC < nWSC; iWSC++) {
122  for (iWan = 0; iWan < nWan; iWan++) {
123  for (jWan = 0; jWan < nWan; jWan++) {
124  ierr = fscanf(fp, "%d%d%d%d%d%lf%lf",
125  &indx_tot[iWSC][0], &indx_tot[iWSC][1], &indx_tot[iWSC][2],
126  &iWan0, &jWan0,
127  &dtmp[0], &dtmp[1]);
128  /*
129  Compute Euclid length
130  */
131  for (ii = 0; ii < 3; ii++) {
132  dR[ii] = 0.0;
133  for (jj = 0; jj < 3; jj++)
134  dR[ii] += StdI->direct[jj][ii] * (StdI->tau[jWan][jj] - StdI->tau[iWan][jj] + indx_tot[iWSC][jj]);
135  }
136  length = sqrt(dR[0] * dR[0] + dR[1] * dR[1] + dR[2] * dR[2]);
137  if (length > cutoff_length && cutoff_length > 0.0) {
138  dtmp[0] = 0.0;
139  dtmp[1] = 0.0;
140  }
141  if (abs(indx_tot[iWSC][0]) > cutoff_R[0] ||
142  abs(indx_tot[iWSC][1]) > cutoff_R[1] ||
143  abs(indx_tot[iWSC][2]) > cutoff_R[2]) {
144  dtmp[0] = 0.0;
145  dtmp[1] = 0.0;
146  }
147  if (iWan0 <= StdI->NsiteUC && jWan0 <= StdI->NsiteUC)
148  Mat_tot[iWSC][iWan0 - 1][jWan0 - 1] = std::complex<double>(dtmp[0], dtmp[1]);
149  }
150  }
154  for (jWSC = 0; jWSC < iWSC; jWSC++) {
155  if (
156  indx_tot[iWSC][0] == -indx_tot[jWSC][0] &&
157  indx_tot[iWSC][1] == -indx_tot[jWSC][1] &&
158  indx_tot[iWSC][2] == -indx_tot[jWSC][2]
159  )
160  for (iWan = 0; iWan < StdI->NsiteUC; iWan++) {
161  for (jWan = 0; jWan < StdI->NsiteUC; jWan++) {
162  Mat_tot[iWSC][iWan][jWan] = 0.0;
163  }
164  }
165  }/*for (jWSC = 0; jWSC < iWSC; jWSC++)*/
166  if (indx_tot[iWSC][0] == 0 &&
167  indx_tot[iWSC][1] == 0 &&
168  indx_tot[iWSC][2] == 0)
169  for (iWan = 0; iWan < StdI->NsiteUC; iWan++) {
170  for (jWan = 0; jWan < iWan; jWan++) {
171  Mat_tot[iWSC][iWan][jWan] = 0.0;
172  }
173  }
174  }/*for (iWSC = 0; iWSC < nWSC; iWSC++)*/
175  fclose(fp);
179  fprintf(stdout, "\n EFFECTIVE terms:\n");
180  fprintf(stdout, " R0 R1 R2 band_i band_f Hamiltonian\n");
181  NtUJ[itUJ] = 0;
182  for (iWSC = 0; iWSC < nWSC; iWSC++) {
183  for (iWan = 0; iWan < StdI->NsiteUC; iWan++) {
184  for (jWan = 0; jWan < StdI->NsiteUC; jWan++) {
185  if (cutoff < abs(Mat_tot[iWSC][iWan][jWan])) {
186  fprintf(stdout, " %5d%5d%5d%5d%5d%12.6f%12.6f\n",
187  indx_tot[iWSC][0], indx_tot[iWSC][1], indx_tot[iWSC][2], iWan, jWan,
188  real(Mat_tot[iWSC][iWan][jWan]), imag(Mat_tot[iWSC][iWan][jWan]));
189  NtUJ[itUJ] += 1;
190  }
191  }
192  }
193  }
194  fprintf(stdout, " Total number of EFFECTIVE term = %d\n", NtUJ[itUJ]);
195  tUJ[itUJ] = (std::complex<double> *)malloc(sizeof(std::complex<double>) * NtUJ[itUJ]);
196  tUJindx[itUJ] = (int **)malloc(sizeof(int*) * NtUJ[itUJ]);
197  for (ii = 0; ii < NtUJ[itUJ]; ii++) tUJindx[itUJ][ii] = (int *)malloc(sizeof(int) * 5);
202  NtUJ[itUJ] = 0;
203  for (iWSC = 0; iWSC < nWSC; iWSC++) {
204  for (iWan = 0; iWan < StdI->NsiteUC; iWan++) {
205  for (jWan = 0; jWan < StdI->NsiteUC; jWan++) {
206  if (cutoff < abs(Mat_tot[iWSC][iWan][jWan])) {
207  for (ii = 0; ii < 3; ii++) tUJindx[itUJ][NtUJ[itUJ]][ii] = indx_tot[iWSC][ii];
208  tUJindx[itUJ][NtUJ[itUJ]][3] = iWan;
209  tUJindx[itUJ][NtUJ[itUJ]][4] = jWan;
210  tUJ[itUJ][NtUJ[itUJ]] = Mat_tot[iWSC][iWan][jWan];
211  NtUJ[itUJ] += 1;
212  }
213  }/*for (jWan = 0; jWan < StdI->NsiteUC; jWan++)*/
214  }/*for (iWan = 0; iWan < StdI->NsiteUC; iWan++)*/
215  }/*for (iWSC = 0; iWSC < nWSC; iWSC++)*/
216 
217  for (iWSC = 0; iWSC < nWSC; iWSC++) {
218  for (iWan = 0; iWan < nWan; iWan++) {
219  free(Mat_tot[iWSC][iWan]);
220  }
221  free(Mat_tot[iWSC]);
222  free(indx_tot[iWSC]);
223  }
224  free(Mat_tot);
225  free(indx_tot);
226 }/*static int read_W90(struct StdIntList *StdI, char *model)*/

◆ StdFace_Wannier90()

void StdFace_Wannier90 ( struct StdIntList *  StdI)

Setup a Hamiltonian for the Wannier90 *_hr.dat.

Author
Mitsuaki Kawamura (The University of Tokyo)

(1) Compute the shape of the super-cell and sites in the super-cell

(2) check & store parameters of Hamiltonian

(3) Set local spin flag (StdIntList::locspinflag) and the number of sites (StdIntList::nsite)

(4) Compute the upper limit of the number of Transfer & Interaction and malloc them.

(4.5) For spin system, compute super exchange interaction.

(5) Set Transfer & Interaction

Parameters
[in,out]StdI

Definition at line 424 of file Wannier90.cpp.

References geometry_W90(), PrintUHFinitial(), read_density_matrix(), read_W90(), StdFace_Coulomb(), StdFace_exit(), StdFace_FindSite(), StdFace_GeneralJ(), StdFace_Hopping(), StdFace_HubbardLocal(), StdFace_InitSite(), StdFace_MagField(), StdFace_MallocInteractions(), StdFace_NotUsed_d(), StdFace_PrintGeometry(), StdFace_PrintVal_d(), StdFace_PrintVal_i(), and StdFace_PrintXSF().

Referenced by StdFace_main().

427 {
428  int isite, jsite, ispin, ntransMax, nintrMax;
429  int iL, iW, iH, kCell, it, ii;
430  double Jtmp[3][3] = { {0.0} };
431  FILE *fp;
432  std::complex<double> Cphase, DenMat0;
433  double dR[3], *Uspin;
434  int NtUJ[3];
435  std::complex<double> **tUJ, *****DenMat;
436  int ***tUJindx;
437  char filename[263];
441  fp = fopen("lattice.xsf", "w");
442 
443  StdFace_PrintVal_d("phase0", &StdI->phase[0], 0.0);
444  StdFace_PrintVal_d("phase1", &StdI->phase[1], 0.0);
445  StdFace_PrintVal_d("phase2", &StdI->phase[2], 0.0);
446  StdI->NsiteUC = 1;
447  StdFace_InitSite(StdI, fp, 3);
448  fprintf(stdout, "\n @ Wannier90 Geometry \n\n");
449  geometry_W90(StdI);
450  StdFace_PrintVal_i("DoubleCounting", &StdI->double_counting, 1);
451 
452  tUJ = (std::complex<double> **)malloc(sizeof(std::complex<double>*) * 3);
453  tUJindx = (int ***)malloc(sizeof(int**) * 3);
454 
455  /*
456  Read Hopping
457  */
458  fprintf(stdout, "\n @ Wannier90 hopping \n\n");
459  StdFace_PrintVal_d("cutoff_t", &StdI->cutoff_t, 1.0e-8);
460  StdFace_PrintVal_d("cutoff_length_t", &StdI->cutoff_length_t, -1.0);
461  sprintf(filename, "%s_hr.dat", StdI->CDataFileHead);
462  read_W90(StdI, filename,
463  StdI->cutoff_t, StdI->cutoff_tR, StdI->cutoff_length_t,
464  0, NtUJ, tUJindx, tUJ);
465  /*
466  Read Coulomb
467  */
468  fprintf(stdout, "\n @ Wannier90 Coulomb \n\n");
469  StdFace_PrintVal_d("cutoff_u", &StdI->cutoff_u, 1.0e-8);
470  StdFace_PrintVal_d("cutoff_length_U", &StdI->cutoff_length_U, -1.0);
471  sprintf(filename, "%s_ur.dat", StdI->CDataFileHead);
472  read_W90(StdI, filename,
473  StdI->cutoff_u, StdI->cutoff_UR, StdI->cutoff_length_U,
474  1, NtUJ, tUJindx, tUJ);
475  /*
476  Read Hund
477  */
478  fprintf(stdout, "\n @ Wannier90 Hund \n\n");
479  StdFace_PrintVal_d("cutoff_j", &StdI->cutoff_j, 1.0e-8);
480  StdFace_PrintVal_d("cutoff_length_J", &StdI->cutoff_length_J, -1.0);
481  sprintf(filename, "%s_jr.dat", StdI->CDataFileHead);
482  read_W90(StdI, filename,
483  StdI->cutoff_j, StdI->cutoff_JR, StdI->cutoff_length_J,
484  2, NtUJ, tUJindx, tUJ);
485  /*
486  Read Density matrix
487  */
488  fprintf(stdout, "\n @ Wannier90 Density-matrix \n\n");
489  sprintf(filename, "%s_dr.dat", StdI->CDataFileHead);
490  DenMat = read_density_matrix(StdI, filename);
494  fprintf(stdout, "\n @ Hamiltonian \n\n");
495  StdFace_NotUsed_d("K", StdI->K);
496  StdFace_PrintVal_d("h", &StdI->h, 0.0);
497  StdFace_PrintVal_d("Gamma", &StdI->Gamma, 0.0);
498  StdFace_NotUsed_d("U", StdI->U);
499 
500  if (strcmp(StdI->model, "spin") == 0 ) {
501  StdFace_PrintVal_i("2S", &StdI->S2, 1);
502  }/*if (strcmp(StdI->model, "spin") == 0 )*/
503  else if (strcmp(StdI->model, "hubbard") == 0) {
504  StdFace_PrintVal_d("mu", &StdI->mu, 0.0);
505  }
506  else{
507  printf("wannier + Kondo is not available !\n");
508  StdFace_exit(-1);
509  }/*if (model != "spin")*/
510  fprintf(stdout, "\n @ Numerical conditions\n\n");
515  StdI->nsite = StdI->NsiteUC * StdI->NCell;
516  StdI->locspinflag = (int *)malloc(sizeof(int) * StdI->nsite);
517 
518  if(strcmp(StdI->model, "spin") == 0 )
519  for (isite = 0; isite < StdI->nsite; isite++) StdI->locspinflag[isite] = StdI->S2;
520  else if(strcmp(StdI->model, "hubbard") == 0 )
521  for (isite = 0; isite < StdI->nsite; isite++) StdI->locspinflag[isite] = 0;
525  if (strcmp(StdI->model, "spin") == 0 ) {
526  ntransMax = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/);
527  nintrMax = StdI->NCell * (StdI->NsiteUC/*D*/ + NtUJ[0]/*J*/ + NtUJ[1] + NtUJ[2])
528  * (3 * StdI->S2 + 1) * (3 * StdI->S2 + StdI->NsiteUC);
529  }
530  else if (strcmp(StdI->model, "hubbard") == 0) {
531  ntransMax = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + NtUJ[0] * 2/*t*/
532  + NtUJ[1] * 2 * 3/*DC(U)*/ + NtUJ[2] * 2 * 2/*DC(J)*/);
533  nintrMax = StdI->NCell * (NtUJ[1] + NtUJ[2] + StdI->NsiteUC);
534  }
535 
536  StdFace_MallocInteractions(StdI, ntransMax, nintrMax);
540  if (strcmp(StdI->model, "spin") == 0) {
541  Uspin = (double *)malloc(sizeof(double) * StdI->NsiteUC);
542  for (it = 0; it < NtUJ[1]; it++)
543  if (tUJindx[1][it][0] == 0 && tUJindx[1][it][1] == 0 && tUJindx[1][it][2] == 0
544  && tUJindx[1][it][3] == tUJindx[1][it][4])
545  Uspin[tUJindx[1][it][3]] = real(tUJ[1][it]);
546  }/*if (strcmp(StdI->model, "spin") == 0)*/
550  for (kCell = 0; kCell < StdI->NCell; kCell++){
551 
552  iW = StdI->Cell[kCell][0];
553  iL = StdI->Cell[kCell][1];
554  iH = StdI->Cell[kCell][2];
555  /*
556  Local term 1
557  */
558  if (strcmp(StdI->model, "spin") == 0) {
559  for (isite = StdI->NsiteUC*kCell; isite < StdI->NsiteUC*(kCell + 1); isite++) {
560  StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, isite);
561  }
562  }/*if (strcmp(StdI->model, "spin") == 0 )*/
563  else {
564  for (isite = StdI->NsiteUC*kCell; isite < StdI->NsiteUC*(kCell + 1); isite++) {
565  StdFace_HubbardLocal(StdI, StdI->mu, -StdI->h, -StdI->Gamma, 0.0, isite);
566  }
567  }/*if (strcmp(StdI->model, "spin") != 0 )*/
568  /*
569  Hopping
570  */
571  for (it = 0; it < NtUJ[0]; it++) {
572  /*
573  Local term
574  */
575  if (tUJindx[0][it][0] == 0 && tUJindx[0][it][1] == 0 && tUJindx[0][it][2] == 0
576  && tUJindx[0][it][3] == tUJindx[0][it][4])
577  {
578  if (strcmp(StdI->model, "hubbard") == 0) {
579  isite = StdI->NsiteUC*kCell + tUJindx[0][it][3];
580  for (ispin = 0; ispin < 2; ispin++) {
581  StdI->trans[StdI->ntrans] = -tUJ[0][it];
582  StdI->transindx[StdI->ntrans][0] = isite;
583  StdI->transindx[StdI->ntrans][1] = ispin;
584  StdI->transindx[StdI->ntrans][2] = isite;
585  StdI->transindx[StdI->ntrans][3] = ispin;
586  StdI->ntrans = StdI->ntrans + 1;
587  }/*for (ispin = 0; ispin < 2; ispin++)*/
588  }/*if (strcmp(StdI->model, "hubbrad") == 0 )*/
589  }/*Local term*/
590  else {
591  /*
592  Non-local term
593  */
594  StdFace_FindSite(StdI, iW, iL, iH,
595  tUJindx[0][it][0], tUJindx[0][it][1], tUJindx[0][it][2],
596  tUJindx[0][it][3], tUJindx[0][it][4], &isite, &jsite, &Cphase, dR);
597  if (strcmp(StdI->model, "spin") == 0) {
598  for (ii = 0; ii < 3; ii++)
599  Jtmp[ii][ii] = 2.0 * real(tUJ[0][it] * conj(tUJ[0][it]))
600  * (1.0 / Uspin[tUJindx[0][it][3]] + 1.0 / Uspin[tUJindx[0][it][4]]);
601  StdFace_GeneralJ(StdI, Jtmp, StdI->S2, StdI->S2, isite, jsite);
602  }/*if (strcmp(StdI->model, "spin") == 0 )*/
603  else {
604  StdFace_Hopping(StdI, - Cphase * tUJ[0][it], jsite, isite, dR);
605  }
606  }/*Non-local term*/
607  }/*for (it = 0; it < NtUJ[0]; it++)*/
608  /*
609  Coulomb integral (U)
610  */
611  for (it = 0; it < NtUJ[1]; it++) {
612  /*
613  Local term
614  */
615  if (tUJindx[1][it][0] == 0 && tUJindx[1][it][1] == 0 && tUJindx[1][it][2] == 0
616  && tUJindx[1][it][3] == tUJindx[1][it][4])
617  {
618  StdI->Cintra[StdI->NCintra] = real(tUJ[1][it]);
619  StdI->CintraIndx[StdI->NCintra][0] = StdI->NsiteUC*kCell + tUJindx[1][it][3];
620  StdI->NCintra += 1;
621  /*
622  Double-counting correction @f$0.5 U_{0ii} D_{0ii}@f$
623  */
624  if (StdI->double_counting == 1) {
625  isite = StdI->NsiteUC*kCell + tUJindx[1][it][3];
626  for (ispin = 0; ispin < 2; ispin++) {
627  DenMat0 = DenMat[0][0][0][tUJindx[1][it][3]][tUJindx[1][it][3]];
628  StdI->trans[StdI->ntrans] = 0.5*real(tUJ[1][it])*DenMat0;
629  StdI->transindx[StdI->ntrans][0] = isite;
630  StdI->transindx[StdI->ntrans][1] = ispin;
631  StdI->transindx[StdI->ntrans][2] = isite;
632  StdI->transindx[StdI->ntrans][3] = ispin;
633  StdI->ntrans = StdI->ntrans + 1;
634  }/*for (ispin = 0; ispin < 2; ispin++)*/
635  }
636  }/*Local term*/
637  else {
638  /*
639  Non-local term
640  */
641  StdFace_FindSite(StdI, iW, iL, iH,
642  tUJindx[1][it][0], tUJindx[1][it][1], tUJindx[1][it][2],
643  tUJindx[1][it][3], tUJindx[1][it][4], &isite, &jsite, &Cphase, dR);
644  StdFace_Coulomb(StdI, real(tUJ[1][it]), isite, jsite);
645  /*
646  Double-counting correction
647  */
648  if (StdI->double_counting == 1) {
649  for (ispin = 0; ispin < 2; ispin++) {
650  /*
651  @f$sum_{(R,j)(>0,i)} U_{Rij} D_{0jj} (Local)@f$
652  */
653  DenMat0 = DenMat[0][0][0][tUJindx[1][it][4]][tUJindx[1][it][4]];
654  StdI->trans[StdI->ntrans] = real(tUJ[1][it])*DenMat0;
655  StdI->transindx[StdI->ntrans][0] = isite;
656  StdI->transindx[StdI->ntrans][1] = ispin;
657  StdI->transindx[StdI->ntrans][2] = isite;
658  StdI->transindx[StdI->ntrans][3] = ispin;
659  StdI->ntrans = StdI->ntrans + 1;
660  /*
661  @f$sum_{(R,j)(>0,i)} U_{Rij} D_{0jj} (Local)@f$
662  */
663  DenMat0 = DenMat[0][0][0][tUJindx[1][it][3]][tUJindx[1][it][3]];
664  StdI->trans[StdI->ntrans] = real(tUJ[1][it])*DenMat0;
665  StdI->transindx[StdI->ntrans][0] = jsite;
666  StdI->transindx[StdI->ntrans][1] = ispin;
667  StdI->transindx[StdI->ntrans][2] = jsite;
668  StdI->transindx[StdI->ntrans][3] = ispin;
669  StdI->ntrans = StdI->ntrans + 1;
670  }/*for (ispin = 0; ispin < 2; ispin++)*/
671  /*
672  @f$-0.5U_{Rij} D_{Rjj}@f$
673  */
674  DenMat0 = DenMat[tUJindx[1][it][0]][tUJindx[1][it][1]][tUJindx[1][it][2]]
675  [tUJindx[1][it][3]][tUJindx[1][it][4]];
676  StdFace_Hopping(StdI, -0.5*Cphase * real(tUJ[1][it])*DenMat0, jsite, isite, dR);
677  }/*if (StdI->double_counting == 1)*/
678  }/*Non-local term*/
679  }/*for (it = 0; it < NtUJ[0]; it++)*/
680  /*
681  Hund coupling (J)
682  */
683  for (it = 0; it < NtUJ[2]; it++) {
684  /*
685  Local term should not be computed
686  */
687  if (tUJindx[2][it][0] != 0 || tUJindx[2][it][1] != 0 || tUJindx[2][it][2] != 0
688  || tUJindx[2][it][3] != tUJindx[2][it][4])
689  {
690  StdFace_FindSite(StdI, iW, iL, iH,
691  tUJindx[2][it][0], tUJindx[2][it][1], tUJindx[2][it][2],
692  tUJindx[2][it][3], tUJindx[2][it][4], &isite, &jsite, &Cphase, dR);
693 
694  StdI->Hund[StdI->NHund] = real(tUJ[2][it]);
695  StdI->HundIndx[StdI->NHund][0] = isite;
696  StdI->HundIndx[StdI->NHund][1] = jsite;
697  StdI->NHund += 1;
698 
699  if (strcmp(StdI->model, "hubbard") == 0) {
700  StdI->Ex[StdI->NEx] = real(tUJ[2][it]);
701  StdI->ExIndx[StdI->NEx][0] = isite;
702  StdI->ExIndx[StdI->NEx][1] = jsite;
703  StdI->NEx += 1;
704 
705  StdI->PairHopp[StdI->NPairHopp] = real(tUJ[2][it]);
706  StdI->PHIndx[StdI->NPairHopp][0] = isite;
707  StdI->PHIndx[StdI->NPairHopp][1] = jsite;
708  StdI->NPairHopp += 1;
709  /*
710  Double-counting correction
711  */
712  if (StdI->double_counting == 1) {
713  for (ispin = 0; ispin < 2; ispin++) {
714  /*
715  @f$- \frac{1}{2}sum_{(R,j)(>0,i)} J_{Rij} D_{0jj}@f$
716  */
717  DenMat0 = DenMat[0][0][0][tUJindx[2][it][4]][tUJindx[2][it][4]];
718  StdI->trans[StdI->ntrans] = -0.5*real(tUJ[2][it]) *DenMat0;
719  StdI->transindx[StdI->ntrans][0] = isite;
720  StdI->transindx[StdI->ntrans][1] = ispin;
721  StdI->transindx[StdI->ntrans][2] = isite;
722  StdI->transindx[StdI->ntrans][3] = ispin;
723  StdI->ntrans = StdI->ntrans + 1;
724  /*
725  @f$- \frac{1}{2}sum_{(R,j)(>0,i)} J_{Rij} D_{0jj}@f$
726  */
727  DenMat0 = DenMat[0][0][0][tUJindx[2][it][3]][tUJindx[2][it][3]];
728  StdI->trans[StdI->ntrans] = -0.5*real(tUJ[2][it]) *DenMat0;
729  StdI->transindx[StdI->ntrans][0] = jsite;
730  StdI->transindx[StdI->ntrans][1] = ispin;
731  StdI->transindx[StdI->ntrans][2] = jsite;
732  StdI->transindx[StdI->ntrans][3] = ispin;
733  StdI->ntrans = StdI->ntrans + 1;
734  }/*for (ispin = 0; ispin < 2; ispin++)*/
735  /*
736  @f$J_{Rij} (D_{Rjj}+2{\rm Re}[D_{Rjj])@f$
737  */
738  DenMat0 = DenMat[tUJindx[2][it][0]][tUJindx[2][it][1]][tUJindx[2][it][2]]
739  [tUJindx[2][it][3]][tUJindx[2][it][4]];
740  StdFace_Hopping(StdI,
741  0.5*Cphase * real(tUJ[2][it])*(DenMat0 + 2.0*real(DenMat0)), jsite, isite, dR);
742  }/*if (StdI->double_counting == 1)*/
743  }
744  else {
745 #if defined(_mVMC)
746  StdI->Ex[StdI->NEx] = real(tUJ[2][it]);
747 #else
748  StdI->Ex[StdI->NEx] = -real(tUJ[2][it]);
749 #endif
750  StdI->ExIndx[StdI->NEx][0] = isite;
751  StdI->ExIndx[StdI->NEx][1] = jsite;
752  StdI->NEx += 1;
753  }
754  }/*Non-local term*/
755  }/*for (it = 0; it < NtUJ[0]; it++)*/
756  }/*for (kCell = 0; kCell < StdI->NCell; kCell++)*/
757 
758  fclose(fp);
759  PrintUHFinitial(StdI, NtUJ, DenMat, tUJindx);
760  StdFace_PrintXSF(StdI);
761  StdFace_PrintGeometry(StdI);
762 
763  for (it = 0; it < NtUJ[0]; it++) free(tUJindx[0][it]);
764  free(tUJindx[0]);
765  free(tUJ[0]);
766  for (it = 0; it < NtUJ[1]; it++) free(tUJindx[1][it]);
767  free(tUJindx[1]);
768  free(tUJ[1]);
769  for (it = 0; it < NtUJ[2]; it++) free(tUJindx[2][it]);
770  free(tUJindx[2]);
771  free(tUJ[2]);
772  if (strcmp(StdI->model, "spin") == 0) free(Uspin);
773 
774 }/*void StdFace_Wannier90*/
read_W90
static void read_W90(struct StdIntList *StdI, char *filename, double cutoff, int *cutoff_R, double cutoff_length, int itUJ, int *NtUJ, int ***tUJindx, std::complex< double > **tUJ)
Read Wannier90 hamiltonian file (*_hr)
Definition: Wannier90.cpp:77
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_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_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
geometry_W90
static void geometry_W90(struct StdIntList *StdI)
Read Geometry file for wannier90.
Definition: Wannier90.cpp:32
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_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_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_MallocInteractions
void StdFace_MallocInteractions(struct StdIntList *StdI, int ntransMax, int nintrMax)
Malloc Arrays for interactions.
Definition: StdFace_ModelUtil.cpp:1204
PrintUHFinitial
static void PrintUHFinitial(struct StdIntList *StdI, int *NtUJ, std::complex< double > *****DenMat, int ***tUJindx)
Print the initial guess of UHF.
Definition: Wannier90.cpp:333
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_InitSite
void StdFace_InitSite(struct StdIntList *StdI, FILE *fp, int dim)
Initialize the super-cell where simulation is performed.
Definition: StdFace_ModelUtil.cpp:636
read_density_matrix
static std::complex< double > ***** read_density_matrix(struct StdIntList *StdI, char *filename)
Read RESPACK Density-matrix file (*_dr.dat)
Definition: Wannier90.cpp:231
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