HΦ  3.1.0
Functions
CalcSpectrumByTPQ.c File Reference

Calculate spectrum function for the TPQ state.
Note: This method is trial and cannot be used in the release mode. More...

#include "CalcSpectrumByTPQ.h"
#include "Lanczos_EigenValue.h"
#include "FileIO.h"
#include "wrapperMPI.h"
#include "vec12.h"
#include "mfmemory.h"

Go to the source code of this file.

Functions

int ReadTPQData (struct EDMainCalStruct *X, double *ene, double *temp, double *specificHeat)
 Read TPQ data at "X->Bind.Large.itr" step in SS_rand file. More...
 
int GetCalcSpectrumTPQ (double complex dcomega, double dtemp, double dspecificheat, double ene, double *tmp_E, int Nsite, int idim_max, double complex *dc_tmpSpec)
 Calculate spectrum function from the TPQ state. More...
 
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) More...
 

Detailed Description

Calculate spectrum function for the TPQ state.
Note: This method is trial and cannot be used in the release mode.

Version
1.2
Author
Kazuyoshi Yoshimi (The University of Tokyo)

Definition in file CalcSpectrumByTPQ.c.

Function Documentation

◆ CalcSpectrumByTPQ()

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)

Parameters
X[in,out] CalcStruct list for getting and pushing calculation information
tmp_v1[in] Normalized excited state.
dnorm[in] Norm of the excited state before normalization.
Nomega[in] Total number of frequencies.
dcSpectrum[out] Calculated spectrum.
dcomega[in] Target frequencies.
Return values
0normally finished
-1unnormally finished
Version
1.2
Author
Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 131 of file CalcSpectrumByTPQ.c.

References alpha, beta, EDMainCalStruct::Bind, c_CalcSpectrumFromTridiagonalEnd, c_CalcSpectrumFromTridiagonalStart, c_GetTridiagonalEnd, c_GetTridiagonalStart, c_InputSpectrumRecalcvecEnd, c_InputSpectrumRecalcvecStart, c_OutputSpectrumRecalcvecEnd, c_OutputSpectrumRecalcvecStart, DefineList::CDataFileHead, cFileNameOutputRestartVec, cFileNameTimeKeep, BindStruct::Check, childfopenALL(), D_FileNameMax, BindStruct::Def, exitMPI(), FALSE, GetCalcSpectrumTPQ(), CheckList::idim_max, DefineList::iFlgCalcSpec, Lanczos_GetTridiagonalMatrixComponents(), DefineList::Lanczos_max, DefineList::Lanczos_restart, myrank, DefineList::Nsite, DefineList::nvec, OutputTMComponents(), ReadTMComponents(), ReadTPQData(), stdoutMPI, TimeKeeper(), TRUE, v0, v1, vec12(), and X.

Referenced by CalcSpectrum().

