CalcSpectrum.c File Reference

File for givinvg functions of calculating spectrum. More...

#include "mltply.h"
#include "CalcSpectrum.h"
#include "CalcSpectrumByLanczos.h"
#include "CalcSpectrumByBiCG.h"
#include "CalcSpectrumByTPQ.h"
#include "CalcSpectrumByFullDiag.h"
#include "CalcTime.h"
#include "SingleEx.h"
#include "PairEx.h"
#include "wrapperMPI.h"
#include "FileIO.h"
#include "mfmemory.h"
#include "readdef.h"
#include "sz.h"
#include "check.h"
#include "diagonalcalc.h"

Go to the source code of this file.

Functions
int	OutputSpectrum (struct EDMainCalStruct *X, int Nomega, double complex *dcSpectrum, double complex *dcomega)
	Output spectrum. More...
int	CalcSpectrum (struct EDMainCalStruct *X)
	A main function to calculate spectrum. More...
int	GetExcitedState (struct BindStruct *X, double complex *tmp_v0, double complex *tmp_v1)
	Parent function to calculate the excited state. More...
int	SetOmega (struct DefineList *X)
	Set target frequencies. More...
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). The original lists before the excitation are given by list_xxx_org. More...

Detailed Description

File for givinvg functions of calculating spectrum.

Version

1.1

Author

Kazuyoshi Yoshimi (The University of Tokyo)

Definition in file CalcSpectrum.c.

Function Documentation

CalcSpectrum()

int CalcSpectrum

(

struct EDMainCalStruct *

X

)

A main function to calculate spectrum.

Parameters

X	[in,out] CalcStruct list for getting and pushing calculation information input: iFlgSpecOmegaOrg, dcOmegaMax, dcOmegaMin, iNOmega etc. output: dcOmegaOrg, iFlagListModified.

Return values

0	normally finished
-1	unnormally finished

Version

1.1

Author

Kazuyoshi Yoshimi (The University of Tokyo)

Youhei Yamaji (The University of Tokyo)

Definition at line 90 of file CalcSpectrum.c.

References EDMainCalStruct::Bind, c_CalcExcitedStateEnd, c_CalcExcitedStateStart, c_CalcSpectrumEnd, c_CalcSpectrumStart, c_InputEigenVectorEnd, c_InputEigenVectorStart, CalcSpectrumByBiCG(), CalcSpectrumByFullDiag(), CalcSpectrumByLanczos(), CalcSpectrumByTPQ(), cFileNameTimeKeep, BindStruct::Check, childfopenALL(), D_FileNameMax, DefineList::dcOmegaMax, DefineList::dcOmegaMin, DefineList::dcOmegaOrg, BindStruct::Def, diagonalcalc(), exitMPI(), FALSE, GetExcitedState(), GetFileNameByKW(), DefineList::iCalcType, CheckList::idim_max, CheckList::idim_maxOrg, DefineList::iFlagListModified, DefineList::iFlgCalcSpec, DefineList::iFlgSpecOmegaOrg, DefineList::iNOmega, LargeList::itr, BindStruct::Large, list_1_org, list_2_1_org, list_2_2_org, MakeExcitedList(), myrank, NormMPI_dc(), DefineList::NPairExcitationOperator, DefineList::NSingleExcitationOperator, OutputSpectrum(), SetOmega(), StartTimer(), stdoutMPI, StopTimer(), TimeKeeper(), TRUE, v0, v1, v1Org, vg, and X.

