17 #include "CalcSpectrum.h" 18 #include "CalcSpectrumByLanczos.h" 19 #include "CalcSpectrumByBiCG.h" 20 #include "CalcSpectrumByTPQ.h" 21 #include "CalcSpectrumByFullDiag.h" 25 #include "wrapperMPI.h" 31 #include "diagonalcalc.h" 52 double complex *dcSpectrum,
53 double complex *dcomega)
65 for (i = 0; i < Nomega; i++) {
66 fprintf(fp,
"%.10lf %.10lf %.10lf %.10lf \n",
68 creal(dcSpectrum[i]), cimag(dcSpectrum[i]));
96 unsigned long int i_max = 0;
98 int iFlagListModified =
FALSE;
104 double complex OmegaMax, OmegaMin;
105 double complex *dcSpectrum;
106 double complex *dcomega;
111 fprintf(stderr,
"Error: Fail to set Omega.\n");
122 c_malloc1(dcSpectrum, Nomega);
123 c_malloc1(dcomega, Nomega);
126 for (i = 0; i < Nomega; i++) {
127 dcomega[i] = (OmegaMax - OmegaMin) / Nomega * i + OmegaMin;
129 fprintf(
stdoutMPI,
"\nFrequency range:\n");
130 fprintf(
stdoutMPI,
" Omega Max. : %15.5e %15.5e\n", creal(OmegaMax), cimag(OmegaMax));
131 fprintf(
stdoutMPI,
" Omega Min. : %15.5e %15.5e\n", creal(OmegaMin), cimag(OmegaMin));
132 fprintf(
stdoutMPI,
" Num. of Omega : %d\n", Nomega);
135 fprintf(stderr,
"Error: Any excitation operators are not defined.\n");
157 fprintf(
stdoutMPI,
" Start: An Eigenvector is inputted in CalcSpectrum.\n");
160 strcat(defname,
"_rank_%d.dat");
162 sprintf(sdt, defname,
myrank);
166 fprintf(stderr,
"Error: A file of Inputvector does not exist.\n");
170 byte_size = fread(&i_stp,
sizeof(i_stp), 1, fp);
172 byte_size = fread(&i_max,
sizeof(i_max), 1, fp);
174 fprintf(stderr,
"Error: myrank=%d, i_max=%ld\n",
myrank, i_max);
175 fprintf(stderr,
"Error: A file of Inputvector is incorrect.\n");
178 byte_size = fread(
v1Org,
sizeof(complex
double), i_max + 1, fp);
181 if (byte_size == 0) printf(
"byte_size: %d \n", (
int)byte_size);
186 fprintf(
stdoutMPI,
" End: An Inputcector is inputted in CalcSpectrum.\n\n");
189 fprintf(
stdoutMPI,
" Start: Calculating an excited Eigenvector.\n");
198 if (fabs(dnorm) < pow(10.0, -15)) {
199 fprintf(stderr,
"Warning: Norm of an excitation vector becomes 0.\n");
200 fprintf(
stdoutMPI,
" End: Calculating an excited Eigenvector.\n\n");
202 fprintf(
stdoutMPI,
" End: Calculating a spectrum.\n\n");
204 for (i = 0; i < Nomega; i++) {
211 #pragma omp parallel for default(none) private(i) shared(v1, v0) firstprivate(i_max, dnorm, X) 213 v1[i] =
v0[i] / dnorm;
216 fprintf(
stdoutMPI,
" End: Calculating an excited Eigenvector.\n\n");
221 if (iFlagListModified ==
TRUE) {
232 fprintf(
stdoutMPI,
" Start: Calculating a spectrum.\n\n");
259 fprintf(stderr,
" Error: TPQ is not supported for calculating spectrum mode.\n");
280 fprintf(stderr,
" Error: The selected calculation type is not supported for calculating spectrum mode.\n");
284 fprintf(
stdoutMPI,
" End: Calculating a spectrum.\n\n");
290 c_free1(dcSpectrum, Nomega);
291 c_free1(dcomega, Nomega);
305 double complex *tmp_v0,
306 double complex *tmp_v1
309 if(
X->Def.NSingleExcitationOperator > 0 &&
X->Def.NPairExcitationOperator > 0){
310 fprintf(stderr,
"Error: Both single and pair excitation operators exist.\n");
315 if(
X->Def.NSingleExcitationOperator > 0){
320 else if(
X->Def.NPairExcitationOperator >0){
343 double E1, E2, E3, E4, Emax;
344 long unsigned int iline_countMax=2;
345 long unsigned int iline_count=2;
348 if(
X->iFlgSpecOmegaMax ==
TRUE &&
X->iFlgSpecOmegaMin ==
TRUE){
352 if (
X->iCalcType == Lanczos ||
X->iCalcType == FullDiag) {
356 fprintf(
stdoutMPI,
"Error: xx_Lanczos_Step.dat does not exist.\n");
361 while (
fgetsMPI(ctmp, 256, fp) != NULL) {
364 iline_countMax = iline_count;
370 while (
fgetsMPI(ctmp, 256, fp) != NULL) {
371 sscanf(ctmp,
"stp=%d %lf %lf %lf %lf %lf\n",
379 if (iline_count == iline_countMax)
break;
383 fprintf(
stdoutMPI,
"Error: Lanczos step must be greater than 4 for using spectrum calculation.\n");
392 fprintf(
stdoutMPI,
"Error: xx_energy.dat does not exist.\n");
397 sscanf(ctmp,
" Energy %lf \n", &E1);
401 if(
X->iFlgSpecOmegaMax ==
FALSE){
402 X->dcOmegaMax= Emax*(double)
X->Nsite;
404 if(
X->iFlgSpecOmegaMin ==
FALSE){
426 *iFlgListModifed =
FALSE;
433 X->Check.idim_maxOrg =
X->Check.idim_max;
434 X->Check.idim_maxMPIOrg =
X->Check.idim_maxMPI;
436 if (
X->Def.NSingleExcitationOperator > 0) {
437 switch (
X->Def.iCalcModel) {
440 case HubbardNConserved:
444 *iFlgListModifed =
TRUE;
450 }
else if (
X->Def.NPairExcitationOperator > 0) {
451 switch (
X->Def.iCalcModel) {
454 case HubbardNConserved:
460 if (
X->Def.PairExcitationOperator[0][1] !=
X->Def.PairExcitationOperator[0][3]) {
461 *iFlgListModifed =
TRUE;
469 if (*iFlgListModifed ==
TRUE) {
484 for(j =0; j<
X->Large.SizeOflist_2_1; j++){
487 for(j =0; j<
X->Large.SizeOflist_2_2; j++){
497 if (
X->Def.NSingleExcitationOperator > 0) {
498 switch (
X->Def.iCalcModel) {
501 case HubbardNConserved:
502 if (
X->Def.SingleExcitationOperator[0][2] == 1) {
503 X->Def.Ne =
X->Def.NeMPI + 1;
506 X->Def.Ne =
X->Def.NeMPI - 1;
512 if (
X->Def.SingleExcitationOperator[0][2] == 1) {
513 X->Def.Ne =
X->Def.NeMPI + 1;
514 if (
X->Def.SingleExcitationOperator[0][1] == 0) {
515 X->Def.Nup =
X->Def.NupOrg + 1;
516 X->Def.Ndown=
X->Def.NdownOrg;
518 X->Def.Nup=
X->Def.NupOrg;
519 X->Def.Ndown =
X->Def.NdownOrg + 1;
522 X->Def.Ne =
X->Def.NeMPI - 1;
523 if (
X->Def.SingleExcitationOperator[0][1] == 0) {
524 X->Def.Nup =
X->Def.NupOrg - 1;
525 X->Def.Ndown=
X->Def.NdownOrg;
528 X->Def.Nup=
X->Def.NupOrg;
529 X->Def.Ndown =
X->Def.NdownOrg - 1;
537 }
else if (
X->Def.NPairExcitationOperator > 0) {
538 X->Def.Ne=
X->Def.NeMPI;
539 switch (
X->Def.iCalcModel) {
542 case HubbardNConserved:
547 if (
X->Def.PairExcitationOperator[0][1] !=
X->Def.PairExcitationOperator[0][3]) {
548 if (
X->Def.PairExcitationOperator[0][1] == 0) {
549 X->Def.Nup =
X->Def.NupOrg + 1;
550 X->Def.Ndown =
X->Def.NdownOrg - 1;
552 X->Def.Nup =
X->Def.NupOrg - 1;
553 X->Def.Ndown =
X->Def.NdownOrg + 1;
558 if (
X->Def.PairExcitationOperator[0][1] !=
X->Def.PairExcitationOperator[0][3]) {
559 if (
X->Def.iFlgGeneralSpin ==
FALSE) {
560 if (
X->Def.PairExcitationOperator[0][1] == 0) {
561 X->Def.Nup =
X->Def.NupOrg - 1;
562 X->Def.Ndown =
X->Def.NdownOrg + 1;
564 X->Def.Nup =
X->Def.NupOrg + 1;
565 X->Def.Ndown =
X->Def.NdownOrg - 1;
569 X->Def.Total2Sz =
X->Def.Total2SzMPI+2*(
X->Def.PairExcitationOperator[0][1]-
X->Def.PairExcitationOperator[0][3]);
578 X->Def.Nsite=
X->Def.NsiteMPI;
596 if(
X->Def.iCalcModel==HubbardNConserved){
597 X->Def.iCalcModel=Hubbard;
601 if (*iFlgListModifed ==
TRUE) {
602 for(j=1; j<=
X->Check.idim_maxOrg; j++){
603 fprintf(stdout,
"Debug1: myrank=%d, list_1_org[ %ld] = %ld\n",
myrank, j,
list_1_org[j]+
myrank*
X->Def.OrgTpow[2*
X->Def.NsiteMPI-1]);
606 for(j=1; j<=
X->Check.idim_max; j++){
607 fprintf(stdout,
"Debug2: myrank=%d, list_1[ %ld] = %ld\n",
myrank, j,
list_1[j]+
myrank* 64);
unsigned int NSingleExcitationOperator
Number of single excitaion operator for spectrum.
int childfopenALL(const char *_cPathChild, const char *_cmode, FILE **_fp)
All processes open file in output/ directory.
void exitMPI(int errorcode)
MPI Abortation wrapper.
struct DefineList Def
Definision of system (Hamiltonian) etc.
int MakeExcitedList(struct BindStruct *X, int *iFlgListModifed)
Make the lists for the excited state; list_1, list_2_1 and list_2_2 (for canonical ensemble)...
void StartTimer(int n)
function for initializing elapse time [start]
long unsigned int * list_2_1_org
unsigned long int idim_max
The dimension of the Hilbert space of this process.
int CalcSpectrumByFullDiag(struct EDMainCalStruct *X, int Nomega, double complex *dcSpectrum, double complex *dcomega)
Compute the Green function with the Lehmann representation and FD .
int sz(struct BindStruct *X, long unsigned int *list_1_, long unsigned int *list_2_1_, long unsigned int *list_2_2_)
generating Hilbert space
void StopTimer(int n)
function for calculating elapse time [elapse time=StartTimer-StopTimer]
const char * c_InputEigenVectorEnd
int iFlagListModified
When the Hilbert space of excited state differs from the original one.
double complex dcOmegaMax
Upper limit of the frequency for the spectrum.
const char * c_CalcSpectrumEnd
struct LargeList Large
Variables for Matrix-Vector product.
int childfopenMPI(const char *_cPathChild, const char *_cmode, FILE **_fp)
Only the root process open file in output/ directory.
const char * c_CalcExcitedStateEnd
const char * cFileNameLanczosStep
int check(struct BindStruct *X)
A program to check size of dimension for Hilbert-space.
int setmem_large(struct BindStruct *X)
Set size of memories for Hamiltonian (Ham, L_vec), vectors(vg, v0, v1, v2, vec, alpha, beta), lists (list_1, list_2_1, list_2_2, list_Diagonal) and Phys(BindStruct.PhysList) struct in the case of Full Diag mode.
long unsigned int * list_1buf_org
int GetPairExcitedState(struct BindStruct *X, double complex *tmp_v0, double complex *tmp_v1)
Calculating the pair excited state by the pair operator; , where indicates a creation (anti-creat...
const char * c_CalcSpectrumStart
int GetExcitedState(struct BindStruct *X, double complex *tmp_v0, double complex *tmp_v1)
Parent function to calculate the excited state.
double NormMPI_dc(unsigned long int idim, double complex *_v1)
Compute norm of process-distributed vector .
int iErrCodeMem
Error Message in HPhiMain.c.
long unsigned int * list_2_2_org
const char * cFileNameCalcDynamicalGreen
int CalcSpectrumByBiCG(struct EDMainCalStruct *X, double complex *vrhs, double complex *v2, double complex *v4, int Nomega, double complex *dcSpectrum, double complex *dcomega)
A main function to calculate spectrum by BiCG method In this function, the library is used...
const char * c_InputEigenVectorStart
int diagonalcalc(struct BindStruct *X)
Calculate diagonal components and obtain the list, list_diagonal.
int CalcSpectrumByLanczos(struct EDMainCalStruct *X, double complex *tmp_v1, double dnorm, int Nomega, double complex *dcSpectrum, double complex *dcomega)
A main function to calculate spectrum by continued fraction expansions.
unsigned long int idim_maxOrg
The local Hilbert-space dimention of original state for the spectrum.
int CalcSpectrumByTPQ(struct EDMainCalStruct *X, double complex *tmp_v1, double dnorm, int Nomega, double complex *dcSpectrum, double complex *dcomega)
A main function to calculate spectrum by TPQ (Note: This method is trial)
const char * c_CalcExcitedStateStart
long unsigned int * list_1_org
long unsigned int * list_2_1
int GetlistSize(struct BindStruct *X)
Set size of lists for the canonical ensemble.
void FinalizeMPI()
MPI Finitialization wrapper.
int CalcSpectrum(struct EDMainCalStruct *X)
A main function to calculate spectrum.
int SetOmega(struct DefineList *X)
Set target frequencies.
int iNOmega
Number of frequencies for spectrum.
long unsigned int * list_1
int iFlgSpecOmegaOrg
Whether DefineList::dcOmegaOrg is input or not.
unsigned int NPairExcitationOperator
Number of pair excitaion operator for spectrum.
const char * cFileNameTimeKeep
long unsigned int * list_2_2
double complex dcOmegaOrg
Origin limit of the frequency for the spectrum.
double complex dcOmegaMin
Lower limit of the frequency for the spectrum.
int GetSingleExcitedState(struct BindStruct *X, double complex *tmp_v0, double complex *tmp_v1)
Calculation of single excited state Target System: Hubbard, Kondo.
int myrank
Process ID, defined in InitializeMPI()
int iFlgCalcSpec
Input parameter CalcSpec in teh CalcMod file.
int OutputSpectrum(struct EDMainCalStruct *X, int Nomega, double complex *dcSpectrum, double complex *dcomega)
Output spectrum.
int GetFileNameByKW(int iKWidx, char **FileName)
function of getting file name labeled by the keyword
char * fgetsMPI(char *InputString, int maxcount, FILE *fp)
MPI file I/O (get a line, fgets) wrapper. Only the root node (myrank = 0) reads and broadcast string...
const char * cFileNameEnergy_Lanczos
struct CheckList Check
Size of the Hilbert space.
char * CDataFileHead
Read from Calcmod in readdef.h. Header of output file such as Green's function.
struct BindStruct Bind
Binded struct.
Definision of system (Hamiltonian) etc.
int TimeKeeper(struct BindStruct *X, const char *cFileName, const char *cTimeKeeper_Message, const char *cWriteType)
Functions for writing a time log.
int iCalcType
Switch for calculation type. 0:Lanczos, 1:TPQCalc, 2:FullDiag.
FILE * stdoutMPI
File pointer to the standard output defined in InitializeMPI()