139 {
140  char sdt[D_FileNameMax];
141  unsigned long int i, i_max;
142  FILE *fp;
143  int iret;
144  unsigned long int liLanczosStp = X->Bind.Def.Lanczos_max;
145  unsigned long int liLanczosStp_vec=0;
146  double dene, dtemp, dspecificHeat;
147  double *tmp_E;
148  double complex dctmp_Spectrum;
149  int stp;
150  size_t byte_size;
151 
152  //Read Ene, temp, C
153  if(ReadTPQData(X, &dene, &dtemp, &dspecificHeat)!=TRUE){
154  return FALSE;
155  }
156 
157  if(X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents_VEC ||
158  X->Bind.Def.iFlgCalcSpec == RECALC_INOUT_TMComponents_VEC) {
159  fprintf(stdoutMPI, " Start: Input vectors for recalculation.\n");
160  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_InputSpectrumRecalcvecStart, "a");
161 
162  sprintf(sdt, cFileNameOutputRestartVec, X->Bind.Def.CDataFileHead, myrank);
163  if (childfopenALL(sdt, "rb", &fp) != 0) {
164  exitMPI(-1);
165  }
166  byte_size = fread(&liLanczosStp_vec, sizeof(liLanczosStp_vec),1,fp);
167  byte_size = fread(&i_max, sizeof(long int), 1, fp);
168  if(i_max != X->Bind.Check.idim_max){
169  fprintf(stderr, "Error: A size of Inputvector is incorrect.\n");
170  exitMPI(-1);
171  }
172  byte_size = fread(v0, sizeof(complex double), X->Bind.Check.idim_max + 1, fp);
173  byte_size = fread(v1, sizeof(complex double), X->Bind.Check.idim_max + 1, fp);
174  fclose(fp);
175  fprintf(stdoutMPI, " End: Input vectors for recalculation.\n");
176  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_InputSpectrumRecalcvecEnd, "a");
177  if (byte_size == 0) printf("byte_size: %d \n", (int)byte_size);
178  }
179 
180  //Read diagonal components
181  if(X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents ||
182  X->Bind.Def.iFlgCalcSpec ==RECALC_FROM_TMComponents_VEC||
183  X->Bind.Def.iFlgCalcSpec == RECALC_INOUT_TMComponents_VEC)
184  {
185  int iFlgTMComp=0;
186  if(X->Bind.Def.iFlgCalcSpec == RECALC_INOUT_TMComponents_VEC ||
187  X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents_VEC)
188  {
189  iFlgTMComp=0;
190  }
191  else{
192  iFlgTMComp=1;
193  }
194  iret=ReadTMComponents(&(X->Bind), &dnorm, &liLanczosStp, iFlgTMComp);
195  if(iret !=TRUE){
196  fprintf(stdoutMPI, " Error: Fail to read TMcomponents\n");
197  return FALSE;
198  }
199 
200  if(X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents){
201  X->Bind.Def.Lanczos_restart=0;
202  }
203  else if(X->Bind.Def.iFlgCalcSpec == RECALC_INOUT_TMComponents_VEC||
204  X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents_VEC){
205  if(liLanczosStp_vec !=liLanczosStp){
206  fprintf(stdoutMPI, " Error: Input files for vector and TMcomponents are incoorect.\n");
207  fprintf(stdoutMPI, " Error: Input vector %ld th stps, TMcomponents %ld th stps.\n", liLanczosStp_vec, liLanczosStp);
208  return FALSE;
209  }
210  X->Bind.Def.Lanczos_restart=liLanczosStp;
211  liLanczosStp = liLanczosStp+X->Bind.Def.Lanczos_max;
212  }
213  }
214 
215  // calculate ai, bi
216  if (X->Bind.Def.iFlgCalcSpec == RECALC_NOT ||
217  X->Bind.Def.iFlgCalcSpec == RECALC_OUTPUT_TMComponents_VEC ||
218  X->Bind.Def.iFlgCalcSpec == RECALC_FROM_TMComponents_VEC ||
219  X->Bind.Def.iFlgCalcSpec == RECALC_INOUT_TMComponents_VEC
220  )
221  {
222  fprintf(stdoutMPI, " Start: Calculate tridiagonal matrix components.\n");
223  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_GetTridiagonalStart, "a");
224  // Functions in Lanczos_EigenValue
225  iret = Lanczos_GetTridiagonalMatrixComponents(&(X->Bind), alpha, beta, tmp_v1, &(liLanczosStp));
226  if (iret != TRUE) {
227  //Error Message will be added.
228  return FALSE;
229  }
230  fprintf(stdoutMPI, " End: Calculate tridiagonal matrix components.\n\n");
231  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_GetTridiagonalEnd, "a");
232  OutputTMComponents(&(X->Bind), alpha,beta, dnorm, liLanczosStp);
233  }//X->Bind.Def.iFlgCalcSpec == RECALC_NOT || RECALC_FROM_TMComponents_VEC
234 
235  stp=liLanczosStp;
236  d_malloc1(tmp_E,stp+1);
237  X->Bind.Def.nvec= stp;
238  vec12(alpha,beta,stp,tmp_E, &(X->Bind));
239  fprintf(stdoutMPI, " Start: Caclulate spectrum from tridiagonal matrix components.\n");
240  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_CalcSpectrumFromTridiagonalStart, "a");
241  for( i = 0 ; i < Nomega; i++) {
242  dctmp_Spectrum=0;
243  iret = GetCalcSpectrumTPQ(dcomega[i], dtemp, dspecificHeat, dene, tmp_E, X->Bind.Def.Nsite, stp, &dctmp_Spectrum);
244  if (iret != TRUE) {
245  //ReAlloc alpha, beta and Set alpha_start and beta_start in Lanczos_EigenValue
246  return FALSE;
247  }
248  dcSpectrum[i] = dnorm * dctmp_Spectrum;
249  }
250  fprintf(stdoutMPI, " End: Caclulate spectrum from tridiagonal matrix components.\n\n");
251  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_CalcSpectrumFromTridiagonalEnd, "a");
252 
253  d_free1(tmp_E,stp+1);
254  //output vectors for recalculation
255  if(X->Bind.Def.iFlgCalcSpec==RECALC_OUTPUT_TMComponents_VEC ||
256  X->Bind.Def.iFlgCalcSpec==RECALC_INOUT_TMComponents_VEC){
257  fprintf(stdoutMPI, " Start: Output vectors for recalculation.\n");
258  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_OutputSpectrumRecalcvecStart, "a");
259 
260  sprintf(sdt, cFileNameOutputRestartVec, X->Bind.Def.CDataFileHead, myrank);
261  if(childfopenALL(sdt, "wb", &fp)!=0){
262  exitMPI(-1);
263  }
264  fwrite(&liLanczosStp, sizeof(liLanczosStp),1,fp);
265  fwrite(&X->Bind.Check.idim_max, sizeof(X->Bind.Check.idim_max),1,fp);
266  fwrite(v0, sizeof(complex double),X->Bind.Check.idim_max+1, fp);
267  fwrite(v1, sizeof(complex double),X->Bind.Check.idim_max+1, fp);
268  fclose(fp);
269 
270  fprintf(stdoutMPI, " End: Output vectors for recalculation.\n");
271  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_OutputSpectrumRecalcvecEnd, "a");
272  }
273 
274  return TRUE;
275 }
childfopenALL
int childfopenALL(const char *_cPathChild, const char *_cmode, FILE **_fp)
All processes open file in output/ directory.
Definition: FileIO.c:50
exitMPI
void exitMPI(int errorcode)
MPI Abortation wrapper.
Definition: wrapperMPI.c:86
BindStruct::Def
struct DefineList Def
Definision of system (Hamiltonian) etc.
Definition: struct.h:409
CheckList::idim_max
unsigned long int idim_max
The dimension of the Hilbert space of this process.
Definition: struct.h:302
v1
double complex * v1
Definition: global.h:35
ReadTPQData
int ReadTPQData(struct EDMainCalStruct *X, double *ene, double *temp, double *specificHeat)
Read TPQ data at "X->Bind.Large.itr" step in SS_rand file.
Definition: CalcSpectrumByTPQ.c:41
D_FileNameMax
#define D_FileNameMax
Definition: global.h:23
Lanczos_GetTridiagonalMatrixComponents
int Lanczos_GetTridiagonalMatrixComponents(struct BindStruct *X, double *_alpha, double *_beta, double complex *tmp_v1, unsigned long int *liLanczos_step)
Calculate tridiagonal matrix components by Lanczos method.
Definition: Lanczos_EigenValue.c:321
cFileNameOutputRestartVec
const char * cFileNameOutputRestartVec
Definition: global.c:81
DefineList::nvec
unsigned int nvec
Read from Calcmod in readdef.h.
Definition: struct.h:46
TRUE
#define TRUE
Definition: global.h:26
c_OutputSpectrumRecalcvecEnd
const char * c_OutputSpectrumRecalcvecEnd
Definition: LogMessage.c:71
DefineList::Nsite
unsigned int Nsite
Number of sites in the INTRA process region.
Definition: struct.h:56
DefineList::Lanczos_restart
int Lanczos_restart
Number of iterations performed in the restart computation.
Definition: struct.h:75
c_InputSpectrumRecalcvecStart
const char * c_InputSpectrumRecalcvecStart
Definition: LogMessage.c:72
c_GetTridiagonalStart
const char * c_GetTridiagonalStart
Definition: LogMessage.c:66
c_CalcSpectrumFromTridiagonalEnd
const char * c_CalcSpectrumFromTridiagonalEnd
Definition: LogMessage.c:69
alpha
double * alpha
Definition: global.h:44
FALSE
#define FALSE
Definition: global.h:25
GetCalcSpectrumTPQ
int GetCalcSpectrumTPQ(double complex dcomega, double dtemp, double dspecificheat, double ene, double *tmp_E, int Nsite, int idim_max, double complex *dc_tmpSpec)
Calculate spectrum function from the TPQ state.
Definition: CalcSpectrumByTPQ.c:93
cFileNameTimeKeep
const char * cFileNameTimeKeep
Definition: global.c:23
c_InputSpectrumRecalcvecEnd
const char * c_InputSpectrumRecalcvecEnd
Definition: LogMessage.c:73
beta
double * beta
Definition: global.h:44
vec12
void vec12(double alpha[], double beta[], unsigned int ndim, double tmp_E[], struct BindStruct *X)
Diagonalize a tri-diagonal matrix and store eigenvectors into vec.
Definition: vec12.c:31
OutputTMComponents
int OutputTMComponents(struct BindStruct *X, double *_alpha, double *_beta, double _dnorm, int liLanczosStp)
Output tridiagonal matrix components obtained by the Lanczos method.
Definition: Lanczos_EigenValue.c:691
X
struct EDMainCalStruct X
Definition: struct.h:431
myrank
int myrank
Process ID, defined in InitializeMPI()
Definition: global.h:162
DefineList::iFlgCalcSpec
int iFlgCalcSpec
Input parameter CalcSpec in teh CalcMod file.
Definition: struct.h:215
c_CalcSpectrumFromTridiagonalStart
const char * c_CalcSpectrumFromTridiagonalStart
Definition: LogMessage.c:68
c_GetTridiagonalEnd
const char * c_GetTridiagonalEnd
Definition: LogMessage.c:67
ReadTMComponents
int ReadTMComponents(struct BindStruct *X, double *_dnorm, unsigned long int *_i_max, int iFlg)
Read tridiagonal matrix components obtained by the Lanczos method. .
Definition: Lanczos_EigenValue.c:621
v0
double complex * v0
Definition: global.h:34
c_OutputSpectrumRecalcvecStart
const char * c_OutputSpectrumRecalcvecStart
Definition: LogMessage.c:70
BindStruct::Check
struct CheckList Check
Size of the Hilbert space.
Definition: struct.h:410
DefineList::CDataFileHead
char * CDataFileHead
Read from Calcmod in readdef.h. Header of output file such as Green&#39;s function.
Definition: struct.h:42
EDMainCalStruct::Bind
struct BindStruct Bind
Binded struct.
Definition: struct.h:419
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::Lanczos_max
unsigned int Lanczos_max
Maximum number of iterations.
Definition: struct.h:74
stdoutMPI
FILE * stdoutMPI
File pointer to the standard output defined in InitializeMPI()
Definition: global.h:164

