update benchmarks, add coments

Author: nkavokine

benchmark/dimer/dimer.py

@@ -0,0 +1,109 @@
+# Two impurity sites coupled to two bath sites.
+# Data in dimer.ref.h5 is obtained by running this script on 800 cores. 
+
+from triqs.gf import *
+from triqs.gf.tools import *
+from triqs.gf.gf_factories import make_gf_from_fourier
+from triqs.operators.util import U_matrix_kanamori, h_int_density
+import h5
+import triqs.utility.mpi as mpi
+from triqs.utility.h5diff import h5diff
+from triqs.operators import n
+
+from triqs_ctseg import Solver
+
+import numpy as np
+from numpy import linalg 
+
+# Numerical values
+beta = 10.                      # Inverse temperature
+mu = 0.0                        # Chemical potential
+eps = np.array([0.0, 0.1])         # Impurity site energies
+eps_bath = np.array([0.27, -0.4])  # Bath site energies
+U = 1.                          # Density-density interaction
+J = 0.2                         # Hund's coupling
+V = 1                           # Bath-impurity coupling
+block_names = ['up', 'dn']
+n_orb = len(eps)
+n_orb_bath = len(eps_bath)
+n_tau = 10001
+gf_struct = [(s, n_orb) for s in block_names] # Green's function structure
+
+########################################################################
+# --------------- Construct Delta(tau) and h_loc0 ---------------
+
+# Non-interacting impurity Hamiltonian in matrix representation
+h_0_mat = np.diag(eps - mu)
+
+# Bath Hamiltonian in matrix representation
+h_bath_mat = np.diag(eps_bath)
+
+# Coupling matrix
+V_mat = np.matrix([[V, V],
+                [V, V]])
+
+# Total non-interacting Hamiltonian 
+h_tot_mat = np.block([[h_0_mat, V_mat],
+                   [V_mat.H, h_bath_mat]])
+
+# Non-interacting impurity Green's function
+n_iw = int(10 * beta)
+iw_mesh = MeshImFreq(beta, 'Fermion', n_iw)
+G0_iw = BlockGf(mesh=iw_mesh, gf_struct=gf_struct)
+for bl, iw in product(block_names, iw_mesh):
+    G0_iw[bl][iw] = linalg.inv(iw.value * np.eye(2*n_orb) - h_tot_mat)[:n_orb, :n_orb]
+
+# Get Delta(iw) and h_loc0 from G0(iw)
+def get_h0_Delta(G0_iw):
+    h0_lst, Delta_iw = [], G0_iw.copy()
+    for bl in G0_iw.indices:
+        Delta_iw[bl] << iOmega_n - inverse(G0_iw[bl])
+        tail, err = fit_hermitian_tail(Delta_iw[bl])
+        Delta_iw[bl] << Delta_iw[bl] - tail[0]
+        h0_lst.append(tail[0].real)
+    return h0_lst, Delta_iw
+
+h0_lst, Delta_iw = get_h0_Delta(G0_iw)
+h_loc0 = sum(h0_lst[i][0,0] * n(bl, 0) + h0_lst[i][1,1] * n(bl, 1) for i, bl in enumerate(block_names))
+
+# Fourier-transform to get Delta(tau)
+tau_mesh = MeshImTime(beta, 'Fermion', n_tau)
+Delta_tau = BlockGf(mesh = tau_mesh, gf_struct = gf_struct)
+for name, g0 in Delta_iw:
+    known_moments = make_zero_tail(Delta_iw[name], 1)
+    tail, err = fit_hermitian_tail(Delta_iw[name], known_moments)
+    Delta_tau[name] << make_gf_from_fourier(Delta_iw[name], tau_mesh, tail).real
+###########################################################
+
+# --------- Construct the CTSEG solver ----------
+constr_params = {
+    'beta': beta,
+    'gf_struct': gf_struct,
+    'n_tau': n_tau,
+}
+S = Solver(**constr_params)
+
+# Input Delta(tau)
+S.Delta_tau << Delta_tau
+
+# Impurity interaction Hamiltonian
+Umat, Upmat = U_matrix_kanamori(n_orb, U_int=U, J_hund=J)
+h_int = h_int_density(block_names, n_orb, Umat, Upmat, off_diag=True)
+
+# --------- Solve parameters ----------
+solve_params = {
+    'h_int': h_int,
+    'h_loc0': h_loc0,
+    'n_warmup_cycles': 5000,
+    'n_cycles': 200000000,
+    'length_cycle': 100,
+}
+
+# ---------- Solve ----------
+S.solve(**solve_params)
+
+# --------- Save output ---------
+if mpi.is_master_node():
+    with h5.HDFArchive("dimer.out.h5", 'w') as A:
+        A['G_tau'] = S.results.G_tau
+

benchmark/move_double/ctint.ref.h5

@@ -4,9 +4,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### Compare `spin_spin.py` result with `ctint`\n",
+    "### Comparison between `CTSEG` and `CTINT` reference. \n",
     "\n",