Referenced by main().

                   {
    char sdt[D_FileNameMax];
    char *defname;
    unsigned long int i;
    unsigned long int i_max = 0;
    int i_stp;
    int iFlagListModified = FALSE;
    FILE *fp;
    double dnorm;

    //ToDo: Nomega should be given as a parameter
    int Nomega;
    double complex OmegaMax, OmegaMin;
    double complex *dcSpectrum;
    double complex *dcomega;
    size_t byte_size;

    //set omega
    if (SetOmega(&(X->Bind.Def)) != TRUE) {
        fprintf(stderr, "Error: Fail to set Omega.\n");
        exitMPI(-1);
    } else {
        if (X->Bind.Def.iFlgSpecOmegaOrg == FALSE) {
            X->Bind.Def.dcOmegaOrg = I*(X->Bind.Def.dcOmegaMax - X->Bind.Def.dcOmegaMin) / (double) X->Bind.Def.iNOmega;
        }
    }
    /*
     Set & malloc omega grid
    */
    Nomega = X->Bind.Def.iNOmega;
    c_malloc1(dcSpectrum, Nomega);
    c_malloc1(dcomega, Nomega);
    OmegaMax = X->Bind.Def.dcOmegaMax + X->Bind.Def.dcOmegaOrg;
    OmegaMin = X->Bind.Def.dcOmegaMin + X->Bind.Def.dcOmegaOrg;
    for (i = 0; i < Nomega; i++) {
        dcomega[i] = (OmegaMax - OmegaMin) / Nomega * i + OmegaMin;
    }
    fprintf(stdoutMPI, "\nFrequency range:\n");
    fprintf(stdoutMPI, "  Omega Max. : %15.5e %15.5e\n", creal(OmegaMax), cimag(OmegaMax));
    fprintf(stdoutMPI, "  Omega Min. : %15.5e %15.5e\n", creal(OmegaMin), cimag(OmegaMin));
    fprintf(stdoutMPI, "  Num. of Omega : %d\n", Nomega);

    if (X->Bind.Def.NSingleExcitationOperator == 0 && X->Bind.Def.NPairExcitationOperator == 0) {
        fprintf(stderr, "Error: Any excitation operators are not defined.\n");
        exitMPI(-1);
    }
    //Make New Lists
    if (MakeExcitedList(&(X->Bind), &iFlagListModified) == FALSE) {
        return FALSE;
    }
    X->Bind.Def.iFlagListModified=iFlagListModified;

    //Set Memory
    c_malloc1(v1Org, X->Bind.Check.idim_maxOrg+1);
    for(i=0; i<X->Bind.Check.idim_maxOrg+1; i++){
      v1Org[i]=0;
    }
    
    //Make excited state
    StartTimer(6100);
    if (X->Bind.Def.iFlgCalcSpec == RECALC_NOT ||
        X->Bind.Def.iFlgCalcSpec == RECALC_OUTPUT_TMComponents_VEC ||
       (X->Bind.Def.iFlgCalcSpec == RECALC_INOUT_TMComponents_VEC && X->Bind.Def.iCalcType == CG)) {
        //input eigen vector
      StartTimer(6101);
        fprintf(stdoutMPI, "  Start: An Eigenvector is inputted in CalcSpectrum.\n");
        TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_InputEigenVectorStart, "a");
        GetFileNameByKW(KWSpectrumVec, &defname);
        strcat(defname, "_rank_%d.dat");
//    sprintf(sdt, cFileNameInputEigen, X->Bind.Def.CDataFileHead, X->Bind.Def.k_exct - 1, myrank);
        sprintf(sdt, defname, myrank);
        childfopenALL(sdt, "rb", &fp);

        if (fp == NULL) {
            fprintf(stderr, "Error: A file of Inputvector does not exist.\n");
            return -1;
        }

        byte_size = fread(&i_stp, sizeof(i_stp), 1, fp);
        X->Bind.Large.itr = i_stp; //For TPQ
        byte_size = fread(&i_max, sizeof(i_max), 1, fp);
        if (i_max != X->Bind.Check.idim_maxOrg) {
            fprintf(stderr, "Error: myrank=%d, i_max=%ld\n", myrank, i_max);
            fprintf(stderr, "Error: A file of Inputvector is incorrect.\n");
            return -1;
        }
        byte_size = fread(v1Org, sizeof(complex double), i_max + 1, fp);
        fclose(fp);
        StopTimer(6101);
        if (byte_size == 0) printf("byte_size: %d \n", (int)byte_size);

        for (i = 0; i <= X->Bind.Check.idim_max; i++) {
            v0[i] = 0;
        }
        fprintf(stdoutMPI, "  End:   An Inputcector is inputted in CalcSpectrum.\n\n");
        TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_InputEigenVectorEnd, "a");
        TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_CalcExcitedStateStart, "a");
        fprintf(stdoutMPI, "  Start: Calculating an excited Eigenvector.\n");

      //mltply Operator
        StartTimer(6102);
        GetExcitedState(&(X->Bind), v0, v1Org);
        StopTimer(6102);
      
      //calculate norm
        dnorm = NormMPI_dc(X->Bind.Check.idim_max, v0);
        if (fabs(dnorm) < pow(10.0, -15)) {
            fprintf(stderr, "Warning: Norm of an excitation vector becomes 0.\n");
            fprintf(stdoutMPI, "  End:   Calculating an excited Eigenvector.\n\n");
            TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_CalcExcitedStateEnd, "a");
            fprintf(stdoutMPI, "  End:  Calculating a spectrum.\n\n");
            TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_CalcSpectrumEnd, "a");
            for (i = 0; i < Nomega; i++) {
              dcSpectrum[i]=0;
            }
            OutputSpectrum(X, Nomega, dcSpectrum, dcomega);
            return TRUE;
        }
        //normalize vector
