Merge pull request #21 from mili1247/3.3.x

move_double implementation (to solve issue #19)

Author: nkavokine

benchmark/dimer/dimer.out.h5

@@ -0,0 +1,81 @@
+# Multi-orbital impurity with diagonal Delta(tau). 
+# Number of orbitals can be changed. 
+from triqs.gf import *
+from triqs.gf.descriptors import Function
+from triqs.utility import mpi
+from triqs.operators import n
+from h5 import *
+import numpy as np
+from triqs.utility.h5diff import h5diff
+from triqs_ctseg import Solver
+
+# Number of orbitals
+n_orb = 3
+
+# Numerical values
+beta = 20.0 # inverse temperature
+mu   = 0.2 # chemical potential
+eps = 0.3 # hybridization levels
+V = 1. # hybridization strengths
+U    = 2.0 # same orbital interaction 
+Up   = 1.0 # different orbital interaction 
+J    = 0.1 # Hund coupling
+n_tau = 1001
+n_tau_bosonic = 1001
+
+# Solver construction parameters
+gf_struct = [(f'up{(l+2)//2}', 1) if l % 2 == 0 else (f'down{(l+1)//2}', 1) for l in range(2 * n_orb)]
+constr_params = {
+    "gf_struct": gf_struct,
+    "beta": beta,
+    "n_tau": n_tau,
+    "n_tau_bosonic": n_tau_bosonic
+}
+
+# Construct solver
+S = Solver(**constr_params)
+
+# Interaction Hamiltonian
+h_int = 0
+for l in range(n_orb):
+    h_int += U * n(f'up{l+1}', 0) * n(f'down{l+1}', 0)
+for l1 in range(n_orb):
+    for l2 in range(l1+1, n_orb):
+        h_int += (Up - J) * (n(f'up{l1+1}', 0) * n(f'up{l2+1}', 0) + n(f'down{l1+1}', 0) * n(f'down{l2+1}', 0))
+        h_int += Up * (n(f'up{l1+1}', 0) * n(f'down{l2+1}', 0) + n(f'down{l1+1}', 0) * n(f'up{l2+1}', 0))
+
+# Local Hamiltonian
+h_loc0 = -mu * sum(n(f"up{i + 1}", 0) + n(f"down{i + 1}", 0) for i in range(n_orb) )
+
+# Hybridization Delta(tau)
+for name, block in S.Delta_tau:
+    Delta_iw = GfImFreq(indices=[0], beta=beta, n_points=n_tau//2)
+    Delta_iw << V**2 * inverse(iOmega_n - eps)
+    block << Fourier(Delta_iw)
+
+# Solve parameters
+solve_params = {
+    "h_int": h_int,
+    "h_loc0": h_loc0,
+    "length_cycle": 50,
+    "n_warmup_cycles": 10000,
+    "n_cycles": 2000000,
+    "measure_F_tau": True,
+    "measure_nn_tau": True,
+    "measure_nn_static": True,
+    "move_double_insert_segment": True,
+    "move_double_remove_segment": True,
+    }
+
+# Solve
+S.solve(**solve_params)
+
+# Save and compare to reference
+if mpi.is_master_node():
+    with HDFArchive("multiorb.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

@@ -0,0 +1,76 @@
+# Single orbital with dynamical spin-spin interactions. 
+# Can be compared with reference obtained with CTINT (ctint.ref.h5)
+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,
+    "move_double_insert_segment": True,
+    "move_double_remove_segment": True,
+    }
+
+# Solve
+S.solve(**solve_params)
+
+# Save and compare to reference
+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
+