◆ GetCalcSpectrumTPQ()

int GetCalcSpectrumTPQ ( double complex  dcomega,
double  dtemp,
double  dspecificheat,
double  ene,
double *  tmp_E,
int  Nsite,
int  idim_max,
double complex *  dc_tmpSpec 
)

Calculate spectrum function from the TPQ state.

Parameters
dcomega[in] Target frequencies.
dtemp[in] Temperature corresponding to the target TPQ state.
dspecificheat[in] Specific heat.
ene[in] Energy for the target TPQ state.
tmp_E[in] Energies included in the excited TPQ state obtained by the continued fraction expansions.
Nsite[in] Total number of sites.
idim_max[in] Dimension of the Hilbert space.
dc_tmpSpec[out] Calculated spectrum.
Return values
FALSEfail to calculate spectrum.
TRUEsucsceed to calculate spectrum.

Definition at line 93 of file CalcSpectrumByTPQ.c.

References FALSE, TRUE, and vec.

Referenced by CalcSpectrumByTPQ().

95 {
96  int l;
97  double tmp_dcSpec;
98  double factor, pre_factor;
99  pre_factor=2.0*dtemp*dtemp*dspecificheat;
100  factor=M_PI*pre_factor;
101  factor=1.0/sqrt(factor);
102  tmp_dcSpec=0;
103 
104  if(cimag(dcomega)>0) {
105  for (l = 1; l <= idim_max; l++) {
106  //TODO: Check omega is real ?
107  //fprintf(stdoutMPI, "Debug: %lf, %lf\n", creal(dcomega) - tmp_E[l] + ene, pre_factor);
108  tmp_dcSpec += (double)(vec[l][1] * conj(vec[l][1])) * exp(-pow((creal(dcomega) - tmp_E[l] + ene),2)/(pre_factor));
109  }
110  }
111  else{
112  fprintf(stderr, " an imarginary part of omega must be positive.\n");
113  return FALSE;
114  }
115  tmp_dcSpec *=factor;
116  *dc_tmpSpec = tmp_dcSpec;
117  return TRUE;
118 }
vec
double complex ** vec
Definition: global.h:45
TRUE
#define TRUE
Definition: global.h:26
FALSE
#define FALSE
Definition: global.h:25

