164 lines
5.8 KiB
C++
164 lines
5.8 KiB
C++
/**********************************************************
|
|
|
|
This software is part of J.-S. Caux's ABACUS library.
|
|
|
|
Copyright (c) J.-S. Caux.
|
|
|
|
-----------------------------------------------------------
|
|
|
|
File: Smoothen_RAW_into_SF_LiebLin_Scaled.cc
|
|
|
|
Purpose: from a .raw file, produces .dsf (dynamical sf) file
|
|
with gaussian width a function of momentum, and a .ssf (static sf) file.
|
|
|
|
***********************************************************/
|
|
|
|
#include "ABACUS.h"
|
|
|
|
using namespace std;
|
|
using namespace ABACUS;
|
|
|
|
namespace ABACUS {
|
|
|
|
DP Smoothen_RAW_into_SF_LiebLin_Scaled (string prefix, DP L, int N, int iKmin, int iKmax, int DiK,
|
|
DP ommin, DP ommax, int Nom, DP width, DP normalization)
|
|
{
|
|
// DiK is the (half-)window in iK which is averaged over. A single iK means DiK == 0.
|
|
|
|
// ommax is omega max for .dsf file, Nom is the number of omega slots used.
|
|
|
|
// width is measured in units of the level spacing (for each momentum independently).
|
|
|
|
// 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_SF_LiebLin_Scaled");
|
|
RecMat<DP> DSFS(Nom, iKmax - iKmin + 1);
|
|
Vect_DP SSF(0.0, iKmax - iKmin + 1);
|
|
|
|
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;
|
|
Vect_DP gwidth (iKmax - iKmin + 1);
|
|
|
|
// At fixed momentum k, we take the bandwidth as bw = e_1 (k) - e_2 (k) calculated in TG limit,
|
|
// i.e. bw = k^2 + 2\pi \rho k - (-k^2 + 2\pi\rho k) = 2 k^2 for 0 < k < 2k_F = 2\pi \rho
|
|
// and bw = 2 * (2\pi \rho)^2 for 2k_F < k.
|
|
|
|
// For iK <= N, there are iK states in bw.
|
|
// For iK > N, there are N states.
|
|
|
|
for (iK = iKmin; iK <= iKmax; ++iK)
|
|
// Make sure the width does not become lower than the omegaout raster:
|
|
gwidth[iK - iKmin] = ABACUS::max(2.0 * (ommax - ommin)/Nom, width * 2.0
|
|
* ABACUS::min( pow(twoPI * ABACUS::max(abs(iK),1)/L, 2.0),
|
|
pow(twoPI * N/L, 2.0))/ABACUS::min(N, ABACUS::max(abs(iK), 1)));
|
|
|
|
|
|
Vect_DP big_gwidth_used (iKmax - iKmin + 1);
|
|
for (iK = iKmin; iK <= iKmax; ++iK)
|
|
big_gwidth_used[iK - iKmin] = 10.0 * gwidth[iK - iKmin]; // neglect terms having gaussian < exp(-50)
|
|
Vect_DP oneovertwowidthsq (iKmax - iKmin + 1);
|
|
for (iK = iKmin; iK <= iKmax; ++iK)
|
|
oneovertwowidthsq[iK - iKmin] = 1.0/(2.0 * gwidth[iK - iKmin] * gwidth[iK - iKmin]);
|
|
|
|
// Reset proper normalization:
|
|
Vect_DP normalization_used (iKmax - iKmin + 1);
|
|
for (iK = iKmin; iK <= iKmax; ++iK)
|
|
normalization_used[iK - iKmin] = normalization * 1.0/(sqrt(twoPI) * gwidth[iK - iKmin]); // Gaussian factor
|
|
|
|
DP FFsq = 0.0;
|
|
while (RAW_infile.peek() != EOF) {
|
|
RAW_infile >> omega >> iK >> FF >> dev >> label;
|
|
if (iK >= iKmin && iK <= iKmax && fabs(omega) > 1.0e-8) { // remove connected part of DSF
|
|
FFsq = FF * FF;
|
|
SSF[iK - iKmin] += FFsq;
|
|
for (int iomega = 0; iomega < Nom; ++iomega)
|
|
if (big_gwidth_used[iK - iKmin] > (d_omega = fabs(omegaout[iomega] - omega)))
|
|
for (int deltaiK = -DiK; deltaiK <= DiK; ++deltaiK)
|
|
if (iK + deltaiK >= iKmin && iK + deltaiK <= iKmax)
|
|
DSFS[iomega][iK + deltaiK - iKmin] += Kweight[abs(deltaiK)] * FFsq * normalization_used[iK + deltaiK - iKmin]
|
|
* exp(-d_omega*d_omega * oneovertwowidthsq[iK + deltaiK - iKmin]);
|
|
}
|
|
}
|
|
RAW_infile.close();
|
|
|
|
// Reset proper normalization:
|
|
for (int iK = 0; iK < iKmax - iKmin + 1; ++iK) {
|
|
SSF[iK] *= normalization/twoPI;
|
|
}
|
|
|
|
// Output to .dsfs and .ssf files
|
|
stringstream DSFS_stringstream; string DSFS_string;
|
|
DSFS_stringstream << prefix;
|
|
if (DiK > 0) DSFS_stringstream << "_DiK_" << DiK;
|
|
DSFS_stringstream << "_ommin_"<< ommin << "_ommax_" << ommax << "_Nom_" << Nom << "_w_" << width << ".dsfs";
|
|
DSFS_string = DSFS_stringstream.str(); const char* DSFS_Cstr = DSFS_string.c_str();
|
|
|
|
ofstream DSFS_outfile;
|
|
DSFS_outfile.open(DSFS_Cstr);
|
|
DSFS_outfile.precision(12);
|
|
|
|
for (int iomega = 0; iomega < Nom; ++iomega) {
|
|
if (iomega > 0) DSFS_outfile << endl;
|
|
for (int iK = 0; iK < iKmax - iKmin + 1; ++iK)
|
|
DSFS_outfile << DSFS[iomega][iK] << "\t";
|
|
}
|
|
DSFS_outfile.close();
|
|
|
|
|
|
stringstream SSF_stringstream; string SSF_string;
|
|
SSF_stringstream << prefix;
|
|
SSF_stringstream << ".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) {
|
|
SSF_outfile << iK + iKmin << "\t" << SSF[iK] << endl;
|
|
}
|
|
SSF_outfile.close();
|
|
|
|
|
|
// Check sums:
|
|
DP sumdsf = 0.0;
|
|
for (int iK = 0; iK < iKmax - iKmin + 1; ++iK) {
|
|
for (int iomega = 0; iomega < Nom; ++iomega)
|
|
sumdsf += DSFS[iomega][iK];
|
|
}
|
|
sumdsf /= (iKmax - iKmin + 1) * Nom;
|
|
|
|
return(sumdsf);
|
|
|
|
}
|
|
|
|
} // namespace ABACUS
|