HΦ  3.2.0
expec_energy_flct.c
Go to the documentation of this file.
1 /* HPhi - Quantum Lattice Model Simulator */
2 /* Copyright (C) 2015 The University of Tokyo */
3 
4 /* This program is free software: you can redistribute it and/or modify */
5 /* it under the terms of the GNU General Public License as published by */
6 /* the Free Software Foundation, either version 3 of the License, or */
7 /* (at your option) any later version. */
8 
9 /* This program is distributed in the hope that it will be useful, */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
12 /* GNU General Public License for more details. */
13 
14 /* You should have received a copy of the GNU General Public License */
15 /* along with this program. If not, see <http://www.gnu.org/licenses/>. */
16 
17 #include "bitcalc.h"
18 #include "mltplyCommon.h"
19 #include "mltply.h"
20 #include "expec_energy_flct.h"
21 #include "wrapperMPI.h"
22 #include "CalcTime.h"
23 
34 int expec_energy_flct(struct BindStruct *X){
35 
36  long unsigned int i,j;
37  long unsigned int irght,ilft,ihfbit;
38  double complex dam_pr,dam_pr1;
39  long unsigned int i_max;
40 
41  switch(X->Def.iCalcType){
42  case Lanczos:
43  fprintf(stdoutMPI, "%s", cLogExpecEnergyStart);
44  TimeKeeper(X, cFileNameTimeKeep, cExpecStart, "a");
45  break;
46  case TPQCalc:
47  case TimeEvolution:
48 #ifdef _DEBUG
49  fprintf(stdoutMPI, "%s", cLogExpecEnergyStart);
50  TimeKeeperWithStep(X, cFileNameTimeKeep, cTPQExpecStart, "a", step_i);
51 #endif
52  break;
53  case FullDiag:
54  case CG:
55  break;
56  default:
57  return -1;
58  }
59 
60  i_max=X->Check.idim_max;
61  if(GetSplitBitByModel(X->Def.Nsite, X->Def.iCalcModel, &irght, &ilft, &ihfbit)!=0){
62  return -1;
63  }
64 
65  X->Large.i_max = i_max;
66  X->Large.irght = irght;
67  X->Large.ilft = ilft;
68  X->Large.ihfbit = ihfbit;
69  X->Large.mode = M_ENERGY;
70  X->Phys.energy=0.0;
71 
72  int nCalcFlct;
73  if(X->Def.iCalcType == Lanczos){
74  nCalcFlct=4301;
75  }
76  else if (X->Def.iCalcType == TPQCalc){
77  nCalcFlct=3201;
78  }
79  else{//For FullDiag
80  nCalcFlct=5301;
81  }
82  StartTimer(nCalcFlct);
83 
84  switch(X->Def.iCalcModel){
85  case HubbardGC:
86  expec_energy_flct_HubbardGC(X);
87  break;
88  case KondoGC:
89  case Hubbard:
90  case Kondo:
91  expec_energy_flct_Hubbard(X);
92  break;
93 
94  case SpinGC:
95  if(X->Def.iFlgGeneralSpin == FALSE) {
96  expec_energy_flct_HalfSpinGC(X);
97  }
98  else{//for generalspin
99  expec_energy_flct_GeneralSpinGC(X);
100  }
101  break;/*case SpinGC*/
102  /* SpinGCBoost */
103  case Spin:
104  /*
105  if(X->Def.iFlgGeneralSpin == FALSE){
106  expec_energy_flct_HalfSpin(X);
107  }
108  else{
109  expec_energy_flct_GeneralSpin(X);
110  }
111  */
112  X->Phys.doublon = 0.0;
113  X->Phys.doublon2 = 0.0;
114  X->Phys.num = X->Def.NsiteMPI;
115  X->Phys.num2 = X->Def.NsiteMPI*X->Def.NsiteMPI;
116  X->Phys.Sz = 0.5 * (double)X->Def.Total2SzMPI;
117  X->Phys.Sz2 = X->Phys.Sz * X->Phys.Sz;
118  break;
119  default:
120  return -1;
121  }
122 
123  StopTimer(nCalcFlct);
124 
125 #pragma omp parallel for default(none) private(i) shared(v1,v0) firstprivate(i_max)
126  for(i = 1; i <= i_max; i++){
127  v1[i]=v0[i];
128  v0[i]=0.0;
129  }
130 
131  int nCalcExpec;
132  if(X->Def.iCalcType == Lanczos){
133  nCalcExpec=4302;
134  }
135  else if (X->Def.iCalcType == TPQCalc){
136  nCalcExpec=3202;
137  }
138  else{//For FullDiag
139  nCalcExpec=5302;
140  }
141  StartTimer(nCalcExpec);
142  mltply(X, v0, v1); // v0+=H*v1
143  StopTimer(nCalcExpec);
144 /* switch -> SpinGCBoost */
145 
146  dam_pr=0.0;
147  dam_pr1=0.0;
148  #pragma omp parallel for default(none) reduction(+:dam_pr, dam_pr1) private(j) shared(v0, v1)firstprivate(i_max)
149  for(j=1;j<=i_max;j++){
150  dam_pr += conj(v1[j])*v0[j]; // E = <v1|H|v1>=<v1|v0>
151  dam_pr1 += conj(v0[j])*v0[j]; // E^2 = <v1|H*H|v1>=<v0|v0>
152  //v0[j]=v1[j]; v1-> orginal v0=H*v1
153  }
154  dam_pr = SumMPI_dc(dam_pr);
155  dam_pr1 = SumMPI_dc(dam_pr1);
156  //fprintf(stdoutMPI, "Debug: ene=%lf, var=%lf\n", creal(dam_pr), creal(dam_pr1));
157 
158  X->Phys.energy = dam_pr;
159  X->Phys.var = dam_pr1;
160 
161  switch(X->Def.iCalcType) {
162  case Lanczos:
163  fprintf(stdoutMPI, "%s", cLogExpecEnergyEnd);
164  TimeKeeper(X, cFileNameTimeKeep, cExpecEnd, "a");
165  break;
166  case TPQCalc:
167  case TimeEvolution:
168 #ifdef _DEBUG
169  fprintf(stdoutMPI, "%s", cLogExpecEnergyEnd);
170  TimeKeeperWithStep(X, cFileNameTimeKeep, cTPQExpecEnd, "a", step_i);
171 #endif
172  break;
173  default:
174  break;
175  }
176 
177  return 0;
178 }
179 
185 int expec_energy_flct_HubbardGC(struct BindStruct *X) {
186  long unsigned int j;
187  long unsigned int isite1;
188  long unsigned int is1_up_a, is1_up_b;
189  long unsigned int is1_down_a, is1_down_b;
190  int bit_up, bit_down, bit_D;
191  long unsigned int ibit_up, ibit_down, ibit_D;
192  double D, tmp_D, tmp_D2;
193  double N, tmp_N, tmp_N2;
194  double Sz, tmp_Sz, tmp_Sz2;
195  double tmp_v02;
196  long unsigned int i_max;
197  unsigned int l_ibit1, u_ibit1, i_32;
198  i_max=X->Check.idim_max;
199 
200  i_32 = 0xFFFFFFFF; //2^32 - 1
201  // tentative doublon
202  tmp_D = 0.0;
203  tmp_D2 = 0.0;
204  tmp_N = 0.0;
205  tmp_N2 = 0.0;
206  tmp_Sz = 0.0;
207  tmp_Sz2 = 0.0;
208 
209 //[s] for bit count
210  is1_up_a = 0;
211  is1_up_b = 0;
212  is1_down_a = 0;
213  is1_down_b = 0;
214  for (isite1 = 1; isite1 <= X->Def.NsiteMPI; isite1++) {
215  if (isite1 > X->Def.Nsite) {
216  is1_up_a += X->Def.Tpow[2 * isite1 - 2];
217  is1_down_a += X->Def.Tpow[2 * isite1 - 1];
218  } else {
219  is1_up_b += X->Def.Tpow[2 * isite1 - 2];
220  is1_down_b += X->Def.Tpow[2 * isite1 - 1];
221  }
222  }
223 //[e]
224 #pragma omp parallel for reduction(+:tmp_D,tmp_D2,tmp_N,tmp_N2,tmp_Sz,tmp_Sz2) default(none) shared(v0,list_1) \
225  firstprivate(i_max, X,myrank,is1_up_a,is1_down_a,is1_up_b,is1_down_b,i_32) \
226  private(j, tmp_v02,D,N,Sz,isite1,bit_up,bit_down,bit_D,u_ibit1,l_ibit1,ibit_up,ibit_down,ibit_D)
227  for (j = 1; j <= i_max; j++) {
228  tmp_v02 = conj(v0[j]) * v0[j];
229  bit_up = 0;
230  bit_down = 0;
231  bit_D = 0;
232 // isite1 > X->Def.Nsite
233  ibit_up = (unsigned long int) myrank & is1_up_a;
234  u_ibit1 = ibit_up >> 32;
235  l_ibit1 = ibit_up & i_32;
236  bit_up += pop(u_ibit1);
237  bit_up += pop(l_ibit1);
238 
239  ibit_down = (unsigned long int) myrank & is1_down_a;
240  u_ibit1 = ibit_down >> 32;
241  l_ibit1 = ibit_down & i_32;
242  bit_down += pop(u_ibit1);
243  bit_down += pop(l_ibit1);
244 
245  ibit_D = (ibit_up) & (ibit_down >> 1);
246  u_ibit1 = ibit_D >> 32;
247  l_ibit1 = ibit_D & i_32;
248  bit_D += pop(u_ibit1);
249  bit_D += pop(l_ibit1);
250 
251 // isite1 <= X->Def.Nsite
252  ibit_up = (unsigned long int) (j - 1) & is1_up_b;
253  u_ibit1 = ibit_up >> 32;
254  l_ibit1 = ibit_up & i_32;
255  bit_up += pop(u_ibit1);
256  bit_up += pop(l_ibit1);
257 
258  ibit_down = (unsigned long int) (j - 1) & is1_down_b;
259  u_ibit1 = ibit_down >> 32;
260  l_ibit1 = ibit_down & i_32;
261  bit_down += pop(u_ibit1);
262  bit_down += pop(l_ibit1);
263 
264  ibit_D = (ibit_up) & (ibit_down >> 1);
265  u_ibit1 = ibit_D >> 32;
266  l_ibit1 = ibit_D & i_32;
267  bit_D += pop(u_ibit1);
268  bit_D += pop(l_ibit1);
269 
270  D = bit_D;
271  N = bit_up + bit_down;
272  Sz = bit_up - bit_down;
273 
274  tmp_D += tmp_v02 * D;
275  tmp_D2 += tmp_v02 * D * D;
276  tmp_N += tmp_v02 * N;
277  tmp_N2 += tmp_v02 * N * N;
278  tmp_Sz += tmp_v02 * Sz;
279  tmp_Sz2 += tmp_v02 * Sz * Sz;
280  }
281  tmp_D = SumMPI_d(tmp_D);
282  tmp_D2 = SumMPI_d(tmp_D2);
283  tmp_N = SumMPI_d(tmp_N);
284  tmp_N2 = SumMPI_d(tmp_N2);
285  tmp_Sz = SumMPI_d(tmp_Sz);
286  tmp_Sz2 = SumMPI_d(tmp_Sz2);
287 
288  X->Phys.doublon = tmp_D;
289  X->Phys.doublon2 = tmp_D2;
290  X->Phys.num = tmp_N;
291  X->Phys.num2 = tmp_N2;
292  X->Phys.Sz = tmp_Sz*0.5;
293  X->Phys.Sz2 = tmp_Sz2*0.25;
294  X->Phys.num_up = 0.5*(tmp_N+tmp_Sz);
295  X->Phys.num_down = 0.5*(tmp_N-tmp_Sz);
296 
297  return 0;
298 }
299 
305 int expec_energy_flct_Hubbard(struct BindStruct *X){
306  long unsigned int j;
307  long unsigned int isite1;
308  long unsigned int is1_up_a,is1_up_b;
309  long unsigned int is1_down_a,is1_down_b;
310  int bit_up,bit_down,bit_D;
311 
312  long unsigned int ibit_up,ibit_down,ibit_D;
313  double D,tmp_D,tmp_D2;
314  double N,tmp_N,tmp_N2;
315  double Sz,tmp_Sz, tmp_Sz2;
316  double tmp_v02;
317  long unsigned int i_max,tmp_list_1;
318  unsigned int l_ibit1,u_ibit1,i_32;
319  i_max=X->Check.idim_max;
320 
321  i_32 = (unsigned int)(pow(2,32)-1);
322 
323  tmp_D = 0.0;
324  tmp_D2 = 0.0;
325  tmp_N = 0.0;
326  tmp_N2 = 0.0;
327  tmp_Sz = 0.0;
328  tmp_Sz2 = 0.0;
329 
330  //[s] for bit count
331  is1_up_a = 0;
332  is1_up_b = 0;
333  is1_down_a = 0;
334  is1_down_b = 0;
335  for(isite1=1;isite1<=X->Def.NsiteMPI;isite1++){
336  if(isite1 > X->Def.Nsite){
337  is1_up_a += X->Def.Tpow[2*isite1 - 2];
338  is1_down_a += X->Def.Tpow[2*isite1 - 1];
339  }else{
340  is1_up_b += X->Def.Tpow[2*isite1 - 2];
341  is1_down_b += X->Def.Tpow[2*isite1 - 1];
342  }
343  }
344 //[e]
345 #pragma omp parallel for reduction(+:tmp_D,tmp_D2,tmp_N,tmp_N2,tmp_Sz,tmp_Sz2) default(none) shared(v0,list_1) \
346  firstprivate(i_max, X,myrank,is1_up_a,is1_down_a,is1_up_b,is1_down_b,i_32) \
347  private(j, tmp_v02,D,N,Sz,isite1,tmp_list_1,bit_up,bit_down,bit_D,u_ibit1,l_ibit1,ibit_up,ibit_down,ibit_D)
348  for(j = 1; j <= i_max; j++) {
349  tmp_v02 = conj(v0[j]) * v0[j];
350  bit_up = 0;
351  bit_down = 0;
352  bit_D = 0;
353  tmp_list_1 = list_1[j];
354 // isite1 > X->Def.Nsite
355  ibit_up = (unsigned long int) myrank & is1_up_a;
356  u_ibit1 = ibit_up >> 32;
357  l_ibit1 = ibit_up & i_32;
358  bit_up += pop(u_ibit1);
359  bit_up += pop(l_ibit1);
360 
361  ibit_down = (unsigned long int) myrank & is1_down_a;
362  u_ibit1 = ibit_down >> 32;
363  l_ibit1 = ibit_down & i_32;
364  bit_down += pop(u_ibit1);
365  bit_down += pop(l_ibit1);
366 
367  ibit_D = (ibit_up) & (ibit_down >> 1);
368  u_ibit1 = ibit_D >> 32;
369  l_ibit1 = ibit_D & i_32;
370  bit_D += pop(u_ibit1);
371  bit_D += pop(l_ibit1);
372 
373 // isite1 <= X->Def.Nsite
374  ibit_up = (unsigned long int) tmp_list_1 & is1_up_b;
375  u_ibit1 = ibit_up >> 32;
376  l_ibit1 = ibit_up & i_32;
377  bit_up += pop(u_ibit1);
378  bit_up += pop(l_ibit1);
379 
380  ibit_down = (unsigned long int) tmp_list_1 & is1_down_b;
381  u_ibit1 = ibit_down >> 32;
382  l_ibit1 = ibit_down & i_32;
383  bit_down += pop(u_ibit1);
384  bit_down += pop(l_ibit1);
385 
386  ibit_D = (ibit_up) & (ibit_down >> 1);
387  u_ibit1 = ibit_D >> 32;
388  l_ibit1 = ibit_D & i_32;
389  bit_D += pop(u_ibit1);
390  bit_D += pop(l_ibit1);
391 
392  D = bit_D;
393  N = bit_up + bit_down;
394  Sz = bit_up - bit_down;
395 
396  tmp_D += tmp_v02 * D;
397  tmp_D2 += tmp_v02 * D * D;
398  tmp_N += tmp_v02 * N;
399  tmp_N2 += tmp_v02 * N * N;
400  tmp_Sz += tmp_v02 * Sz;
401  tmp_Sz2 += tmp_v02 * Sz * Sz;
402  }
403 
404 
405  tmp_D = SumMPI_d(tmp_D);
406  tmp_D2 = SumMPI_d(tmp_D2);
407  tmp_N = SumMPI_d(tmp_N);
408  tmp_N2 = SumMPI_d(tmp_N2);
409  tmp_Sz = SumMPI_d(tmp_Sz);
410  tmp_Sz2 = SumMPI_d(tmp_Sz2);
411 
412  X->Phys.doublon = tmp_D;
413  X->Phys.doublon2 = tmp_D2;
414  X->Phys.num = tmp_N;
415  X->Phys.num2 = tmp_N2;
416  X->Phys.Sz = tmp_Sz*0.5;
417  X->Phys.Sz2 = tmp_Sz2*0.25;
418  X->Phys.num_up = 0.5*(tmp_N+tmp_Sz);
419  X->Phys.num_down = 0.5*(tmp_N-tmp_Sz);
420  return 0;
421 }
422 
428 int expec_energy_flct_HalfSpinGC(struct BindStruct *X){
429  long unsigned int j;
430  long unsigned int isite1;
431  long unsigned int is1_up_a,is1_up_b;
432 
433  long unsigned int ibit1;
434  double Sz,tmp_Sz, tmp_Sz2;
435  double tmp_v02;
436  long unsigned int i_max;
437  unsigned int l_ibit1,u_ibit1,i_32;
438  i_max=X->Check.idim_max;
439 
440  i_32 = 0xFFFFFFFF; //2^32 - 1
441 
442  // tentative doublon
443  tmp_Sz = 0.0;
444  tmp_Sz2 = 0.0;
445 
446 //[s] for bit count
447  is1_up_a = 0;
448  is1_up_b = 0;
449  for(isite1=1;isite1<=X->Def.NsiteMPI;isite1++){
450  if(isite1 > X->Def.Nsite){
451  is1_up_a += X->Def.Tpow[isite1 - 1];
452  }else{
453  is1_up_b += X->Def.Tpow[isite1 - 1];
454  }
455  }
456 //[e]
457 #pragma omp parallel for reduction(+:tmp_Sz,tmp_Sz2)default(none) shared(v0) \
458  firstprivate(i_max,X,myrank,i_32,is1_up_a,is1_up_b) private(j,Sz,ibit1,isite1,tmp_v02,u_ibit1,l_ibit1)
459  for(j = 1; j <= i_max; j++){
460  tmp_v02 = conj(v0[j])*v0[j];
461  Sz = 0.0;
462 
463 // isite1 > X->Def.Nsite
464  ibit1 = (unsigned long int) myrank & is1_up_a;
465  u_ibit1 = ibit1 >> 32;
466  l_ibit1 = ibit1 & i_32;
467  Sz += pop(u_ibit1);
468  Sz += pop(l_ibit1);
469 // isite1 <= X->Def.Nsite
470  ibit1 = (unsigned long int) (j-1)&is1_up_b;
471  u_ibit1 = ibit1 >> 32;
472  l_ibit1 = ibit1 & i_32;
473  Sz += pop(u_ibit1);
474  Sz += pop(l_ibit1);
475  Sz = 2*Sz-X->Def.NsiteMPI;
476 
477  tmp_Sz += Sz*tmp_v02;
478  tmp_Sz2 += Sz*Sz*tmp_v02;
479  }
480  tmp_Sz = SumMPI_d(tmp_Sz);
481  tmp_Sz2 = SumMPI_d(tmp_Sz2);
482 
483  X->Phys.doublon = 0.0;
484  X->Phys.doublon2 = 0.0;
485  X->Phys.num = X->Def.NsiteMPI;
486  X->Phys.num2 = X->Def.NsiteMPI*X->Def.NsiteMPI;
487  X->Phys.Sz = tmp_Sz*0.5;
488  X->Phys.Sz2 = tmp_Sz2*0.25;
489  X->Phys.num_up = 0.5*(X->Def.NsiteMPI+tmp_Sz);
490  X->Phys.num_down = 0.5*(X->Def.NsiteMPI-tmp_Sz);
491 
492  return 0;
493 }
494 
500 int expec_energy_flct_GeneralSpinGC(struct BindStruct *X){
501  long unsigned int j;
502  long unsigned int isite1;
503 
504  double Sz,tmp_Sz, tmp_Sz2;
505  double tmp_v02;
506  long unsigned int i_max;
507  i_max=X->Check.idim_max;
508 
509  // tentative doublon
510  tmp_Sz = 0.0;
511  tmp_Sz2 = 0.0;
512 
513  for(j = 1; j <= i_max; j++){
514  tmp_v02 = conj(v0[j])*v0[j];
515  Sz = 0.0;
516  for(isite1=1;isite1<=X->Def.NsiteMPI;isite1++){
517  //prefactor 0.5 is added later.
518  if(isite1 > X->Def.Nsite){
519  Sz += GetLocal2Sz(isite1, myrank, X->Def.SiteToBit, X->Def.Tpow);
520  }else{
521  Sz += GetLocal2Sz(isite1, j-1, X->Def.SiteToBit, X->Def.Tpow);
522  }
523  }
524  tmp_Sz += Sz*tmp_v02;
525  tmp_Sz2 += Sz*Sz*tmp_v02;
526  }
527 
528  tmp_Sz = SumMPI_d(tmp_Sz);
529  tmp_Sz2 = SumMPI_d(tmp_Sz2);
530 
531  X->Phys.doublon = 0.0;
532  X->Phys.doublon2 = 0.0;
533  X->Phys.num = X->Def.NsiteMPI;
534  X->Phys.num2 = X->Def.NsiteMPI*X->Def.NsiteMPI;
535  X->Phys.Sz = tmp_Sz*0.5;
536  X->Phys.Sz2 = tmp_Sz2*0.25;
537  X->Phys.num_up = 0.5*(X->Def.NsiteMPI+tmp_Sz);
538  X->Phys.num_down = 0.5*(X->Def.NsiteMPI-tmp_Sz);
539 
540  return 0;
541 }
542 
548 int expec_energy_flct_HalfSpin(struct BindStruct *X){
549  long unsigned int j;
550  long unsigned int isite1;
551  long unsigned int is1_up_a,is1_up_b;
552 
553  long unsigned int ibit1;
554  double Sz,tmp_Sz, tmp_Sz2;
555  double tmp_v02;
556  long unsigned int i_max, tmp_list_1;
557  unsigned int l_ibit1,u_ibit1,i_32;
558  i_max=X->Check.idim_max;
559 
560  i_32 = 0xFFFFFFFF; //2^32 - 1
561 
562  // tentative doublon
563  tmp_Sz = 0.0;
564  tmp_Sz2 = 0.0;
565 
566 //[s] for bit count
567  is1_up_a = 0;
568  is1_up_b = 0;
569  for(isite1=1;isite1<=X->Def.NsiteMPI;isite1++){
570  if(isite1 > X->Def.Nsite){
571  is1_up_a += X->Def.Tpow[isite1 - 1];
572  }else{
573  is1_up_b += X->Def.Tpow[isite1 - 1];
574  }
575  }
576 //[e]
577 #pragma omp parallel for reduction(+:tmp_Sz,tmp_Sz2)default(none) shared(v0, list_1) \
578  firstprivate(i_max,X,myrank,i_32,is1_up_a,is1_up_b) private(j,Sz,ibit1,isite1,tmp_v02,u_ibit1,l_ibit1, tmp_list_1)
579  for(j = 1; j <= i_max; j++){
580  tmp_v02 = conj(v0[j])*v0[j];
581  Sz = 0.0;
582  tmp_list_1 = list_1[j];
583 
584 // isite1 > X->Def.Nsite
585  ibit1 = (unsigned long int) myrank & is1_up_a;
586  u_ibit1 = ibit1 >> 32;
587  l_ibit1 = ibit1 & i_32;
588  Sz += pop(u_ibit1);
589  Sz += pop(l_ibit1);
590 // isite1 <= X->Def.Nsite
591  ibit1 = (unsigned long int) tmp_list_1 &is1_up_b;
592  u_ibit1 = ibit1 >> 32;
593  l_ibit1 = ibit1 & i_32;
594  Sz += pop(u_ibit1);
595  Sz += pop(l_ibit1);
596  Sz = 2*Sz-X->Def.NsiteMPI;
597 
598  tmp_Sz += Sz*tmp_v02;
599  tmp_Sz2 += Sz*Sz*tmp_v02;
600  }
601  tmp_Sz = SumMPI_d(tmp_Sz);
602  tmp_Sz2 = SumMPI_d(tmp_Sz2);
603 
604  X->Phys.doublon = 0.0;
605  X->Phys.doublon2 = 0.0;
606  X->Phys.num = X->Def.NsiteMPI;
607  X->Phys.num2 = X->Def.NsiteMPI*X->Def.NsiteMPI;
608  X->Phys.Sz = tmp_Sz*0.5;
609  X->Phys.Sz2 = tmp_Sz2*0.25;
610  X->Phys.num_up = 0.5*(X->Def.NsiteMPI+tmp_Sz);
611  X->Phys.num_down = 0.5*(X->Def.NsiteMPI-tmp_Sz);
612 
613  return 0;
614 }
615 
621 int expec_energy_flct_GeneralSpin(struct BindStruct *X){
622  long unsigned int j;
623  long unsigned int isite1;
624 
625  double Sz,tmp_Sz, tmp_Sz2;
626  double tmp_v02;
627  long unsigned int i_max, tmp_list1;
628  i_max=X->Check.idim_max;
629 
630  // tentative doublon
631  tmp_Sz = 0.0;
632  tmp_Sz2 = 0.0;
633 
634 #pragma omp parallel for reduction(+:tmp_Sz,tmp_Sz2)default(none) shared(v0, list_1) \
635  firstprivate(i_max,X,myrank) private(j,Sz,isite1,tmp_v02, tmp_list1)
636  for(j = 1; j <= i_max; j++){
637  tmp_v02 = conj(v0[j])*v0[j];
638  Sz = 0.0;
639  tmp_list1 = list_1[j];
640  for(isite1=1;isite1<=X->Def.NsiteMPI;isite1++){
641  //prefactor 0.5 is added later.
642  if(isite1 > X->Def.Nsite){
643  Sz += GetLocal2Sz(isite1, myrank, X->Def.SiteToBit, X->Def.Tpow);
644  }else{
645  Sz += GetLocal2Sz(isite1, tmp_list1, X->Def.SiteToBit, X->Def.Tpow);
646  }
647  }
648  tmp_Sz += Sz*tmp_v02;
649  tmp_Sz2 += Sz*Sz*tmp_v02;
650  }
651 
652  tmp_Sz = SumMPI_d(tmp_Sz);
653  tmp_Sz2 = SumMPI_d(tmp_Sz2);
654 
655  X->Phys.doublon = 0.0;
656  X->Phys.doublon2 = 0.0;
657  X->Phys.num = X->Def.NsiteMPI;
658  X->Phys.num2 = X->Def.NsiteMPI*X->Def.NsiteMPI;
659  X->Phys.Sz = tmp_Sz*0.5;
660  X->Phys.Sz2 = tmp_Sz2*0.25;
661  X->Phys.num_up = 0.5*(X->Def.NsiteMPI+tmp_Sz);
662  X->Phys.num_down = 0.5*(X->Def.NsiteMPI-tmp_Sz);
663 
664  return 0;
665 }
cExpecEnd
const char * cExpecEnd
Definition: LogMessage.c:103
BindStruct::Def
struct DefineList Def
Definision of system (Hamiltonian) etc.
Definition: struct.h:410
mltply
int mltply(struct BindStruct *X, double complex *tmp_v0, double complex *tmp_v1)
Parent function of multiplying the wavefunction by the Hamiltonian. . First, the calculation of diago...
Definition: mltply.c:56
StartTimer
void StartTimer(int n)
function for initializing elapse time [start]
Definition: time.c:71
cTPQExpecEnd
const char * cTPQExpecEnd
Definition: LogMessage.c:105
SumMPI_dc
double complex SumMPI_dc(double complex norm)
MPI wrapper function to obtain sum of Double complex across processes.
Definition: wrapperMPI.c:205
LargeList::ihfbit
long unsigned int ihfbit
Used for Ogata-Lin ???
Definition: struct.h:345
CheckList::idim_max
unsigned long int idim_max
The dimension of the Hilbert space of this process.
Definition: struct.h:303
PhysList::num_up
double num_up
Expectation value of the number of up-spin electtrons.
Definition: struct.h:368
StopTimer
void StopTimer(int n)
function for calculating elapse time [elapse time=StartTimer-StopTimer]
Definition: time.c:83
v1
double complex * v1
Definition: global.h:35
cTPQExpecStart
const char * cTPQExpecStart
Definition: LogMessage.c:104
expec_energy_flct_GeneralSpin
int expec_energy_flct_GeneralSpin(struct BindStruct *X)
Calculate expected values of energies and physical quantities for General-Spin model.
Definition: expec_energy_flct.c:621
BindStruct::Large
struct LargeList Large
Variables for Matrix-Vector product.
Definition: struct.h:412
SumMPI_d
double SumMPI_d(double norm)
MPI wrapper function to obtain sum of Double across processes.
Definition: wrapperMPI.c:222
PhysList::num_down
double num_down
Expectation value of the number of down-spin electtrons.
Definition: struct.h:369
BindStruct::Phys
struct PhysList Phys
Physical quantities.
Definition: struct.h:413
expec_energy_flct_HubbardGC
int expec_energy_flct_HubbardGC(struct BindStruct *X)
Calculate expected values of energies and physical quantities for Hubbard GC model.
Definition: expec_energy_flct.c:185
PhysList::Sz2
double Sz2
Expectation value of the Square of total Sz.
Definition: struct.h:361
LargeList::mode
int mode
multiply or expectation value.
Definition: struct.h:330
GetLocal2Sz
int GetLocal2Sz(const unsigned int isite, const long unsigned int org_bit, const long int *SiteToBit, const long unsigned int *Tpow)
get 2sz at a site for general spin
Definition: bitcalc.c:449
PhysList::num
double num
Expectation value of the Number of electrons.
Definition: struct.h:358
BindStruct
Bind.
Definition: struct.h:409
DefineList::Nsite
unsigned int Nsite
Number of sites in the INTRA process region.
Definition: struct.h:56
pop
int pop(unsigned int x)
calculating number of 1-bits in x (32 bit) This method is introduced in S.H. Warren, Hacker$B!G(Bs Delight, second ed., Addison-Wesley, ISBN: 0321842685, 2012.
Definition: bitcalc.c:492
DefineList::Tpow
long unsigned int * Tpow
[2 * DefineList::NsiteMPI] malloc in setmem_def().
Definition: struct.h:90
LargeList::ilft
long unsigned int ilft
Used for Ogata-Lin ???
Definition: struct.h:344
GetSplitBitByModel
int GetSplitBitByModel(const int Nsite, const int iCalcModel, long unsigned int *irght, long unsigned int *ilft, long unsigned int *ihfbit)
function of splitting original bit into right and left spaces.
Definition: bitcalc.c:78
LargeList::i_max
long int i_max
Length of eigenvector.
Definition: struct.h:317
DefineList::NsiteMPI
unsigned int NsiteMPI
Total number of sites, differ from DefineList::Nsite.
Definition: struct.h:57
PhysList::var
double var
Expectation value of the Energy variance.
Definition: struct.h:363
DefineList::iFlgGeneralSpin
int iFlgGeneralSpin
Flag for the general (Sz/=1/2) spin.
Definition: struct.h:86
FALSE
#define FALSE
Definition: global.h:25
list_1
long unsigned int * list_1
Definition: global.h:47
expec_energy_flct_GeneralSpinGC
int expec_energy_flct_GeneralSpinGC(struct BindStruct *X)
Calculate expected values of energies and physical quantities for General-SpinGC model.
Definition: expec_energy_flct.c:500
cLogExpecEnergyStart
const char * cLogExpecEnergyStart
Definition: LogMessage.c:138
cExpecStart
const char * cExpecStart
Definition: LogMessage.c:102
DefineList::SiteToBit
long int * SiteToBit
[DefineList::NsiteMPI] Similar to DefineList::Tpow. For general spin.
Definition: struct.h:94
cFileNameTimeKeep
const char * cFileNameTimeKeep
Definition: global.c:23
expec_energy_flct_Hubbard
int expec_energy_flct_Hubbard(struct BindStruct *X)
Calculate expected values of energies and physical quantities for Hubbard model.
Definition: expec_energy_flct.c:305
expec_energy_flct_HalfSpin
int expec_energy_flct_HalfSpin(struct BindStruct *X)
Calculate expected values of energies and physical quantities for Half-Spin model.
Definition: expec_energy_flct.c:548
X
struct EDMainCalStruct X
Definition: struct.h:432
LargeList::irght
long unsigned int irght
Used for Ogata-Lin ???
Definition: struct.h:343
PhysList::doublon
double doublon
Expectation value of the Doublon.
Definition: struct.h:356
DefineList::iCalcModel
int iCalcModel
Switch for model. 0:Hubbard, 1:Spin, 2:Kondo, 3:HubbardGC, 4:SpinGC, 5:KondoGC, 6:HubbardNConserved.
Definition: struct.h:198
PhysList::Sz
double Sz
Expectation value of the Total Sz.
Definition: struct.h:360
myrank
int myrank
Process ID, defined in InitializeMPI()
Definition: global.h:163
expec_energy_flct
int expec_energy_flct(struct BindStruct *X)
Parent function to calculate expected values of energy and physical quantities.
Definition: expec_energy_flct.c:34
cLogExpecEnergyEnd
const char * cLogExpecEnergyEnd
Definition: LogMessage.c:139
PhysList::energy
double energy
Expectation value of the total energy.
Definition: struct.h:355
expec_energy_flct_HalfSpinGC
int expec_energy_flct_HalfSpinGC(struct BindStruct *X)
Calculate expected values of energies and physical quantities for Half-SpinGC model.
Definition: expec_energy_flct.c:428
v0
double complex * v0
Definition: global.h:34
BindStruct::Check
struct CheckList Check
Size of the Hilbert space.
Definition: struct.h:411
PhysList::num2
double num2
Expectation value of the quare of the number of electrons.
Definition: struct.h:359
step_i
int step_i
Definition: global.h:66
TimeKeeperWithStep
int TimeKeeperWithStep(struct BindStruct *X, const char *cFileName, const char *cTimeKeeper_Message, const char *cWriteType, const int istep)
Functions for writing a time log.
Definition: log.c:78
PhysList::doublon2
double doublon2
Expectation value of the Square of doublon.
Definition: struct.h:357
TimeKeeper
int TimeKeeper(struct BindStruct *X, const char *cFileName, const char *cTimeKeeper_Message, const char *cWriteType)
Functions for writing a time log.
Definition: log.c:42
DefineList::Total2SzMPI
int Total2SzMPI
Total across processes.
Definition: struct.h:70
DefineList::iCalcType
int iCalcType
Switch for calculation type. 0:Lanczos, 1:TPQCalc, 2:FullDiag.
Definition: struct.h:192
stdoutMPI
FILE * stdoutMPI
File pointer to the standard output defined in InitializeMPI()
Definition: global.h:165