-    "We use black line to denote the ctint result. "
+    "Single-orbital impurity with retarded spin-spin interactions. "
    ]
   },
   {
@@ -29,11 +29,7 @@
     "mpl.rcParams['figure.dpi'] = 150\n",
     "plt.rcParams['mathtext.fontset'] = 'cm'\n",
     "plt.rcParams['mathtext.rm'] = 'serif'\n",
-    "plt.rc('font', size=12, family='Times New Roman')\n",
-    "\n",
-    "def normalize_color(color):\n",
-    "    brightness = color[0]**2 + color[1]**2 + color[2]**2\n",
-    "    return (color[0]/brightness, color[1]/brightness, color[2]/brightness, color[3])\n"
+    "plt.rc('font', size=12, family='Times New Roman')"
    ]
   },
   {
@@ -64,9 +60,9 @@
     "with HDFArchive(\"spin_spin.ref.h5\", \"r\") as A:\n",
     "    G_tau_seg = A[\"G_tau\"]\n",
     "\n",
-    "oplot(G_tau_int[\"up\"].real, 'k', linewidth = 2, label=\"ctint\")\n",
-    "oplot(G_tau_seg[\"up\"].real, color=\"aqua\", linewidth = 0.6, label=\"ctseg\")\n",
-    "plt.title(\"Spin Spin\")\n",
+    "oplot(G_tau_int[\"up\"].real, 'k', linewidth = 2, label=\"CTINT\")\n",
+    "oplot(G_tau_seg[\"up\"].real, color=\"aqua\", linewidth = 0.6, label=\"CTSEG\")\n",
+    "plt.title(\"Impurity with spin-spin interactions\")\n",
     "plt.ylabel(r\"$G(\\tau)$\")\n",
     "plt.legend()\n",
     "plt.show()"

benchmark/move_double/multiorb.py

@@ -1,5 +1,5 @@
 # Single orbital with dynamical spin-spin interactions. 
-# Can be compared with reference obtained with CTINT (ctint.ref.h5)
+# Data in spin_spin.ref.h5 is obtained by running this script on 800 cores. 
 from triqs.gf import *
 import triqs.utility.mpi as mpi
 from triqs.gf.descriptors import Function
@@ -57,20 +57,17 @@
     "n_cycles": 1000000,
     "measure_F_tau": True,
     "measure_nn_tau": True,
-    "measure_nn_static": True,
-    "move_double_insert_segment": True,
-    "move_double_remove_segment": True,
+    "measure_nn_static": True
     }
 
 # Solve
 S.solve(**solve_params)
 
-# Save and compare to reference
+# Save data
 if mpi.is_master_node():
     with h5.HDFArchive("spin_spin.out.h5", 'w') as A:
         A['G_tau'] = S.results.G_tau
         A['F_tau'] = S.results.F_tau
         A['nn_tau'] = S.results.nn_tau
         A['nn'] = S.results.nn_static
         A['densities'] = S.results.densities
-

benchmark/move_double/multiorb.ref.h5

@@ -25,7 +25,7 @@ namespace triqs_ctseg::moves {
 
     LOG("\n =================== ATTEMPT DOUBLE INSERT ================ \n");
 
-    // ------------ Choice of segment --------------
+    // ------------ Choice of segments --------------
     // Select insertion colors
     int color_0 = rng(config.n_color());
     int color_1 = rng(config.n_color() - 1);
@@ -90,12 +90,12 @@ namespace triqs_ctseg::moves {
         ln_trace_ratio += -real(wdata.K(double(prop_seg[i].length()))(color, color)); // Correct double counting
     } // color
 
-    // Counting the overlap between the inserting segments
-    // Make the prop_seg[1] as a seglist to use K_overlap()
-    std::vector<segment_t> seglist_temp = {prop_seg[1]};
+    // Overlap between the inserted segments
     ln_trace_ratio += -wdata.U(color_0, color_1) * overlap(prop_seg[1], prop_seg[0]);
-    if (wdata.has_Dt)
+    if (wdata.has_Dt) {
+      std::vector<segment_t> seglist_temp = {prop_seg[1]};
       ln_trace_ratio += K_overlap(seglist_temp, prop_seg[0].tau_c, prop_seg[0].tau_cdag, wdata.K, color_0, color_1);
+    }
 
     double trace_ratio = std::exp(ln_trace_ratio);
 
@@ -116,24 +116,24 @@ namespace triqs_ctseg::moves {
           return 0;
         }
         if (prop_seg[0].tau_cdag == prop_seg[1].tau_cdag or prop_seg[0].tau_c == prop_seg[1].tau_c) {
-          LOG("Two inserting segments have same times. Rejecting.");
+          LOG("Two identical times in inserted segments. Rejecting.");
           return 0;
         }
       }