#pragma omp parallel for default(none) private(i) shared(v1, v0) firstprivate(i_max, dnorm, X)
        for (i = 1; i <= X->Bind.Check.idim_max; i++) {
            v1[i] = v0[i] / dnorm;
        }

        fprintf(stdoutMPI, "  End:   Calculating an excited Eigenvector.\n\n");
        TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_CalcExcitedStateEnd, "a");
    }
    StopTimer(6100);
    //Reset list_1, list_2_1, list_2_2
    if (iFlagListModified == TRUE) {
      free(v1Org);
      free(list_1_org);
      free(list_2_1_org);
      free(list_2_2_org);
    }
    //calculate Diagonal term
    diagonalcalc(&(X->Bind));

  
  int iret=TRUE;
  fprintf(stdoutMPI, "  Start: Calculating a spectrum.\n\n");
  TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_CalcSpectrumStart, "a");
  StartTimer(6200);
  switch (X->Bind.Def.iCalcType) {
    case Lanczos:

      iret = CalcSpectrumByLanczos(X, v1, dnorm, Nomega, dcSpectrum, dcomega);

          if (iret != TRUE) {
            //Error Message will be added.
            return FALSE;
          }

          break;//Lanczos Spectrum

    case CG:

      iret = CalcSpectrumByBiCG(X, v0, v1, vg, Nomega, dcSpectrum, dcomega);

      if (iret != TRUE) {
        //Error Message will be added.
        return FALSE;
      }

      break;//Lanczos Spectrum

    case TPQCalc:
      fprintf(stderr, "  Error: TPQ is not supported for calculating spectrum mode.\n");
      return FALSE;//TPQ is not supprted.
#ifdef _CALCSPEC_TPQ
      iret = CalcSpectrumByTPQ(X, v1, dnorm, Nomega, dcSpectrum, dcomega);
          if (iret != TRUE) {
            //Error Message will be added.
            return FALSE;
          }
#endif

    case FullDiag:
      iret = CalcSpectrumByFullDiag(X, Nomega, dcSpectrum, dcomega);
          break;

          // case CalcSpecByShiftedKlyrov will be added
    default:
      break;
  }
  StopTimer(6200);
  
  if (iret != TRUE) {
    fprintf(stderr, "  Error: The selected calculation type is not supported for calculating spectrum mode.\n");
    return FALSE;
  }
  else {
    fprintf(stdoutMPI, "  End:  Calculating a spectrum.\n\n");
    TimeKeeper(&(X->Bind), cFileNameTimeKeep, c_CalcSpectrumEnd, "a");
    iret = OutputSpectrum(X, Nomega, dcSpectrum, dcomega);
    return TRUE;
  }

  c_free1(dcSpectrum, Nomega);
  c_free1(dcomega, Nomega);

}/*int CalcSpectrum*/

GetExcitedState()

int GetExcitedState	(	struct BindStruct *	X,
		double complex *	tmp_v0,
		double complex *	tmp_v1
	)

Parent function to calculate the excited state.

Parameters

X	[in] Struct to get number of excitation operators.
tmp_v0	[out] Result \( v_0 = H_{ex} v_1 \).
tmp_v1	[in] The original state before excitation \( v_1 \).

Return values

FALSE	Fail to calculate the excited state.
TRUE	Success to calculate the excited state.

Definition at line 303 of file CalcSpectrum.c.

References FALSE, GetPairExcitedState(), GetSingleExcitedState(), TRUE, and X.

Referenced by CalcSpectrum().

