translate anderson and jperp tests into python

Author: nkavokine

c++/triqs_ctseg/work_data.cpp

@@ -52,6 +52,7 @@ namespace triqs_ctseg {
 
     // Print block/index/color correspondence
     if (c.rank() == 0) {
+      spdlog::info("\n");
       for (auto const &color : range(n_color)) {
         spdlog::info("Block: {}    Index: {}    Color: {}", gf_struct[block_number[color]].first, index_in_block[color],
                      color);
@@ -73,8 +74,8 @@ namespace triqs_ctseg {
 
     // Report
     if (c.rank() == 0) {
-      spdlog::info("Interaction matrix: U = {}", U);
-      spdlog::info("Orbital energies: mu - eps = {}", mu);
+      spdlog::info("\n Interaction matrix: U = {} \n", U);
+      spdlog::info("Orbital energies: mu - eps = {} \n", mu);
     }
 
     // Dynamical interactions: convert Block2Gf to matrix Gf of size n_colors
@@ -164,7 +165,7 @@ namespace triqs_ctseg {
 
     // Report
     if (c.rank() == 0) {
-      spdlog::info("Dynamical interactions = {}, Jperp interactions = {}", has_Dt, has_Jperp);
+      spdlog::info("Dynamical interactions = {}, Jperp interactions = {} \n", has_Dt, has_Jperp);
       if (p.measure_F_tau and !rot_inv)
         spdlog::info("WARNING: Cannot measure F(tau) because spin-spin interaction is not rotationally invariant.");
     }

test/c++/Jperp.cpp

@@ -6,6 +6,9 @@ endforeach()
 
 # List of all tests
 set(all_tests
+  ./anderson
+  ./dynamical_U
+  ./Jperp
   ./spin_spin
   ./two_orbitals
   ./multiorb

test/c++/Jperp.ref.h5

@@ -0,0 +1,61 @@
+# Single orbital with Jperp interaction and no hybridization 
+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
+eps = 0.2
+n_tau = 2051
+n_tau_bosonic = 2001
+wp = 1
+g = 4
+
+# 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)
+
+# Jperp interaction
+d0_iw = GfImFreq(indices=[0], beta=beta, statistic="Boson", n_points=n_tau_bosonic//2)
+d0_tau = GfImTime(indices=[0], beta=beta, statistic="Boson", n_points=n_tau_bosonic)
+d0_iw << Function(lambda w: g * wp**2 / (w**2 - wp**2))
+d0_tau << Fourier(d0_iw)
+S.Jperp_tau << d0_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": 10000,
+    "measure_nn_tau": True,
+    "measure_nn_static": True
+    }
+
+# Solve
+S.solve(**solve_params)
+
+# Save and compare to reference
+if mpi.is_master_node():
+    with h5.HDFArchive("Jperp.out.h5", 'w') as A:
+        A['G_tau'] = S.results.G_tau
+        A['nn_tau'] = S.results.nn_tau
+        A['nn'] = S.results.nn_static
+        A['densities'] = S.results.densities
+
+    h5diff("Jperp.out.h5", "Jperp.ref.h5", precision=1e-9) 

test/c++/dynamical_U.cpp

@@ -0,0 +1,62 @@
+# Single orbital with static interaction at half-filling. 
+from triqs.gf import *
+import triqs.utility.mpi as mpi
+from triqs.gf.descriptors import Function
+from triqs.operators import n
+from triqs.gf import iOmega_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
+eps = 0.2
+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)
+
+# Hybridization Delta(tau)
+iw_mesh = MeshImFreq(beta, 'Fermion', n_tau//2)
+delta_iw = GfImFreq(indices = [0], mesh = iw_mesh)
+delta_iw << inverse(iOmega_n - eps)
+S.Delta_tau["up"] << Fourier(delta_iw)
+S.Delta_tau["down"] << Fourier(delta_iw)
+
+# 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": 10000,
+    "measure_F_tau": True,
+    "measure_nn_tau": True,
+    "measure_nn_static": True
+    }
+
+# Solve
+S.solve(**solve_params)
+
+# Save and compare to reference
+if mpi.is_master_node():
+    with h5.HDFArchive("anderson.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
+
+    h5diff("anderson.out.h5", "anderson.ref.h5", precision=1e-9)