136 lines
4.3 KiB
C++
136 lines
4.3 KiB
C++
/**********************************************************
|
|
|
|
This software is part of J.-S. Caux's ABACUS library.
|
|
|
|
Copyright (c) J.-S. Caux.
|
|
|
|
-----------------------------------------------------------
|
|
|
|
File: LiebLin_RAW_File_stats.cc
|
|
|
|
Purpose: Analyzes the distribution of matrix element values in a RAW file,
|
|
to help with optimization of the scanning procedure.
|
|
|
|
***********************************************************/
|
|
|
|
#include "ABACUS.h"
|
|
|
|
using namespace std;
|
|
using namespace ABACUS;
|
|
|
|
|
|
int main(int argc, char* argv[])
|
|
{
|
|
|
|
if (argc != 9) { // provide some info
|
|
|
|
cout << endl << "Welcome to ABACUS\t(copyright J.-S. Caux)." << endl;
|
|
cout << endl << "Usage of LiebLin_RAW_File_Stats executable: " << endl;
|
|
cout << endl << "Provide the following arguments:" << endl << endl;
|
|
cout << "char whichDSF \t\t Which structure factor? Options are: "
|
|
"d for rho rho, g for psi psi{dagger}, o for psi{dagger} psi" << endl;
|
|
cout << "DP c_int \t\t Value of the interaction parameter: use positive real values only" << endl;
|
|
cout << "DP L \t\t\t Length of the system: use positive real values only" << endl;
|
|
cout << "int N \t\t\t Number of particles: use positive integer values only" << endl;
|
|
cout << "int iKmin" << endl << "int iKmax \t\t Min and max momentum integers to scan over: "
|
|
"recommended values: -2*N and 2*N" << endl;
|
|
cout << "DP kBT \t\t Temperature (positive only of course)" << endl;
|
|
cout << "int Aggregate size" << endl;
|
|
}
|
|
|
|
else { // (argc == 9), correct nr of arguments
|
|
char whichDSF = *argv[1];
|
|
DP c_int = atof(argv[2]);
|
|
DP L = atof(argv[3]);
|
|
int N = atoi(argv[4]);
|
|
int iKmin = atoi(argv[5]);
|
|
int iKmax = atoi(argv[6]);
|
|
DP kBT = atof(argv[7]);
|
|
int AgSize = atoi(argv[8]);
|
|
|
|
if (AgSize < 2) ABACUSerror("Give an aggregate size > 1 in LiebLin_RAW_File_Stats.");
|
|
|
|
stringstream RAW_stringstream; string RAW_string;
|
|
Data_File_Name (RAW_stringstream, whichDSF, c_int, L, N, iKmin, iKmax, kBT, 0.0, "");
|
|
RAW_stringstream << ".raw";
|
|
RAW_string = RAW_stringstream.str(); const char* RAW_Cstr = RAW_string.c_str();
|
|
|
|
ifstream RAW_infile;
|
|
RAW_infile.open(RAW_Cstr);
|
|
if (RAW_infile.fail()) {
|
|
cout << RAW_Cstr << endl;
|
|
ABACUSerror("Could not open RAW_infile... ");
|
|
}
|
|
|
|
stringstream STAT_stringstream; string STAT_string;
|
|
Data_File_Name (STAT_stringstream, whichDSF, c_int, L, N, iKmin, iKmax, kBT, 0.0, "");
|
|
STAT_stringstream << ".stat";
|
|
STAT_string = STAT_stringstream.str(); const char* STAT_Cstr = STAT_string.c_str();
|
|
|
|
ofstream STATfile;
|
|
STATfile.open(STAT_Cstr);
|
|
if (STATfile.fail()) {
|
|
cout << STAT_Cstr << endl; ABACUSerror("Could not open STATfile.");
|
|
}
|
|
|
|
LiebLin_Bethe_State AveragingState = Canonical_Saddle_Point_State (c_int, L, N, whichDSF == 'Z' ? 0.0 : kBT);
|
|
|
|
DP Chem_Pot = Chemical_Potential (AveragingState);
|
|
Vect<DP> sumrule_factor(iKmax - iKmin + 1);
|
|
for (int ik = 0; ik < iKmax - iKmin + 1; ++ik)
|
|
sumrule_factor[ik] = Sumrule_Factor (whichDSF, AveragingState, Chem_Pot, ik, ik);
|
|
// Normalize by total number of considered momenta
|
|
DP correction_factor = 1.0/(iKmax - iKmin + 1);
|
|
if (iKmin <= 0 && iKmax >= 0 && whichDSF == 'd') { // correct for fact that RhoRho is 0 at k = 0
|
|
sumrule_factor[0] = 0.0;
|
|
correction_factor = 1.0/(iKmax - iKmin);
|
|
}
|
|
|
|
for (int ik = 0; ik < iKmax - iKmin + 1; ++ik) sumrule_factor[ik] *= correction_factor;
|
|
|
|
DP omega;
|
|
int iK;
|
|
DP ME;
|
|
DP dev;
|
|
string label;
|
|
|
|
int nread = 0;
|
|
|
|
DP srcont = 0.0;
|
|
DP abssrcont = 0.0;
|
|
DP maxsrcont = 0.0;
|
|
DP totsrcont = 0.0;
|
|
DP accumulatedsrcont = 0.0;
|
|
int naccounted = 0;
|
|
|
|
while (RAW_infile.peek() != EOF) {
|
|
|
|
RAW_infile >> omega >> iK >> ME >> dev >> label;
|
|
nread++;
|
|
|
|
if (iK >= iKmin && iK <= iKmax) {
|
|
srcont = omega * ME * ME * sumrule_factor[iK - iKmin];
|
|
abssrcont = fabs(srcont);
|
|
maxsrcont = ABACUS::max(maxsrcont, abssrcont);
|
|
totsrcont += srcont;
|
|
accumulatedsrcont += srcont;
|
|
naccounted++;
|
|
}
|
|
|
|
if (naccounted >= AgSize) {
|
|
STATfile << nread << "\t" << maxsrcont << "\t" << totsrcont/AgSize << "\t"
|
|
<< totsrcont/(AgSize * (maxsrcont > 0.0 ? maxsrcont : 1.0)) << "\t" << accumulatedsrcont << endl;
|
|
naccounted = 0;
|
|
maxsrcont = 0.0;
|
|
totsrcont = 0.0;
|
|
}
|
|
}
|
|
|
|
RAW_infile.close();
|
|
STATfile.close();
|
|
}
|
|
|
|
|
|
return(0);
|
|
}
|