diff --git a/c++/triqs_ctseg/work_data.cpp b/c++/triqs_ctseg/work_data.cpp index 2c6b3a1..4d3c027 100644 --- a/c++/triqs_ctseg/work_data.cpp +++ b/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."); } diff --git a/test/c++/Jperp.cpp b/test/c++/Jperp.cpp deleted file mode 100644 index 23d0b24..0000000 --- a/test/c++/Jperp.cpp +++ /dev/null @@ -1,106 +0,0 @@ -// Copyright (c) 2022-2024 Simons Foundation -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You may obtain a copy of the License at -// https://www.gnu.org/licenses/gpl-3.0.txt -// -// Authors: Nikita Kavokine, Olivier Parcollet, Nils Wentzell - -#include -#include -#include - -using triqs::operators::n; -using namespace triqs_ctseg; - -TEST(CTSEG, Jperp) { - - mpi::communicator c; // Start the mpi - - double beta = 10.0; - double U = 4.0; - double mu = 2.0; - double epsilon = 0.3; - int n_cycles = 10000; - int n_warmup_cycles = 1000; - int length_cycle = 50; - int random_seed = 23488; - int n_iw = 5000; - int n_tau = 10001; - double precision = 1.e-12; - - // Prepare the parameters - constr_params_t param_constructor; - solve_params_t param_solve; - - // Construction parameters - param_constructor.beta = beta; - param_constructor.gf_struct = {{"up", 1}, {"down", 1}}; - param_constructor.n_tau = n_tau; - param_constructor.n_tau_bosonic = n_tau; - - // Create solver instance - solver_core Solver(param_constructor); - - // Solve parameters - param_solve.h_int = U * n("up", 0) * n("down", 0); - param_solve.h_loc0 = -mu * (n("up", 0) + n("down", 0)); - param_solve.n_cycles = n_cycles; - param_solve.n_warmup_cycles = n_warmup_cycles; - param_solve.length_cycle = length_cycle; - param_solve.random_seed = random_seed; - param_solve.measure_nn_tau = true; - param_solve.measure_nn_static = true; - - // Prepare Delta - nda::clef::placeholder<0> om_; - auto Delta_w = gf({beta, Fermion, n_iw}, {1, 1}); - auto Delta_tau = gf({beta, Fermion, param_constructor.n_tau}, {1, 1}); - Delta_w(om_) << 1.0 / (om_ - epsilon) + 1.0 / (om_ + epsilon); - Delta_tau() = fourier(Delta_w); - Solver.Delta_tau()[0] = Delta_tau; - Solver.Delta_tau()[1] = Delta_tau; - - // Prepare j_perp interaction - double l = 1.0; // electron boson coupling - double w0 = 1.0; // screening frequency - auto J0w = gf({beta, Boson, n_iw}, {1, 1}); - J0w(om_) << 4 * l * l * w0 / (om_ * om_ - w0 * w0); - Solver.Jperp_tau() = fourier(J0w); - - // Solve - Solver.solve(param_solve); - - // Save the results - if (c.rank() == 0) { - h5::file out_file("Jperp.out.h5", 'w'); - h5_write(out_file, "G_tau", Solver.results.G_tau); - h5_write(out_file, "nn_tau", Solver.results.nn_tau); - h5_write(out_file, "densities", Solver.results.densities); - } - - // Compare with reference - if (c.rank() == 0) { - h5::file ref_file("Jperp.ref.h5", 'r'); - block_gf G_tau; - block2_gf nn_tau; - std::map> densities; - h5_read(ref_file, "G_tau", G_tau); - h5_read(ref_file, "nn_tau", nn_tau); - h5_read(ref_file, "densities", densities); - EXPECT_ARRAY_NEAR(densities["up"], Solver.results.densities.value()["up"], precision); - EXPECT_ARRAY_NEAR(densities["down"], Solver.results.densities.value()["down"], precision); - EXPECT_BLOCK_GF_NEAR(G_tau, Solver.results.G_tau, precision); - EXPECT_BLOCK2_GF_NEAR(nn_tau, Solver.results.nn_tau.value(), precision); - } -} -MAKE_MAIN; diff --git a/test/c++/Jperp.ref.h5 b/test/c++/Jperp.ref.h5 deleted file mode 100644 index a5ac311..0000000 Binary files a/test/c++/Jperp.ref.h5 and /dev/null differ diff --git a/test/c++/dynamical_U.cpp b/test/c++/dynamical_U.cpp deleted file mode 100644 index 50f6c97..0000000 --- a/test/c++/dynamical_U.cpp +++ /dev/null @@ -1,115 +0,0 @@ -// Copyright (c) 2022-2024 Simons Foundation -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You may obtain a copy of the License at -// https://www.gnu.org/licenses/gpl-3.0.txt -// -// Authors: Nikita Kavokine, Olivier Parcollet, Nils Wentzell - -#include -#include -#include - -using triqs::operators::n; -using namespace triqs_ctseg; - -TEST(CTSEG, Dynamical_U) { - - mpi::communicator c; // Start the mpi - - double beta = 20.0; - double U = 2.0; - double mu = 1.0; - double epsilon = 0.2; - int n_cycles = 10000; - int n_warmup_cycles = 1000; - int length_cycle = 50; - int random_seed = 23488; - int n_iw = 5000; - int n_tau = 10001; - double precision = 1.e-13; - - // Prepare the parameters - constr_params_t param_constructor; - solve_params_t param_solve; - - // Construction parameters - param_constructor.beta = beta; - param_constructor.gf_struct = {{"up", 1}, {"down", 1}}; - param_constructor.n_tau = n_tau; - param_constructor.n_tau_bosonic = n_tau; - - // Create solver instance - solver_core Solver(param_constructor); - - // Solve parameters - param_solve.h_int = U * n("up", 0) * n("down", 0); - param_solve.h_loc0 = -mu * (n("up", 0) + n("down", 0)); - param_solve.n_cycles = n_cycles; - param_solve.n_warmup_cycles = n_warmup_cycles; - param_solve.length_cycle = length_cycle; - param_solve.random_seed = random_seed; - param_solve.measure_F_tau = true; - param_solve.measure_nn_tau = true; - param_solve.measure_nn_static = true; - - // Prepare Delta - nda::clef::placeholder<0> om_; - auto Delta_w = gf({beta, Fermion, n_iw}, {1, 1}); - auto Delta_tau = gf({beta, Fermion, param_constructor.n_tau}, {1, 1}); - Delta_w(om_) << 1.0 / (om_ - epsilon); - Delta_tau() = fourier(Delta_w); - Solver.Delta_tau()[0] = Delta_tau; - Solver.Delta_tau()[1] = Delta_tau; - - // Prepare dynamical interaction - double l = 1.0; // electron boson coupling - double w0 = 1.0; // screening frequency - auto D0w = gf({beta, Boson, n_iw}, {1, 1}); - auto D0t = gf({beta, Boson, param_constructor.n_tau}, {1, 1}); - D0w(om_) << 2 * l * l * w0 / (om_ * om_ - w0 * w0); - D0t() = fourier(D0w); - Solver.D0_tau()(0, 0) = D0t; - Solver.D0_tau()(0, 1) = D0t; - Solver.D0_tau()(1, 0) = D0t; - Solver.D0_tau()(1, 1) = D0t; - - // Solve - Solver.solve(param_solve); - - // Save the results - if (c.rank() == 0) { - h5::file out_file("dynamical_U.out.h5", 'w'); - h5_write(out_file, "G_tau", Solver.results.G_tau); - h5_write(out_file, "F_tau", Solver.results.F_tau); - h5_write(out_file, "nn_tau", Solver.results.nn_tau); - h5_write(out_file, "densities", Solver.results.densities); - } - - // Compare with reference - if (c.rank() == 0) { - h5::file ref_file("dynamical_U.ref.h5", 'r'); - block_gf G_tau, F_tau; - block2_gf nn_tau; - std::map> densities; - h5_read(ref_file, "G_tau", G_tau); - h5_read(ref_file, "F_tau", F_tau); - h5_read(ref_file, "nn_tau", nn_tau); - h5_read(ref_file, "densities", densities); - EXPECT_ARRAY_NEAR(densities["up"], Solver.results.densities.value()["up"], precision); - EXPECT_ARRAY_NEAR(densities["down"], Solver.results.densities.value()["down"], precision); - EXPECT_BLOCK_GF_NEAR(G_tau, Solver.results.G_tau, precision); - EXPECT_BLOCK_GF_NEAR(F_tau, Solver.results.F_tau.value(), precision); - EXPECT_BLOCK2_GF_NEAR(nn_tau, Solver.results.nn_tau.value(), precision); - } -} -MAKE_MAIN; diff --git a/test/c++/dynamical_U.ref.h5 b/test/c++/dynamical_U.ref.h5 deleted file mode 100644 index f8b4dc6..0000000 Binary files a/test/c++/dynamical_U.ref.h5 and /dev/null differ diff --git a/test/python/CMakeLists.txt b/test/python/CMakeLists.txt index 8ac575b..3cf1990 100644 --- a/test/python/CMakeLists.txt +++ b/test/python/CMakeLists.txt @@ -6,6 +6,9 @@ endforeach() # List of all tests set(all_tests + ./anderson + ./dynamical_U + ./Jperp ./spin_spin ./two_orbitals ./multiorb diff --git a/test/python/Jperp.py b/test/python/Jperp.py new file mode 100644 index 0000000..99e7d77 --- /dev/null +++ b/test/python/Jperp.py @@ -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) diff --git a/test/python/Jperp.ref.h5 b/test/python/Jperp.ref.h5 new file mode 100644 index 0000000..0ee9b8d Binary files /dev/null and b/test/python/Jperp.ref.h5 differ diff --git a/test/python/anderson.py b/test/python/anderson.py new file mode 100644 index 0000000..53f5fe9 --- /dev/null +++ b/test/python/anderson.py @@ -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) diff --git a/test/python/anderson.ref.h5 b/test/python/anderson.ref.h5 new file mode 100644 index 0000000..967c5fe Binary files /dev/null and b/test/python/anderson.ref.h5 differ diff --git a/test/python/dynamical_U.py b/test/python/dynamical_U.py new file mode 100644 index 0000000..22ec9bd --- /dev/null +++ b/test/python/dynamical_U.py @@ -0,0 +1,74 @@ +# Single orbital with dynamical interaction 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 +wp = 1 +g = 1 + +# 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) + +# Dynamical 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.D0_tau["up", "up"] << d0_tau +S.D0_tau["down", "down"] << d0_tau +S.D0_tau["up", "down"] << d0_tau +S.D0_tau["down", "up"] << 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_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("dynamical_U.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("dynamical_U.out.h5", "dynamical_U.ref.h5", precision=1e-9) diff --git a/test/python/dynamical_U.ref.h5 b/test/python/dynamical_U.ref.h5 new file mode 100644 index 0000000..31d9270 Binary files /dev/null and b/test/python/dynamical_U.ref.h5 differ