Add linear interpolation for hybridization and new parameter ntau_delta

Add files via upload

test

test2

c++/triqs_cthyb/parameters.cpp

@@ -66,6 +66,7 @@ namespace triqs_cthyb {
     h5_write(grp, "n_iw", cp.n_iw);
     h5_write(grp, "n_tau", cp.n_tau);
     h5_write(grp, "n_l", cp.n_l);
+    h5_write(grp, "n_tau_delta", cp.n_tau_delta);
     h5_write(grp, "delta_interface", cp.delta_interface);
   }
 
@@ -75,6 +76,7 @@ namespace triqs_cthyb {
     h5_read(grp, "n_iw", cp.n_iw);
     h5_read(grp, "n_tau", cp.n_tau);
     h5_read(grp, "n_l", cp.n_l);
+    h5_read(grp, "n_tau_delta", cp.n_tau_delta);
     h5::try_read(grp, "delta_interface", cp.delta_interface);
     triqs::gfs::h5_read_gf_struct(grp, "gf_struct", cp.gf_struct);
   }

c++/triqs_cthyb/parameters.hpp

@@ -49,6 +49,9 @@ namespace triqs_cthyb {
     /// Number of Legendre polynomials for gf<legendre, matrix_valued>
     int n_l = 50;
 
+    // Number of tau points for Delta_tau<imtime, matrix_valued> (-1 is the convention for n_tau_delta = n_tau)
+    int n_tau_delta = -1;
+
     /// Use Delta_tau and h_loc0 as input instead of G0_iw?
     bool delta_interface = false;
 

c++/triqs_cthyb/qmc_data.hpp

@@ -53,7 +53,28 @@ namespace triqs_cthyb {
       delta_block_adaptor &operator=(delta_block_adaptor &&) = default;
 
       det_scalar_t operator()(std::pair<time_pt, int> const &x, std::pair<time_pt, int> const &y) const {
-        det_scalar_t res = delta_block[closest_mesh_pt(double(x.first - y.first))](x.second, y.second);
+        double dtau = double(x.first - y.first);
+        auto ind = delta_block.mesh().to_data_index(closest_mesh_pt(dtau));
+        det_scalar_t res, ya, yb;
+        double xa, xb;
+        bool use_linear_interpolation = false;
+        if (use_linear_interpolation) {
+          xa = delta_block.mesh().to_value(ind);
+          ya = delta_block[ind](x.second, y.second);
+          if (dtau >= xa) {
+            xb = delta_block.mesh().to_value(ind+1);
+            yb = delta_block[ind+1](x.second, y.second);
+          }
+          else {
+            xb = xa;
+            xa = delta_block.mesh().to_value(ind-1);
+            yb = ya;
+            ya = delta_block[ind-1](x.second, y.second);
+          }
+          res = ya + (dtau - xa) * (yb - ya) / (xb - xa);
+        }
+        else
+          res = delta_block[ind](x.second, y.second);
         return (x.first >= y.first ? res : -res); // x,y first are time_pt, wrapping is automatic in the - operation, but need to
                                                   // compute the sign
       }

c++/triqs_cthyb/solver_core.cpp

@@ -61,16 +61,19 @@ namespace triqs_cthyb {
   };
 
   solver_core::solver_core(constr_parameters_t const &p)
-     : beta(p.beta), gf_struct(p.gf_struct), n_iw(p.n_iw), n_tau(p.n_tau), n_l(p.n_l), delta_interface(p.delta_interface), constr_parameters(p) {
+     : beta(p.beta), gf_struct(p.gf_struct), n_iw(p.n_iw), n_tau(p.n_tau), n_l(p.n_l), n_tau_delta(p.n_tau_delta),
+       delta_interface(p.delta_interface), constr_parameters(p) {
 
-    if (p.n_tau < 2 * p.n_iw)
+    if (n_tau_delta == -1) n_tau_delta = n_tau;
+
+    if (n_tau_delta < 2 * p.n_iw && !delta_interface)
       TRIQS_RUNTIME_ERROR
          << "Must use as least twice as many tau points as Matsubara frequencies: n_iw = " << p.n_iw
-         << " but n_tau = " << p.n_tau << ".";
+         << " but n_tau_delta = " << n_tau_delta << ".";
 
     // Allocate single particle greens functions
     if (not delta_interface) _G0_iw = block_gf<imfreq>({beta, Fermion, n_iw}, gf_struct);
-    _Delta_tau = block_gf<imtime>({beta, Fermion, n_tau}, gf_struct);
+    _Delta_tau = block_gf<imtime>({beta, Fermion, n_tau_delta}, gf_struct);
   }
 
   /// -------------------------------------------------------------------------------------------

c++/triqs_cthyb/solver_core.hpp

@@ -43,7 +43,7 @@ namespace triqs_cthyb {
     gf_struct_t gf_struct; // Block structure of the Green function
     many_body_op_t _h_loc; // The local Hamiltonian = h_int + h0
     many_body_op_t _h_loc0; //noninteracting part of the local Hamiltonian
-    int n_iw, n_tau, n_l;
+    int n_iw, n_tau, n_l, n_tau_delta;
     bool delta_interface;
 
     std::vector<matrix_t> _density_matrix; // density matrix, when used in Norm mode

python/triqs_cthyb/solver.py

@@ -34,7 +34,7 @@
 
 class Solver(SolverCore):
 
-    def __init__(self, beta, gf_struct, n_iw=1025, n_tau=10001, n_l=30, delta_interface = False):
+    def __init__(self, beta, gf_struct, n_iw=1025, n_tau=10001, n_l=30, n_tau_delta=-1, delta_interface = False):
         """
         Initialise the solver.
 
@@ -53,6 +53,8 @@ def __init__(self, beta, gf_struct, n_iw=1025, n_tau=10001, n_l=30, delta_interf
                Number of imaginary time points used for the Green's functions.
         n_l : integer, optional
              Number of legendre polynomials to use in accumulations of the Green's functions.
+        n_tau_delta : integer, optional
+             Number of imaginary time points used for the hybridization
         delta_interface: bool, optional
             Are Delta_tau and Delta_infty provided as input instead of G0_iw?
         """
@@ -61,7 +63,7 @@ def __init__(self, beta, gf_struct, n_iw=1025, n_tau=10001, n_l=30, delta_interf
 
         # Initialise the core solver
         SolverCore.__init__(self, beta=beta, gf_struct=gf_struct,
-                            n_iw=n_iw, n_tau=n_tau, n_l=n_l, delta_interface = delta_interface)
+                            n_iw=n_iw, n_tau=n_tau, n_l=n_l, n_tau_delta=n_tau_delta, delta_interface = delta_interface)
 
         mesh = MeshImFreq(beta = beta, S="Fermion", n_max = n_iw)
         self.Sigma_iw = BlockGf(mesh = mesh, gf_struct = gf_struct)
@@ -72,6 +74,7 @@ def __init__(self, beta, gf_struct, n_iw=1025, n_tau=10001, n_l=30, delta_interf
         self.gf_struct = gf_struct
         self.n_iw = n_iw
         self.n_tau = n_tau
+        self.n_tau_delta = n_tau_delta
         self.delta_interface = delta_interface
         self.G_moments = None
         self.Sigma_moments = None

python/triqs_cthyb/solver_core_desc.py

@@ -137,6 +137,8 @@
 +-----------------+-------------------------+---------+-----------------------------------------------------------------+
 | n_l             | int                     | 50      | Number of Legendre polynomials for gf<legendre, matrix_valued>  |
 +-----------------+-------------------------+---------+-----------------------------------------------------------------+
+| n_tau_delta     | int                     | -1      | Number of tau points for Delta_tau<imtime, matrix_valued>       |
++-----------------+-------------------------+---------+-----------------------------------------------------------------+
 | delta_interface | bool                    | false   | Use Delta_tau and h_loc0 as input instead of G0_iw?             |
 +-----------------+-------------------------+---------+-----------------------------------------------------------------+
 """)
@@ -634,6 +636,11 @@
              initializer = """ 50 """,
              doc = r"""Number of Legendre polynomials for gf<legendre, matrix_valued>""")
 
+c.add_member(c_name = "n_tau_delta",
+             c_type = "int",
+             initializer = """ -1 """,
+             doc = r"""Number of tau points for Delta_tau<imtime, matrix_valued>""")
+
 c.add_member(c_name = "delta_interface",
              c_type = "bool",
              initializer = """ false """,