implement double insert / remove segments

Author: mili1247

benchmark/dimer/dimer.out.h5

@@ -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/dimer/dimer.py

@@ -0,0 +1,73 @@
+# Single orbital with dynamical spin-spin interactions. 
+# 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
+from triqs.operators import n
+import h5
+from triqs.utility.h5diff import h5diff
+from triqs_ctseg import SolverCore as Solver
+
+# Numerical values
+beta = 10
+U = 4.0
+mu = U/2
+J = 0.5
+n_tau = 2051
+n_tau_bosonic = 2001
+
+# Solver construction parameters
+gf_struct = [('down', 1), ('up', 1)]
+constr_params = {
+    "gf_struct": gf_struct,
+    "beta": beta,
+    "n_tau": n_tau,
+    "n_tau_bosonic": n_tau_bosonic
+}
+
+# Construct solver
+S = Solver(**constr_params)
+
+# Get inputs from reference file
+with h5.HDFArchive("ctint.ref.h5", 'r') as Af:
+    g0 = Af["dmft_loop/i_001/S/G0_iw/up"]
+    Q_tau = Af["dmft_loop/i_000/Q_tau"]
+
+# Hybridization Delta(tau)
+n_iw = 1025
+Delta = GfImFreq(indices=[0], beta=beta, n_points=n_iw)
+invg0 = GfImFreq(indices=[0], beta=beta, n_points=n_iw)
+invg0 << inverse(g0)
+Delta << iOmega_n + mu - invg0
+S.Delta_tau << Fourier(Delta)
+
+# Spin-spin interaction (D0(tau) and Jperp(tau))
+S.Jperp_tau << -J**2*Q_tau
+S.D0_tau["up", "up"] << -0.25*J**2*Q_tau
+S.D0_tau["down", "down"] << -0.25*J**2*Q_tau
+S.D0_tau["up", "down"] << 0.25*J**2*Q_tau
+S.D0_tau["down", "up"] << 0.25*J**2*Q_tau
+
+# Solve parameters
+solve_params = {
+    "h_int": U*n("up", 0)*n("down", 0),
+    "h_loc0": -mu * (n("up", 0) + n("down", 0)),
+    "length_cycle": 50,
+    "n_warmup_cycles": 1000,
+    "n_cycles": 1000000,
+    "measure_F_tau": True,
+    "measure_nn_tau": True,
+    "measure_nn_static": True
+    }
+
+# Solve
+S.solve(**solve_params)
+
+# 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/dimer/dimer_pyed.ref.h5

@@ -19,6 +19,8 @@
 
 #include "./moves/insert_segment.hpp"
 #include "./moves/remove_segment.hpp"
+#include "./moves/double_insert_segment.hpp"
+#include "./moves/double_remove_segment.hpp"
 #include "./moves/split_segment.hpp"
 #include "./moves/regroup_segment.hpp"
 #include "./moves/move_segment.hpp"