-      // Sort smaller and larger times for tau_cdag and tau_c
+      // Find the ordering of the two tau_cdag and the two tau_c to be inserted. They must be inserted in the correct (inceasing) order. 
       auto tau_cdag_lower = prop_seg[0].tau_cdag > prop_seg[1].tau_cdag ? prop_seg[1].tau_cdag : prop_seg[0].tau_cdag;
       auto tau_cdag_upper = prop_seg[0].tau_cdag < prop_seg[1].tau_cdag ? prop_seg[1].tau_cdag : prop_seg[0].tau_cdag;
       auto tau_c_lower    = prop_seg[0].tau_c    > prop_seg[1].tau_c    ? prop_seg[1].tau_c    : prop_seg[0].tau_c;
       auto tau_c_upper    = prop_seg[0].tau_c    < prop_seg[1].tau_c    ? prop_seg[1].tau_c    : prop_seg[0].tau_c;
-      // Determine block indices corresponding to lower and upper tau values
-      auto bl_cdag_lower  = prop_seg[0].tau_cdag > prop_seg[1].tau_cdag ? bl_idx_1             : bl_idx_0;
-      auto bl_cdag_upper  = prop_seg[0].tau_cdag < prop_seg[1].tau_cdag ? bl_idx_1             : bl_idx_0;
-      auto bl_c_lower     = prop_seg[0].tau_c    > prop_seg[1].tau_c    ? bl_idx_1             : bl_idx_0;
-      auto bl_c_upper     = prop_seg[0].tau_c    < prop_seg[1].tau_c    ? bl_idx_1             : bl_idx_0;
+      // Determine orbital indices (inside the block) corresponding to lower and upper tau values
+      auto idx_cdag_lower  = prop_seg[0].tau_cdag > prop_seg[1].tau_cdag ? bl_idx_1             : bl_idx_0;
+      auto idx_cdag_upper  = prop_seg[0].tau_cdag < prop_seg[1].tau_cdag ? bl_idx_1             : bl_idx_0;
+      auto idx_c_lower     = prop_seg[0].tau_c    > prop_seg[1].tau_c    ? bl_idx_1             : bl_idx_0;
+      auto idx_c_upper     = prop_seg[0].tau_c    < prop_seg[1].tau_c    ? bl_idx_1             : bl_idx_0;
       det_ratio = D.try_insert2(det_lower_bound_x(D, tau_cdag_lower), det_lower_bound_x(D, tau_cdag_upper) + 1, // +1 to avoid duplicates
                                 det_lower_bound_y(D, tau_c_lower)   , det_lower_bound_y(D, tau_c_upper)    + 1, // +1 to avoid duplicates
-                                {tau_cdag_lower, bl_cdag_lower}     , {tau_cdag_upper, bl_cdag_upper}         ,
-                                {tau_c_lower   , bl_c_lower}        , {tau_c_upper   , bl_c_upper}            );
+                                {tau_cdag_lower, idx_cdag_lower}     , {tau_cdag_upper, idx_cdag_upper}         ,
+                                {tau_c_lower   , idx_c_lower}        , {tau_c_upper   , idx_c_upper}            );
     }
     else { // insert on different blocks
       auto &D_0 = wdata.dets[bl_0];

benchmark/move_double/plot.ipynb

@@ -24,7 +24,7 @@ namespace triqs_ctseg::moves {
 
     LOG("\n =================== ATTEMPT DOUBLE REMOVE ================ \n");
 
-    // ------------ Choice of segment --------------
+    // ------------ Choice of segments --------------
     // Select removal colors
     int color_0 = rng(config.n_color());
     int color_1 = rng(config.n_color() - 1);
@@ -59,7 +59,6 @@ namespace triqs_ctseg::moves {
 
     // ------------  Trace ratio  -------------
     // Same as double insert, up to the sign
-    // FIXME : pull it out ?
     double ln_trace_ratio = 0.0;
     for (auto const &[i, color] : itertools::enumerate(colors)) {
       ln_trace_ratio += -wdata.mu(color) * prop_seg[i].length();
@@ -71,18 +70,16 @@ namespace triqs_ctseg::moves {
       if (wdata.has_Dt) ln_trace_ratio -= real(wdata.K(double(prop_seg[i].length()))(color, color));
     } // color
 
-    // Counting the overlap between the removal segments
-    // Make the prop_seg[1] as a seglist to use K_overlap()
-    // Be careful to the sign here!
-    std::vector<segment_t> seglist_temp = {prop_seg[1]};
+    // Overlap between the removed segments  
     ln_trace_ratio += -wdata.U(color_0, color_1) * overlap(prop_seg[1], prop_seg[0]);
-    if (wdata.has_Dt)
+    if (wdata.has_Dt) {
+      std::vector<segment_t> seglist_temp = {prop_seg[1]}; // trick for using K_overlap on a single segment
       ln_trace_ratio += K_overlap(seglist_temp, prop_seg[0].tau_c, prop_seg[0].tau_cdag, wdata.K, color_0, color_1);
+    }
 
     double trace_ratio = std::exp(ln_trace_ratio);
 
     // ------------  Det ratio  ---------------
-    // same code as in double insert. In Insert, it is a true bound, does not insert at same time
     auto &bl_0       = wdata.block_number[color_0];
     auto &bl_1       = wdata.block_number[color_1];
     is_same_block = bl_0 == bl_1;