{
   if(X->Def.NSingleExcitationOperator > 0 && X->Def.NPairExcitationOperator > 0){
    fprintf(stderr, "Error: Both single and pair excitation operators exist.\n");
    return FALSE;
    }


    if(X->Def.NSingleExcitationOperator > 0){
      if(GetSingleExcitedState(X,tmp_v0, tmp_v1)!=TRUE){
        return FALSE;
      }
    }
    else if(X->Def.NPairExcitationOperator >0){
      if(GetPairExcitedState(X,tmp_v0, tmp_v1)!=TRUE){
        return FALSE;
      }
    }

  return TRUE;
}

MakeExcitedList()

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). The original lists before the excitation are given by list_xxx_org.

Parameters

X	[in, out] Struct to get and give information to make the lists for the excited state. Output: iCalcModel (From HubbardNConserved to Hubbard), {Ne, Nup, Ndown, Nsite, Total2Sz} (update for MPI)
iFlgListModifed	[out] If the list is modified due to the excitation, the value becomes TRUE(1), otherwise FALSE(0).

Return values

-1	fail to make lists.
0	sucsess to make lists.

Definition at line 421 of file CalcSpectrum.c.

References cErrLargeMem, check(), exitMPI(), FALSE, FinalizeMPI(), GetlistSize(), iErrCodeMem, list_1, list_1_org, list_1buf_org, list_2_1, list_2_1_org, list_2_2, list_2_2_org, MPIFALSE, myrank, setmem_large(), stdoutMPI, sz(), TRUE, and X.

