From f838450ad6945fbb40462a4d05827a56cd730ee5 Mon Sep 17 00:00:00 2001 From: Carlos Mejuto Zaera Date: Mon, 9 Oct 2023 17:10:30 +0200 Subject: [PATCH] Added functionalities to specify GF frequency grid. --- include/macis/gf/gf.hpp | 24 ++++++++ src/macis/gf/gf.cxx | 123 ++++++++++++++++++++++++++++------------ 2 files changed, 111 insertions(+), 36 deletions(-) diff --git a/include/macis/gf/gf.hpp b/include/macis/gf/gf.hpp index 6ed73b2c..efdf8b6c 100644 --- a/include/macis/gf/gf.hpp +++ b/include/macis/gf/gf.hpp @@ -52,8 +52,32 @@ struct GFSettings { bool writeGF = false; bool print = false; bool saveGFmats = false; + double wmin = -8.; + double wmax = 8.; + size_t nws = 2001; + double eta = 0.1; + std::string w_scale = "lin"; + bool real_g = true; + bool beta = 1.; }; +/** + * @brief Generates the frequency grid over which to evaluate the GF, + * details of which are specified in the settings instance. We + * have implemented linear, Matsubara and logarithmic frequency + * grids. + * + * @param[in] const GFSettings& settings: Settings defining the frequency + * grid. + * + * @returns std::vector >: Vector of complex frequencies + * building the grid. + * + * @author Carlos Mejuto Zaera + * @date 09/10/2023 + */ +std::vector > GetGFFreqGrid( const GFSettings& settings ); + /** * @brief Gets fermionic sign incurred by inserting an electron of spin * up in a given orbital to a given determinant. diff --git a/src/macis/gf/gf.cxx b/src/macis/gf/gf.cxx index 8dff6304..61e4c536 100644 --- a/src/macis/gf/gf.cxx +++ b/src/macis/gf/gf.cxx @@ -9,46 +9,97 @@ #include "macis/gf/gf.hpp" namespace macis { -void write_GF(const std::vector > > &GF, - const std::vector > &ws, - const std::vector &GF_orbs, const std::vector &todelete, - const bool is_part) { - using dbl = std::numeric_limits; - size_t nfreqs = ws.size(); - int GFmat_size = GF_orbs.size() - todelete.size(); - - if(GF_orbs.size() > 1) { - std::string fname = is_part ? "LanGFMatrix_ADD.dat" : "LanGFMatrix_SUB.dat"; - std::ofstream ofile(fname); - ofile.precision(dbl::max_digits10); - for(int iii = 0; iii < nfreqs; iii++) { - ofile << std::scientific << real(ws[iii]) << " " << imag(ws[iii]) << " "; - for(int jjj = 0; jjj < GFmat_size; jjj++) { - for(int lll = 0; lll < GFmat_size; lll++) - ofile << std::scientific << real(GF[iii][jjj * GFmat_size + lll]) - << " " << imag(GF[iii][jjj * GFmat_size + lll]) << " "; + + std::vector > GetGFFreqGrid( const GFSettings& settings) + { + double wmin = settings.wmin; + double wmax = settings.wmax; + double eta = settings.eta; + double beta = settings.beta; + size_t nws = settings.nws; + bool real_g = settings.real_g; + + std::complex fact(0.,0.); + if( real_g ) + fact = std::complex(1., 0.); + else + { + fact = std::complex( 0., 1.); + eta = 0.; + } + + std::vector > ws( nws, std::complex(0.,0.) ); + + std::string scale = settings.w_scale; + if( scale == "lin" ) + { + for( int i = 0; i < nws; i++ ) + ws[i] =fact * ( wmin + (wmax - wmin) / double(nws-1.) * double(i) + + std::complex(0., eta) ); + } + else if( scale == "matsubara" ) + { + if( real_g == true ) + throw(std::runtime_error("Error in GetGFFreqGrid! Asked for 'real' Matsubara frequency grid.")); + for( int i = 0; i < nws; i++ ) + ws[i] = std::complex(0., (2. * double(i) + 1.) * M_PI / beta); + } + else if( scale == "log" ) + { + if( ( wmin < 0. && wmax > 0.) || ( wmin > 0. && wmax < 0. ) ) + throw( std::runtime_error( "Error in GetGFFreqGrid! Requested grid touches or passes by 0." ) ); + for( int i = 0; i < nws; i++ ) + { + double step = std::log(wmin) + ( std::log(wmax) - std::log(wmin) ) / double(nws-1.) * double(i); + ws[i] = fact * std::exp( step ) + std::complex( 0., eta ); } - ofile << std::endl; } + else + throw(std::runtime_error("Error in GetGFFreqGrid! Frequency scale passed is not supported. Options are 'lin', 'log' and 'matsubara'.")); + + return ws; + } - std::string fname2 = is_part ? "GFMatrix_OrbitalIndices_ADD.dat" - : "GFMatrix_OrbitalIndices_SUB.dat"; - std::ofstream ofile2(fname2); - for(int iii = 0; iii < GF_orbs.size(); iii++) { - if(std::find(todelete.begin(), todelete.end(), iii) != todelete.end()) - continue; - ofile2 << GF_orbs[iii] << std::endl; + void write_GF(const std::vector > > &GF, + const std::vector > &ws, + const std::vector &GF_orbs, const std::vector &todelete, + const bool is_part) { + using dbl = std::numeric_limits; + size_t nfreqs = ws.size(); + int GFmat_size = GF_orbs.size() - todelete.size(); + + if(GF_orbs.size() > 1) { + std::string fname = is_part ? "LanGFMatrix_ADD.dat" : "LanGFMatrix_SUB.dat"; + std::ofstream ofile(fname); + ofile.precision(dbl::max_digits10); + for(int iii = 0; iii < nfreqs; iii++) { + ofile << std::scientific << real(ws[iii]) << " " << imag(ws[iii]) << " "; + for(int jjj = 0; jjj < GFmat_size; jjj++) { + for(int lll = 0; lll < GFmat_size; lll++) + ofile << std::scientific << real(GF[iii][jjj * GFmat_size + lll]) + << " " << imag(GF[iii][jjj * GFmat_size + lll]) << " "; + } + ofile << std::endl; + } + + std::string fname2 = is_part ? "GFMatrix_OrbitalIndices_ADD.dat" + : "GFMatrix_OrbitalIndices_SUB.dat"; + std::ofstream ofile2(fname2); + for(int iii = 0; iii < GF_orbs.size(); iii++) { + if(std::find(todelete.begin(), todelete.end(), iii) != todelete.end()) + continue; + ofile2 << GF_orbs[iii] << std::endl; + } + } else { + std::string fname = is_part ? "LanGF_ADD_" : "LanGF_SUB_"; + fname += std::to_string(GF_orbs[0] + 1) + "_" + + std::to_string(GF_orbs[0] + 1) + ".dat"; + std::ofstream ofile(fname); + ofile.precision(dbl::max_digits10); + for(int iii = 0; iii < nfreqs; iii++) + ofile << std::scientific << real(ws[iii]) << " " << imag(ws[iii]) << " " + << real(GF[iii][0]) << " " << imag(GF[iii][0]) << std::endl; } - } else { - std::string fname = is_part ? "LanGF_ADD_" : "LanGF_SUB_"; - fname += std::to_string(GF_orbs[0] + 1) + "_" + - std::to_string(GF_orbs[0] + 1) + ".dat"; - std::ofstream ofile(fname); - ofile.precision(dbl::max_digits10); - for(int iii = 0; iii < nfreqs; iii++) - ofile << std::scientific << real(ws[iii]) << " " << imag(ws[iii]) << " " - << real(GF[iii][0]) << " " << imag(GF[iii][0]) << std::endl; } -} } // namespace macis