Added functionalities to specify GF frequency grid.

wavefunction91 · Oct 9, 2023 · f838450 · f838450
1 parent b797e30
commit f838450
Show file tree

Hide file tree

Showing 2 changed files with 111 additions and 36 deletions.
diff --git 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<std::complex<double> >: Vector of complex frequencies
+ *          building the grid.
+ *
+ * @author Carlos Mejuto Zaera
+ * @date 09/10/2023
+ */
+std::vector<std::complex<double> > 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
@@ -9,46 +9,97 @@
 #include "macis/gf/gf.hpp"
 
 namespace macis {
-void write_GF(const std::vector<std::vector<std::complex<double> > > &GF,
-              const std::vector<std::complex<double> > &ws,
-              const std::vector<int> &GF_orbs, const std::vector<int> &todelete,
-              const bool is_part) {
-  using dbl = std::numeric_limits<double>;
-  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<std::complex<double> > 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<double> fact(0.,0.);
+    if( real_g )
+      fact = std::complex<double>(1., 0.);
+    else
+    {
+      fact = std::complex<double>( 0., 1.);
+      eta = 0.;
+    }
+
+    std::vector<std::complex<double> > ws( nws, std::complex<double>(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<double>(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<double>(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<double>( 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<std::vector<std::complex<double> > > &GF,
+                const std::vector<std::complex<double> > &ws,
+                const std::vector<int> &GF_orbs, const std::vector<int> &todelete,
+                const bool is_part) {
+    using dbl = std::numeric_limits<double>;
+    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