399 lines
14 KiB
C++
399 lines
14 KiB
C++
/**********************************************************
|
|
|
|
This software is part of J.-S. Caux's ABACUS library.
|
|
|
|
Copyright (c) J.-S. Caux.
|
|
|
|
-----------------------------------------------------------
|
|
|
|
File: Smoothen_RAW_into_SF.cc
|
|
|
|
Purpose: from a .raw file, produces .dsf (dynamical sf) and .ssf (static sf) files.
|
|
|
|
***********************************************************/
|
|
|
|
#include "ABACUS.h"
|
|
|
|
using namespace std;
|
|
using namespace ABACUS;
|
|
|
|
namespace ABACUS {
|
|
|
|
DP Smoothen_RAW_into_SF (string prefix, int iKmin, int iKmax, int DiK, DP devmax,
|
|
DP ommin, DP ommax, int Nom, DP gwidth, DP normalization, DP denom_sum_K)
|
|
{
|
|
// ommax is omega max for .dsf file, Nom is the number of omega slots used.
|
|
|
|
// gwidth is the width of the smoothing Gaussian, defined as
|
|
// exp(-omega^2/(2 * gwidth^2))
|
|
|
|
// DiK is the (half-)window in iK which is averaged over. Averaging over a single iK means DiK == 0.
|
|
|
|
// Open the original raw file:
|
|
stringstream RAW_stringstream; string RAW_string;
|
|
RAW_stringstream << prefix << ".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... ");
|
|
}
|
|
|
|
if (iKmax - iKmin + 1 < 0) ABACUSerror("Improper iKmin, iKmax in Smoothen_RAW_into_SF");
|
|
RecMat<DP> DSF(Nom, iKmax - iKmin + 1);
|
|
Vect_DP SSF(0.0, iKmax - iKmin + 1);
|
|
Vect_DP ASF(0.0, Nom);
|
|
|
|
DP omega;
|
|
int iK;
|
|
DP FF;
|
|
DP dev;
|
|
string label;
|
|
|
|
// Momenta: average over 2*DiK + 1 entries. Weigh them linearly decreasing away from central one.
|
|
// Setting central one to value 1 + DiK,
|
|
// total weight is 1 + DiK + 2* \sum_1^DiK n = 1 + DiK + DiK (DiK + 1) = (DiK + 1)^2.
|
|
// Weight given is thus abs(DiK + 1 - (iK - iK'))/(DiK + 1)^2 for abs(iK - iK') <= DiK.
|
|
Vect_DP Kweight(DiK + 1);
|
|
for (int i = 0; i < DiK + 1; ++i) Kweight[i] = (DiK + 1.0 - i)/pow(DiK + 1.0, 2);
|
|
|
|
Vect_DP omegaout (Nom);
|
|
for (int i = 0; i < Nom; ++i) omegaout[i] = ommin + (0.5 + i) * (ommax - ommin)/Nom;
|
|
|
|
DP d_omega;
|
|
DP big_gwidth_used = 10.0 * gwidth; // neglect terms having gaussian < exp(-50)
|
|
DP oneovertwowidthsq = 1.0/(2.0 * gwidth * gwidth);
|
|
DP SFfactor = 1.0;
|
|
|
|
while (RAW_infile.peek() != EOF) {
|
|
RAW_infile >> omega >> iK >> FF >> dev >> label;
|
|
if (iK >= iKmin && iK <= iKmax && dev < devmax && fabs(omega) > 1.0e-8) { // remove connected part of DSF
|
|
for (int deltaiK = -DiK; deltaiK <= DiK; ++deltaiK)
|
|
if (iK + deltaiK >= iKmin && iK + deltaiK <= iKmax)
|
|
SSF[iK + deltaiK - iKmin] += Kweight[abs(deltaiK)] * FF * FF;
|
|
for (int iomega = 0; iomega < Nom; ++iomega)
|
|
if (big_gwidth_used > (d_omega = fabs(omegaout[iomega] - omega))) {
|
|
SFfactor = FF * FF * exp(-d_omega*d_omega * oneovertwowidthsq);
|
|
ASF[iomega] += SFfactor;
|
|
if (fabs(omega) > 1.0e-12) // exclude the delta function contribution coming from diagonal term, if present
|
|
for (int deltaiK = -DiK; deltaiK <= DiK; ++deltaiK)
|
|
if (iK + deltaiK >= iKmin && iK + deltaiK <= iKmax)
|
|
DSF[iomega][iK + deltaiK - iKmin] += Kweight[abs(deltaiK)] * SFfactor;
|
|
}
|
|
}
|
|
}
|
|
RAW_infile.close();
|
|
|
|
// Reset proper normalization:
|
|
DP normalization_used = normalization * 1.0/(sqrt(twoPI) * gwidth); // Gaussian factor
|
|
|
|
for (int iK = 0; iK < iKmax - iKmin + 1; ++iK) {
|
|
SSF[iK] *= normalization/twoPI; // twoPI from integral over omega: \int d\omega/2\pi
|
|
for (int iomega = 0; iomega < Nom; ++iomega)
|
|
DSF[iomega][iK] *= normalization_used;
|
|
}
|
|
|
|
DP ASFnormalization = normalization_used/denom_sum_K;
|
|
for (int iomega = 0; iomega < Nom; ++iomega) ASF[iomega] *= ASFnormalization;
|
|
|
|
|
|
// Output to .dsf, .ssf and .asf files
|
|
|
|
stringstream DSF_stringstream; string DSF_string;
|
|
DSF_stringstream << prefix;
|
|
if (DiK > 0) DSF_stringstream << "_DiK_" << DiK;
|
|
DSF_stringstream << "_devmax_" << devmax << "_ommin_"<< ommin << "_ommax_" << ommax << "_Nom_" << Nom << "_w_" << gwidth << ".dsf";
|
|
DSF_string = DSF_stringstream.str(); const char* DSF_Cstr = DSF_string.c_str();
|
|
|
|
ofstream DSF_outfile;
|
|
DSF_outfile.open(DSF_Cstr);
|
|
DSF_outfile.precision(12);
|
|
|
|
for (int iomega = 0; iomega < Nom; ++iomega) {
|
|
if (iomega > 0) DSF_outfile << endl;
|
|
for (int iK = 0; iK < iKmax - iKmin + 1; ++iK)
|
|
DSF_outfile << DSF[iomega][iK] << "\t";
|
|
}
|
|
DSF_outfile.close();
|
|
|
|
stringstream SSF_stringstream; string SSF_string;
|
|
SSF_stringstream << prefix;
|
|
if (DiK > 0) SSF_stringstream << "_DiK_" << DiK;
|
|
SSF_stringstream << "_devmax_" << devmax << ".ssf";
|
|
SSF_string = SSF_stringstream.str(); const char* SSF_Cstr = SSF_string.c_str();
|
|
|
|
ofstream SSF_outfile;
|
|
SSF_outfile.open(SSF_Cstr);
|
|
SSF_outfile.precision(12);
|
|
|
|
for (int iK = 0; iK < iKmax - iKmin + 1; ++iK) {
|
|
if (iK > 0) SSF_outfile << endl;
|
|
SSF_outfile << iK + iKmin << "\t" << SSF[iK];
|
|
}
|
|
SSF_outfile.close();
|
|
|
|
|
|
stringstream ASF_stringstream; string ASF_string;
|
|
ASF_stringstream << prefix << "_devmax_" << devmax;
|
|
ASF_stringstream << "_ommin_"<< ommin << "_ommax_" << ommax << "_Nom_" << Nom << "_w_" << gwidth << ".asf";
|
|
ASF_string = ASF_stringstream.str(); const char* ASF_Cstr = ASF_string.c_str();
|
|
|
|
ofstream ASF_outfile;
|
|
ASF_outfile.open(ASF_Cstr);
|
|
ASF_outfile.precision(12);
|
|
|
|
for (int iomega = 0; iomega < Nom; ++iomega) {
|
|
if (iomega > 0) ASF_outfile << endl;
|
|
ASF_outfile << omegaout[iomega] << "\t" << ASF[iomega];
|
|
}
|
|
ASF_outfile.close();
|
|
|
|
|
|
|
|
// Check sums:
|
|
DP sumdsf = 0.0;
|
|
DP sumssf = 0.0;
|
|
for (int iK = 0; iK < iKmax - iKmin + 1; ++iK) {
|
|
sumssf += SSF[iK];
|
|
for (int iomega = 0; iomega < Nom; ++iomega)
|
|
sumdsf += DSF[iomega][iK];
|
|
}
|
|
sumssf /= (iKmax - iKmin + 1);
|
|
sumdsf /= (iKmax - iKmin + 1) * Nom;
|
|
|
|
return(sumdsf);
|
|
}
|
|
|
|
// This is the same function as above, but now using data for a diagonal ensemble of states
|
|
DP Smoothen_RAW_into_SF (string prefix, Vect<string> rawfilename, Vect<DP> weight, int iKmin, int iKmax, int DiK, DP devmax,
|
|
DP ommin, DP ommax, int Nom, DP gwidth, DP normalization, DP denom_sum_K)
|
|
{
|
|
// ommax is omega max for .dsf file, Nom is the number of omega slots used.
|
|
|
|
// gwidth is the width of the smoothing Gaussian, defined as
|
|
// exp(-omega^2/(2 * gwidth^2))
|
|
|
|
// DiK is the (half-)window in iK which is averaged over. Averaging over a single iK means DiK == 0.
|
|
|
|
if (iKmax - iKmin + 1 < 0) ABACUSerror("Improper iKmin, iKmax in Smoothen_RAW_into_SF");
|
|
RecMat<DP> DSF(Nom, iKmax - iKmin + 1);
|
|
Vect_DP SSF(0.0, iKmax - iKmin + 1);
|
|
Vect_DP ASF(0.0, Nom);
|
|
|
|
DP omega;
|
|
int iK;
|
|
DP FF;
|
|
DP dev;
|
|
string label;
|
|
|
|
// Momenta: average over 2*DiK + 1 entries. Weigh them linearly decreasing away from central one.
|
|
// Setting central one to value 1 + DiK,
|
|
// total weight is 1 + DiK + 2* \sum_1^DiK n = 1 + DiK + DiK (DiK + 1) = (DiK + 1)^2.
|
|
// Weight given is thus abs(DiK + 1 - (iK - iK'))/(DiK + 1)^2 for abs(iK - iK') <= DiK.
|
|
Vect_DP Kweight(DiK + 1);
|
|
for (int i = 0; i < DiK + 1; ++i) Kweight[i] = (DiK + 1.0 - i)/pow(DiK + 1.0, 2);
|
|
|
|
Vect_DP omegaout (Nom);
|
|
for (int i = 0; i < Nom; ++i) omegaout[i] = ommin + (0.5 + i) * (ommax - ommin)/Nom;
|
|
|
|
DP d_omega;
|
|
DP big_gwidth_used = 10.0 * gwidth; // neglect terms having gaussian < exp(-50)
|
|
DP oneovertwowidthsq = 1.0/(2.0 * gwidth * gwidth);
|
|
DP SFfactor = 1.0;
|
|
|
|
for (int ns = 0; ns < weight.size(); ++ns) {
|
|
|
|
// Open the original raw file:
|
|
const char* RAW_Cstr = rawfilename[ns].c_str();
|
|
|
|
ifstream RAW_infile;
|
|
RAW_infile.open(RAW_Cstr);
|
|
if (RAW_infile.fail()) {
|
|
cout << RAW_Cstr << endl;
|
|
ABACUSerror("Could not open RAW_infile... ");
|
|
}
|
|
|
|
while (RAW_infile.peek() != EOF) {
|
|
RAW_infile >> omega >> iK >> FF >> dev >> label;
|
|
if (iK >= iKmin && iK <= iKmax && dev < devmax && fabs(omega) > 1.0e-8) { // remove connected part of DSF)
|
|
SSF[iK - iKmin] += weight[ns] * FF * FF;
|
|
for (int iomega = 0; iomega < Nom; ++iomega)
|
|
if (big_gwidth_used > (d_omega = fabs(omegaout[iomega] - omega))) {
|
|
SFfactor = weight[ns] * FF * FF * exp(-d_omega*d_omega * oneovertwowidthsq);
|
|
ASF[iomega] += SFfactor;
|
|
if (fabs(omega) > 1.0e-12) // exclude the delta function contribution coming from diagonal term, if present
|
|
for (int deltaiK = -DiK; deltaiK <= DiK; ++deltaiK)
|
|
if (iK + deltaiK >= iKmin && iK + deltaiK <= iKmax)
|
|
DSF[iomega][iK + deltaiK - iKmin] += Kweight[abs(deltaiK)] * SFfactor;
|
|
}
|
|
}
|
|
}
|
|
RAW_infile.close();
|
|
} // for ns
|
|
|
|
// Reset proper normalization:
|
|
DP normalization_used = normalization * 1.0/(sqrt(twoPI) * gwidth); // Gaussian factor
|
|
|
|
for (int iK = 0; iK < iKmax - iKmin + 1; ++iK) {
|
|
SSF[iK] *= normalization/twoPI; // twoPI from integral over omega: \int d\omega/2\pi
|
|
for (int iomega = 0; iomega < Nom; ++iomega)
|
|
DSF[iomega][iK] *= normalization_used;
|
|
}
|
|
|
|
DP ASFnormalization = normalization_used/denom_sum_K;
|
|
for (int iomega = 0; iomega < Nom; ++iomega) ASF[iomega] *= ASFnormalization;
|
|
|
|
|
|
// Output to .dsf, .ssf and .asf files
|
|
|
|
stringstream DSF_stringstream; string DSF_string;
|
|
DSF_stringstream << prefix;
|
|
DSF_stringstream << "_ns_" << weight.size();
|
|
if (DiK > 0) DSF_stringstream << "_DiK_" << DiK;
|
|
DSF_stringstream << "_devmax_" << devmax << "_ommin_"<< ommin << "_ommax_" << ommax << "_Nom_" << Nom << "_w_" << gwidth << ".dsf";
|
|
DSF_string = DSF_stringstream.str(); const char* DSF_Cstr = DSF_string.c_str();
|
|
|
|
ofstream DSF_outfile;
|
|
DSF_outfile.open(DSF_Cstr);
|
|
DSF_outfile.precision(12);
|
|
|
|
for (int iomega = 0; iomega < Nom; ++iomega) {
|
|
if (iomega > 0) DSF_outfile << endl;
|
|
for (int iK = 0; iK < iKmax - iKmin + 1; ++iK)
|
|
DSF_outfile << DSF[iomega][iK] << "\t";
|
|
}
|
|
DSF_outfile.close();
|
|
|
|
stringstream SSF_stringstream; string SSF_string;
|
|
SSF_stringstream << prefix;
|
|
SSF_stringstream << "_ns_" << weight.size();
|
|
SSF_stringstream << "_devmax_" << devmax << ".ssf";
|
|
SSF_string = SSF_stringstream.str(); const char* SSF_Cstr = SSF_string.c_str();
|
|
|
|
ofstream SSF_outfile;
|
|
SSF_outfile.open(SSF_Cstr);
|
|
SSF_outfile.precision(12);
|
|
|
|
for (int iK = 0; iK < iKmax - iKmin + 1; ++iK) {
|
|
if (iK > 0) SSF_outfile << endl;
|
|
SSF_outfile << iK + iKmin << "\t" << SSF[iK];
|
|
}
|
|
SSF_outfile.close();
|
|
|
|
|
|
stringstream ASF_stringstream; string ASF_string;
|
|
ASF_stringstream << prefix;
|
|
ASF_stringstream << "_ns_" << weight.size();
|
|
ASF_stringstream << "_devmax_" << devmax << "_ommin_"<< ommin << "_ommax_" << ommax << "_Nom_" << Nom << "_w_" << gwidth << ".asf";
|
|
ASF_string = ASF_stringstream.str(); const char* ASF_Cstr = ASF_string.c_str();
|
|
|
|
ofstream ASF_outfile;
|
|
ASF_outfile.open(ASF_Cstr);
|
|
ASF_outfile.precision(12);
|
|
|
|
for (int iomega = 0; iomega < Nom; ++iomega) {
|
|
if (iomega > 0) ASF_outfile << endl;
|
|
ASF_outfile << omegaout[iomega] << "\t" << ASF[iomega];
|
|
}
|
|
ASF_outfile.close();
|
|
|
|
|
|
// Check sums:
|
|
DP sumdsf = 0.0;
|
|
DP sumssf = 0.0;
|
|
for (int iK = 0; iK < iKmax - iKmin + 1; ++iK) {
|
|
sumssf += SSF[iK];
|
|
for (int iomega = 0; iomega < Nom; ++iomega)
|
|
sumdsf += DSF[iomega][iK];
|
|
}
|
|
sumssf /= (iKmax - iKmin + 1);
|
|
sumdsf /= (iKmax - iKmin + 1) * Nom;
|
|
|
|
return(sumdsf);
|
|
}
|
|
|
|
|
|
|
|
|
|
DP Smoothen_RAW_into_ASF (string prefix, int iKmin, int iKmax, DP devmax, DP ommin, DP ommax, int Nom, DP gwidth,
|
|
DP normalization, DP denom_sum_K)
|
|
{
|
|
// ommax is omega max for .asf file, Nom is the number of omega slots used.
|
|
|
|
// gwidth is the width of the smoothing Gaussian, defined as
|
|
// exp(-omega^2/(2 * gwidth^2))
|
|
|
|
// Open the original raw file:
|
|
stringstream RAW_stringstream; string RAW_string;
|
|
RAW_stringstream << prefix << ".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... ");
|
|
}
|
|
|
|
if (iKmax - iKmin + 1 < 0) ABACUSerror("Improper iKmin, iKmax in Smoothen_RAW_into_ASF");
|
|
Vect_DP ASF(Nom);
|
|
|
|
DP omega;
|
|
int iK;
|
|
DP FF;
|
|
DP dev;
|
|
string label;
|
|
|
|
Vect_DP omegaout (Nom);
|
|
for (int i = 0; i < Nom; ++i) omegaout[i] = ommin + (0.5 + i) * (ommax - ommin)/Nom;
|
|
|
|
DP d_omega;
|
|
DP big_gwidth_used = 10.0 * gwidth; // neglect terms having gaussian < exp(-50)
|
|
DP oneovertwowidthsq = 1.0/(2.0 * gwidth * gwidth);
|
|
|
|
while (RAW_infile.peek() != EOF) {
|
|
RAW_infile >> omega >> iK >> FF >> dev >> label;
|
|
if (iK >= iKmin && iK <= iKmax && dev < devmax && fabs(omega) > 1.0e-8) { // remove connected part of DSF)
|
|
for (int iomega = 0; iomega < Nom; ++iomega)
|
|
if (big_gwidth_used > (d_omega = fabs(omegaout[iomega] - omega)))
|
|
ASF[iomega] += FF * FF * exp(-d_omega*d_omega * oneovertwowidthsq);
|
|
}
|
|
}
|
|
RAW_infile.close();
|
|
|
|
// Reset proper normalization:
|
|
DP normalization_used = normalization * 1.0/(sqrt(twoPI) * gwidth); // Gaussian factor
|
|
normalization_used /= denom_sum_K;
|
|
for (int iomega = 0; iomega < Nom; ++iomega)
|
|
ASF[iomega] *= normalization_used;
|
|
|
|
// Output to .asf file
|
|
stringstream ASF_stringstream; string ASF_string;
|
|
ASF_stringstream << prefix;
|
|
//if (iKmax != iKmin) DSF_stringstream << "_iKmin_" << iKmin << "_iKmax_" << iKmax;
|
|
ASF_stringstream << "_devmax_" << devmax << "_ommin_"<< ommin << "_ommax_" << ommax << "_Nom_" << Nom << "_w_" << gwidth << ".asf";
|
|
ASF_string = ASF_stringstream.str(); const char* ASF_Cstr = ASF_string.c_str();
|
|
|
|
ofstream ASF_outfile;
|
|
ASF_outfile.open(ASF_Cstr);
|
|
ASF_outfile.precision(12);
|
|
|
|
for (int iomega = 0; iomega < Nom; ++iomega) {
|
|
if (iomega > 0) ASF_outfile << endl;
|
|
ASF_outfile << ASF[iomega] << "\t";
|
|
}
|
|
ASF_outfile.close();
|
|
|
|
DP sumasf = 0.0;
|
|
for (int iomega = 0; iomega < Nom; ++iomega)
|
|
sumasf += ASF[iomega];
|
|
|
|
sumasf /= Nom;
|
|
|
|
return(sumasf);
|
|
}
|
|
|
|
} // namespace ABACUS
|