18 #include "mltplyCommon.h" 20 #include "expec_energy_flct.h" 21 #include "wrapperMPI.h" 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;
41 switch(
X->Def.iCalcType){
60 i_max=
X->Check.idim_max;
65 X->Large.i_max = i_max;
66 X->Large.irght = irght;
68 X->Large.ihfbit = ihfbit;
69 X->Large.mode = M_ENERGY;
73 if(
X->Def.iCalcType == Lanczos){
76 else if (
X->Def.iCalcType == TPQCalc){
84 switch(
X->Def.iCalcModel){
95 if(
X->Def.iFlgGeneralSpin ==
FALSE) {
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;
125 #pragma omp parallel for default(none) private(i) shared(v1,v0) firstprivate(i_max) 126 for(i = 1; i <= i_max; i++){
132 if(
X->Def.iCalcType == Lanczos){
135 else if (
X->Def.iCalcType == TPQCalc){
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];
151 dam_pr1 += conj(
v0[j])*
v0[j];
158 X->Phys.energy = dam_pr;
159 X->Phys.var = dam_pr1;
161 switch(
X->Def.iCalcType) {
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;
196 long unsigned int i_max;
197 unsigned int l_ibit1, u_ibit1, i_32;
198 i_max=
X->Check.idim_max;
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];
219 is1_up_b +=
X->Def.Tpow[2 * isite1 - 2];
220 is1_down_b +=
X->Def.Tpow[2 * isite1 - 1];
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];
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);
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);
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);
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);
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);
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);
271 N = bit_up + bit_down;
272 Sz = bit_up - bit_down;
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;
288 X->Phys.doublon = tmp_D;
289 X->Phys.doublon2 = tmp_D2;
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);
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;
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;
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;
321 i_32 = (
unsigned int)(pow(2,32)-1);
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];
340 is1_up_b +=
X->Def.Tpow[2*isite1 - 2];
341 is1_down_b +=
X->Def.Tpow[2*isite1 - 1];
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];
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);
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);
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);
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);
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);
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);
393 N = bit_up + bit_down;
394 Sz = bit_up - bit_down;
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;
412 X->Phys.doublon = tmp_D;
413 X->Phys.doublon2 = tmp_D2;
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);
430 long unsigned int isite1;
431 long unsigned int is1_up_a,is1_up_b;
433 long unsigned int ibit1;
434 double Sz,tmp_Sz, tmp_Sz2;
436 long unsigned int i_max;
437 unsigned int l_ibit1,u_ibit1,i_32;
438 i_max=
X->Check.idim_max;
449 for(isite1=1;isite1<=
X->Def.NsiteMPI;isite1++){
450 if(isite1 >
X->Def.Nsite){
451 is1_up_a +=
X->Def.Tpow[isite1 - 1];
453 is1_up_b +=
X->Def.Tpow[isite1 - 1];
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];
464 ibit1 = (
unsigned long int)
myrank & is1_up_a;
465 u_ibit1 = ibit1 >> 32;
466 l_ibit1 = ibit1 & i_32;
470 ibit1 = (
unsigned long int) (j-1)&is1_up_b;
471 u_ibit1 = ibit1 >> 32;
472 l_ibit1 = ibit1 & i_32;
475 Sz = 2*Sz-
X->Def.NsiteMPI;
477 tmp_Sz += Sz*tmp_v02;
478 tmp_Sz2 += Sz*Sz*tmp_v02;
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);
502 long unsigned int isite1;
504 double Sz,tmp_Sz, tmp_Sz2;
506 long unsigned int i_max;
507 i_max=
X->Check.idim_max;
513 for(j = 1; j <= i_max; j++){
514 tmp_v02 = conj(
v0[j])*
v0[j];
516 for(isite1=1;isite1<=
X->Def.NsiteMPI;isite1++){
518 if(isite1 >
X->Def.Nsite){
521 Sz +=
GetLocal2Sz(isite1, j-1,
X->Def.SiteToBit,
X->Def.Tpow);
524 tmp_Sz += Sz*tmp_v02;
525 tmp_Sz2 += Sz*Sz*tmp_v02;
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);
550 long unsigned int isite1;
551 long unsigned int is1_up_a,is1_up_b;
553 long unsigned int ibit1;
554 double Sz,tmp_Sz, tmp_Sz2;
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;
569 for(isite1=1;isite1<=
X->Def.NsiteMPI;isite1++){
570 if(isite1 >
X->Def.Nsite){
571 is1_up_a +=
X->Def.Tpow[isite1 - 1];
573 is1_up_b +=
X->Def.Tpow[isite1 - 1];
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];
585 ibit1 = (
unsigned long int)
myrank & is1_up_a;
586 u_ibit1 = ibit1 >> 32;
587 l_ibit1 = ibit1 & i_32;
591 ibit1 = (
unsigned long int) tmp_list_1 &is1_up_b;
592 u_ibit1 = ibit1 >> 32;
593 l_ibit1 = ibit1 & i_32;
596 Sz = 2*Sz-
X->Def.NsiteMPI;
598 tmp_Sz += Sz*tmp_v02;
599 tmp_Sz2 += Sz*Sz*tmp_v02;
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);
623 long unsigned int isite1;
625 double Sz,tmp_Sz, tmp_Sz2;
627 long unsigned int i_max, tmp_list1;
628 i_max=
X->Check.idim_max;
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];
640 for(isite1=1;isite1<=
X->Def.NsiteMPI;isite1++){
642 if(isite1 >
X->Def.Nsite){
645 Sz +=
GetLocal2Sz(isite1, tmp_list1,
X->Def.SiteToBit,
X->Def.Tpow);
648 tmp_Sz += Sz*tmp_v02;
649 tmp_Sz2 += Sz*Sz*tmp_v02;
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);
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...
void StartTimer(int n)
function for initializing elapse time [start]
const char * cTPQExpecEnd
double complex SumMPI_dc(double complex norm)
MPI wrapper function to obtain sum of Double complex across processes.
void StopTimer(int n)
function for calculating elapse time [elapse time=StartTimer-StopTimer]
const char * cTPQExpecStart
int expec_energy_flct_GeneralSpin(struct BindStruct *X)
Calculate expected values of energies and physical quantities for General-Spin model.
double SumMPI_d(double norm)
MPI wrapper function to obtain sum of Double across processes.
int expec_energy_flct_HubbardGC(struct BindStruct *X)
Calculate expected values of energies and physical quantities for Hubbard GC model.
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
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.
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.
long unsigned int * list_1
int expec_energy_flct_GeneralSpinGC(struct BindStruct *X)
Calculate expected values of energies and physical quantities for General-SpinGC model.
const char * cLogExpecEnergyStart
const char * cFileNameTimeKeep
int expec_energy_flct_Hubbard(struct BindStruct *X)
Calculate expected values of energies and physical quantities for Hubbard model.
int expec_energy_flct_HalfSpin(struct BindStruct *X)
Calculate expected values of energies and physical quantities for Half-Spin model.
int myrank
Process ID, defined in InitializeMPI()
int expec_energy_flct(struct BindStruct *X)
Parent function to calculate expected values of energy and physical quantities.
const char * cLogExpecEnergyEnd
int expec_energy_flct_HalfSpinGC(struct BindStruct *X)
Calculate expected values of energies and physical quantities for Half-SpinGC model.
int TimeKeeperWithStep(struct BindStruct *X, const char *cFileName, const char *cTimeKeeper_Message, const char *cWriteType, const int istep)
Functions for writing a time log.
int TimeKeeper(struct BindStruct *X, const char *cFileName, const char *cTimeKeeper_Message, const char *cWriteType)
Functions for writing a time log.
FILE * stdoutMPI
File pointer to the standard output defined in InitializeMPI()