Referenced by CalcSpectrum().

  {
    long int j;
    *iFlgListModifed = FALSE;
    //To Get Original space
    if (check(X) == MPIFALSE) {
        FinalizeMPI();
        return -1;
    }

    X->Check.idim_maxOrg = X->Check.idim_max;
    X->Check.idim_maxMPIOrg = X->Check.idim_maxMPI;

    if (X->Def.NSingleExcitationOperator > 0) {
        switch (X->Def.iCalcModel) {
            case HubbardGC:
                break;
            case HubbardNConserved:
            case KondoGC:
            case Hubbard:
            case Kondo:
                *iFlgListModifed = TRUE;
                break;
            case Spin:
            case SpinGC:
                return FALSE;
        }
    } else if (X->Def.NPairExcitationOperator > 0) {
        switch (X->Def.iCalcModel) {
            case HubbardGC:
            case SpinGC:
            case HubbardNConserved:
                break;
            case KondoGC:
            case Hubbard:
            case Kondo:
            case Spin:
                if (X->Def.PairExcitationOperator[0][1] != X->Def.PairExcitationOperator[0][3]) {
                    *iFlgListModifed = TRUE;
                }
                break;
        }
    } else {
        return FALSE;
    }

    if (*iFlgListModifed == TRUE) {
        if(GetlistSize(X)==TRUE) {
            lui_malloc1(list_1_org, X->Check.idim_max + 1);
#ifdef MPI
            lui_malloc1(list_1buf_org, X->Check.idim_maxMPI + 1);
#endif // MPI
            lui_malloc1(list_2_1_org, X->Large.SizeOflist_2_1);
            lui_malloc1(list_2_2_org, X->Large.SizeOflist_2_2);
            if(list_1_org==NULL
               || list_2_1_org==NULL
               || list_2_2_org==NULL
                    )
            {
                return -1;
            }
            for(j =0; j<X->Large.SizeOflist_2_1; j++){
                list_2_1_org[j]=0;
            }
            for(j =0; j<X->Large.SizeOflist_2_2; j++){
                list_2_2_org[j]=0;
            }

        }

        if (sz(X, list_1_org, list_2_1_org, list_2_2_org) != 0) {
            return FALSE;
        }

        if (X->Def.NSingleExcitationOperator > 0) {
            switch (X->Def.iCalcModel) {
                case HubbardGC:
                    break;
                case HubbardNConserved:
                    if (X->Def.SingleExcitationOperator[0][2] == 1) { //cis
                        X->Def.Ne = X->Def.NeMPI + 1;
                    }
                    else{
                        X->Def.Ne = X->Def.NeMPI - 1;
                    }
                    break;
                case KondoGC:
                case Hubbard:
                case Kondo:
                    if (X->Def.SingleExcitationOperator[0][2] == 1) { //cis
                        X->Def.Ne = X->Def.NeMPI + 1;
                        if (X->Def.SingleExcitationOperator[0][1] == 0) {//up
                            X->Def.Nup = X->Def.NupOrg + 1;
                            X->Def.Ndown=X->Def.NdownOrg;
                        } else {//down
                            X->Def.Nup=X->Def.NupOrg;
                            X->Def.Ndown = X->Def.NdownOrg + 1;
                        }
                    } else {//ajt
                        X->Def.Ne = X->Def.NeMPI - 1;
                        if (X->Def.SingleExcitationOperator[0][1] == 0) {//up
                            X->Def.Nup = X->Def.NupOrg - 1;
                            X->Def.Ndown=X->Def.NdownOrg;

                        } else {//down
                            X->Def.Nup=X->Def.NupOrg;
                            X->Def.Ndown = X->Def.NdownOrg - 1;
                        }
                    }
                    break;
                case Spin:
                case SpinGC:
                    return FALSE;
            }
        } else if (X->Def.NPairExcitationOperator > 0) {
            X->Def.Ne=X->Def.NeMPI;
            switch (X->Def.iCalcModel) {
                case HubbardGC:
                case SpinGC:
                case HubbardNConserved:
                    break;
                case KondoGC:
                case Hubbard:
                case Kondo:
                    if (X->Def.PairExcitationOperator[0][1] != X->Def.PairExcitationOperator[0][3]) {
                      if (X->Def.PairExcitationOperator[0][1] == 0) {//up
                        X->Def.Nup = X->Def.NupOrg + 1;
                        X->Def.Ndown = X->Def.NdownOrg - 1;
                      } else {//down
                        X->Def.Nup = X->Def.NupOrg - 1;
                        X->Def.Ndown = X->Def.NdownOrg + 1;
                      }
                    }
                    break;
              case Spin:
                if (X->Def.PairExcitationOperator[0][1] != X->Def.PairExcitationOperator[0][3]) {
                  if (X->Def.iFlgGeneralSpin == FALSE) {
                    if (X->Def.PairExcitationOperator[0][1] == 0) {//down
                      X->Def.Nup = X->Def.NupOrg - 1;
                      X->Def.Ndown = X->Def.NdownOrg + 1;
                    } else {//up
                      X->Def.Nup = X->Def.NupOrg + 1;
                      X->Def.Ndown = X->Def.NdownOrg - 1;
                    }
                  }
                  else{//for general spin
                      X->Def.Total2Sz = X->Def.Total2SzMPI+2*(X->Def.PairExcitationOperator[0][1]-X->Def.PairExcitationOperator[0][3]);
                  }
                }
                break;
            }
        } else {
            return FALSE;
        }
        //Update Infomation
        X->Def.Nsite=X->Def.NsiteMPI;

        if (check(X) == MPIFALSE) {
            FinalizeMPI();
            return FALSE;
        }
    }

    //set memory
    if (setmem_large(X) != 0) {
        fprintf(stdoutMPI, cErrLargeMem, iErrCodeMem);
        exitMPI(-1);
    }

    if (sz(X, list_1, list_2_1, list_2_2) != 0) {
        return FALSE;
    }

    if(X->Def.iCalcModel==HubbardNConserved){
        X->Def.iCalcModel=Hubbard;
    }

#ifdef _DEBUG
  if (*iFlgListModifed == TRUE) {
    for(j=1; j<=X->Check.idim_maxOrg; j++){
        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]);
    }

    for(j=1; j<=X->Check.idim_max; j++){
        fprintf(stdout, "Debug2: myrank=%d, list_1[ %ld] = %ld\n", myrank, j, list_1[j]+myrank* 64);
    }
    }
#endif

    return TRUE;
}

OutputSpectrum()

int OutputSpectrum	(	struct EDMainCalStruct *	X,
		int	Nomega,
		double complex *	dcSpectrum,
		double complex *	dcomega
	)

Output spectrum.

Parameters

X	[in] Read information of the frequency origin.
Nomega	[in] A total number of discrete frequencies.
dcSpectrum	[in] Array of spectrum.
dcomega	[in] Array of discrete frequencies.

Return values

FALSE	Fail to open the output file.
TRUE	Success to output the spectrum.

Definition at line 49 of file CalcSpectrum.c.