◆ ReadTPQData()

int ReadTPQData ( struct EDMainCalStruct *  X,
double *  ene,
double *  temp,
double *  specificHeat 
)

Read TPQ data at "X->Bind.Large.itr" step in SS_rand file.

Parameters
[in]XCalcStruct list for getting and pushing calculation information
[out]eneenergy
[out]temptemperature
[out]specificHeatspecific heat
Return values
TRUEsucceed to read data
FALSEfail to read data

Definition at line 41 of file CalcSpectrumByTPQ.c.

Referenced by CalcSpectrumByTPQ().

46  {
47  FILE *fp;
48  char sdt[D_FileNameMax];
49  char ctmp2[256];
50  double dinv_temp;
51  double dene, dHvar, dn, ddoublon;
52  int istp;
53  sprintf(sdt, cFileNameSSRand, X->Bind.Def.irand);
54  childfopenMPI(sdt, "r", &fp);
55  if(fp==NULL){
56  fprintf(stderr, " Error: SS_rand%d.dat does not exist.\n", X->Bind.Def.irand);
57  fclose(fp);
58  return FALSE;
59  }
60  fgetsMPI(ctmp2, 256, fp);
61  while(fgetsMPI(ctmp2, 256, fp) != NULL) {
62  sscanf(ctmp2, "%lf %lf %lf %lf %lf %d\n",
63  &dinv_temp,
64  &dene,
65  &dHvar,
66  &dn,
67  &ddoublon,
68  &istp
69  );
70  if(istp==X->Bind.Large.itr) break;
71  }
72  fclose(fp);
73 
74  *ene = dene;
75  *temp = 1.0/dinv_temp;
76  *specificHeat = (dHvar-dene*dene)*(dinv_temp*dinv_temp);
77 
78  return TRUE;
79 }
DefineList::irand
int irand
Input keyword TargetTPQRand ???
Definition: struct.h:79
BindStruct::Def
struct DefineList Def
Definision of system (Hamiltonian) etc.
Definition: struct.h:409
LargeList::itr
int itr
Iteration number.
Definition: struct.h:314
D_FileNameMax
#define D_FileNameMax
Definition: global.h:23
BindStruct::Large
struct LargeList Large
Variables for Matrix-Vector product.
Definition: struct.h:411
childfopenMPI
int childfopenMPI(const char *_cPathChild, const char *_cmode, FILE **_fp)
Only the root process open file in output/ directory.
Definition: FileIO.c:27
TRUE
#define TRUE
Definition: global.h:26
FALSE
#define FALSE
Definition: global.h:25
X
struct EDMainCalStruct X
Definition: struct.h:431
cFileNameSSRand
const char * cFileNameSSRand
Definition: global.c:56
fgetsMPI
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...
Definition: wrapperMPI.c:122
EDMainCalStruct::Bind
struct BindStruct Bind
Binded struct.
Definition: struct.h:419