Referenced by CalcSpectrum().

{
  FILE *fp;
  char sdt[D_FileNameMax];
  int i;

  //output spectrum
  sprintf(sdt, cFileNameCalcDynamicalGreen, X->Bind.Def.CDataFileHead);
  if(childfopenMPI(sdt, "w", &fp)!=0){
    return FALSE;
  }

  for (i = 0; i < Nomega; i++) {
    fprintf(fp, "%.10lf %.10lf %.10lf %.10lf \n",
      creal(dcomega[i]-X->Bind.Def.dcOmegaOrg), cimag(dcomega[i]-X->Bind.Def.dcOmegaOrg),
      creal(dcSpectrum[i]), cimag(dcSpectrum[i]));
  }/*for (i = 0; i < Nomega; i++)*/

  fclose(fp);
  return TRUE;
}/*int OutputSpectrum*/

SetOmega()

int SetOmega

(

struct DefineList *

X

)

Set target frequencies.

Parameters

X	[in, out] Struct to give and get the information of target frequencies. Output: dcOmegaMax, dcOmegaMin

Return values

FALSE	Fail to set frequencies.
TRUE	Success to set frequencies.

Definition at line 337 of file CalcSpectrum.c.

References cFileNameEnergy_Lanczos, cFileNameLanczosStep, childfopenMPI(), D_FileNameMax, FALSE, fgetsMPI(), LargeValue, stdoutMPI, TRUE, and X.

Referenced by CalcSpectrum().

 {
  FILE *fp;
  char sdt[D_FileNameMax],ctmp[256];
  int istp=4;
  double E1, E2, E3, E4, Emax;
    long unsigned int iline_countMax=2;
    long unsigned int iline_count=2;


  if(X->iFlgSpecOmegaMax == TRUE && X->iFlgSpecOmegaMin == TRUE){
    return TRUE;
  }
  else{
    if (X->iCalcType == Lanczos || X->iCalcType == FullDiag) {
      sprintf(sdt, cFileNameLanczosStep, X->CDataFileHead);
      childfopenMPI(sdt, "r", &fp);
      if (fp == NULL) {
        fprintf(stdoutMPI, "Error: xx_Lanczos_Step.dat does not exist.\n");
        return FALSE;
      }
      fgetsMPI(ctmp, 256, fp); //1st line is skipped
      fgetsMPI(ctmp, 256, fp); //2nd line is skipped
      while (fgetsMPI(ctmp, 256, fp) != NULL) {
        iline_count++;
      }
      iline_countMax = iline_count;
      iline_count = 2;
      rewind(fp);
      fgetsMPI(ctmp, 256, fp); //1st line is skipped
      fgetsMPI(ctmp, 256, fp); //2nd line is skipped

      while (fgetsMPI(ctmp, 256, fp) != NULL) {
        sscanf(ctmp, "stp=%d %lf %lf %lf %lf %lf\n",
          &istp,
          &E1,
          &E2,
          &E3,
          &E4,
          &Emax);
        iline_count++;
        if (iline_count == iline_countMax) break;
      }
      fclose(fp);
      if (istp < 4) {
        fprintf(stdoutMPI, "Error: Lanczos step must be greater than 4 for using spectrum calculation.\n");
        return FALSE;
      }
    }/*if (X->iCalcType == Lanczos || X->iCalcType == FullDiag)*/
    else
    {
      sprintf(sdt, cFileNameEnergy_Lanczos, X->CDataFileHead);
      childfopenMPI(sdt, "r", &fp);
      if (fp == NULL) {
        fprintf(stdoutMPI, "Error: xx_energy.dat does not exist.\n");
        return FALSE;
      }/*if (fp == NULL)*/
      fgetsMPI(ctmp, 256, fp); //1st line is skipped
      fgetsMPI(ctmp, 256, fp); //1st line is skipped
      sscanf(ctmp, "  Energy  %lf \n", &E1);
      Emax = LargeValue;
    }
    //Read Lanczos_Step
    if(X->iFlgSpecOmegaMax == FALSE){
      X->dcOmegaMax= Emax*(double)X->Nsite;
    }
    if(X->iFlgSpecOmegaMin == FALSE){
      X->dcOmegaMin= E1;
    }
  }/*Omegamax and omegamin is not specified in modpara*/

  return TRUE;
}