From 6477b8eb9221ef85a8097ee7400aaede2b211103 Mon Sep 17 00:00:00 2001 From: Jae-Mo Lihm Date: Mon, 9 Feb 2026 17:14:01 -0500 Subject: [PATCH] Refactor dc_solver to use spin_kind_e and fix spin-polarized double counting The dc_solver constructor now takes spin_kind_e instead of n_sigma, enabling correct spin-dependent behavior for Polarized systems. Previously, per-spin densities were unconditionally overridden to N_total/2, making sFLL/sAMF identical to cFLL/cAMF. - Add density matrix overloads for dc_self_energy and dc_energy - Change dc_energy return type from nda::matrix to double - Add n_sigma_from_spin_kind free function to local_space.hpp - Register spin_kind_e enum in the dc Python module - Remove redundant free functions (double_counting, double_counting_sigma_dc, double_counting_energy_dc) - Add hardcoded-value tests for dc_formulas, density matrix overloads, Polarized vs NonPolarized, and NonColinear --- c++/triqs_modest/double_counting.cpp | 168 +++++++++++---------------- c++/triqs_modest/double_counting.hpp | 46 +++++--- c++/triqs_modest/local_space.hpp | 12 +- python/triqs_modest/utils/dc.cpp | 1 + python/triqs_modest/utils/dc.toml | 2 +- test/c++/dc_solver_test.cpp | 134 ++++++++++++++++++++- test/c++/embed_desc_test.cpp | 2 +- 7 files changed, 241 insertions(+), 124 deletions(-) diff --git a/c++/triqs_modest/double_counting.cpp b/c++/triqs_modest/double_counting.cpp index cff6dad..0494a20 100644 --- a/c++/triqs_modest/double_counting.cpp +++ b/c++/triqs_modest/double_counting.cpp @@ -33,8 +33,6 @@ namespace triqs::modest { 0.5 * U * N_tot * N_tot - 0.25 * ((U + 2.0 * L_orbit * J) / (2.0 * L_orbit + 1.0)) * N_tot * N_tot // Energy - 0.25 * ((U + 2.0 * L_orbit * J) / (2.0 * L_orbit + 1.0)) * Mag * Mag}; } else if (method == "cHeld") { - // double C = double(n_orb) - 1.0; - // double Umean = (U + ((U - 2.0 * J) + (U - 3.0 * J)) * C) / (2.0 * double(n_orb) - 1); double U_mean = (U + (n_orb - 1) * (U - 2 * J) + (n_orb - 1) * (U - 3 * J)) / (2 * n_orb - 1); return { (U_mean * (N_tot - 0.5)), //DC @@ -44,129 +42,99 @@ namespace triqs::modest { throw std::runtime_error("Not implemented! Complain to the developers"); } } - //------------------------------------------------------------------------------------ - - std::pair, 2>, nda::matrix> - //double_counting_from_gf(block2_gf const &GC, double U_int, double J_hund, - double_counting(nda::array, 2> const &density_matrix, double U_int, double J_hund, std::string const method) { - - auto [n_atoms, n_sigma] = density_matrix.shape(); - auto Sigma_DC = nda::array, 2>(n_atoms, n_sigma); - auto E_DC = nda::matrix(n_atoms, n_sigma); - - //for (auto [atom, sigma] : indices(GC)) { - for (auto atom : nda::range(n_atoms)) - for (auto sigma : nda::range(n_sigma)) { - //auto n_orb = GC(atom, sigma).target_shape()[0]; - auto n_orb = density_matrix(atom, sigma).shape()[0]; - Sigma_DC(atom, sigma) = nda::eye(n_orb); - } - - for (auto atom : range(n_atoms)) { - auto Nsigma = nda::zeros(n_sigma); - //for (auto sigma : range(n_sigma)) { Nsigma(sigma) += real(trace(density(GC(atom, sigma)))); } - for (auto sigma : range(n_sigma)) { Nsigma(sigma) += real(trace(density_matrix(atom, sigma))); } - auto Ntot = stdr::fold_left(Nsigma, 0, std::plus<>()); - for (auto sigma : range(n_sigma)) { - if (n_sigma == 2) { - Nsigma(sigma) = Ntot / n_sigma; - } else if (n_sigma == 1) { - Nsigma(sigma) = 0.5 * Ntot; - } - } - - for (auto sigma : range(n_sigma)) { - auto n_orb = density_matrix(atom, sigma).shape()[0]; - n_orb = (n_sigma == 2) ? n_orb : static_cast(n_orb / 2); - auto [DC_val, E_val] = dc_formulas(method, Ntot, Nsigma(sigma), n_orb, U_int, J_hund); - Sigma_DC(atom, sigma) *= DC_val; - E_DC(atom, sigma) = E_val; - } - } - return {Sigma_DC, E_DC}; - } //------------------------------------------------------------------------------------ - std::pair> get_total_density(nda::array, 2> const &density_matrix) { - auto [n_blocks, n_sigma] = density_matrix.shape(); - auto N_per_sigma = nda::zeros(n_sigma); - for (auto sigma : range(n_sigma)) { - for (auto bl : range(n_blocks)) { N_per_sigma(sigma) += real(trace(density_matrix(bl, sigma))); } - } - auto N_total = stdr::fold_left(N_per_sigma, 0, std::plus<>()); - for (auto sigma : range(n_sigma)) { - if (n_sigma == 2) { - N_per_sigma(sigma) = N_total / n_sigma; - } else if (n_sigma == 1) { - N_per_sigma(sigma) = 0.5 * N_total; - } + // Compute total density N_total and per-spin densities N_per_sigma from a density matrix. + // + // The spin_kind determines how per-spin densities are treated: + // - Polarized: N_per_sigma retains the actual per-spin traces from the density matrix. + // This allows spin-dependent DC formulas (sFLL, sAMF) to produce + // different corrections for each spin channel. + // - NonPolarized: N_per_sigma is overridden to N_total/2 for both spins, making + // spin-dependent formulas reduce to their charge-only counterparts. + // - NonColinear: Single spin channel (n_sig=1); N_per_sigma is set to N_total/2 + // so that dc_formulas sees the per-spin density as half the total. + static std::pair> get_total_density(spin_kind_e spin_kind, + nda::array, 2> const &density_matrix) { + auto [n_blocks, n_sig] = density_matrix.shape(); + auto N_per_sigma = nda::zeros(n_sig); + for (auto sigma : range(n_sig)) + for (auto bl : range(n_blocks)) N_per_sigma(sigma) += real(trace(density_matrix(bl, sigma))); + auto N_total = stdr::fold_left(N_per_sigma, 0.0, std::plus<>()); + + switch (spin_kind) { + case spin_kind_e::Polarized: + break; + case spin_kind_e::NonPolarized: + for (auto sigma : range(n_sig)) N_per_sigma(sigma) = 0.5 * N_total; + break; + case spin_kind_e::NonColinear: + N_per_sigma(0) = 0.5 * N_total; + break; } return {N_total, N_per_sigma}; } - //------------------------------------------------------------------------------------ - nda::array, 2> double_counting_sigma_dc(nda::array, 2> const &density_matrix, double U_int, double J_hund, - std::string const method) { + //------------------------------------------------------------------------------------ + dc_solver::dc_solver(spin_kind_e spin_kind, std::string method, double U_int, double J_hund) + : _spin_kind{spin_kind}, n_sigma{n_sigma_from_spin_kind(spin_kind)}, method{std::move(method)}, U_int{U_int}, J_hund{J_hund} {}; - auto [N_total, N_per_sigma] = get_total_density(density_matrix); + //------------------------------------------------------------------------------------ + nda::array, 2> dc_solver::get_density_matrix_from_gf(block_gf const &gimp) { + auto n_blocks = gimp.size() / n_sigma; + auto density_matrix = nda::array, 2>(n_blocks, n_sigma); + for (auto bl : range(n_blocks)) + for (auto sigma : range(n_sigma)) density_matrix(bl, sigma) = density(gimp[bl + sigma * n_blocks]); + return density_matrix; + } - auto [n_blocks, n_sigma] = density_matrix.shape(); + //------------------------------------------------------------------------------------ + std::vector> dc_solver::dc_self_energy(nda::array, 2> const &density_matrix) { + auto [N_total, N_per_sigma] = get_total_density(_spin_kind, density_matrix); + auto [n_blocks, n_sig] = density_matrix.shape(); - auto n_orb = - stdr::fold_left(density_matrix(r_all, 0) | stdv::transform([](auto x) { return x.shape()[0]; }) | tl::to(), 0, std::plus<>()); - n_orb = (n_sigma == 2) ? n_orb : static_cast(n_orb / 2); + long n_orb = 0; + for (auto bl : range(n_blocks)) n_orb += density_matrix(bl, 0).shape()[0]; + if (_spin_kind == spin_kind_e::NonColinear) n_orb /= 2; - auto Sigma_DC = nda::array, 2>(n_blocks, n_sigma); + auto Sigma_DC = std::vector>(n_blocks * n_sig); for (auto bl : range(n_blocks)) - for (auto sigma : range(n_sigma)) - Sigma_DC(bl, sigma) = std::get<0>(dc_formulas(method, N_total, N_per_sigma(sigma), n_orb, U_int, J_hund)) - * nda::eye(density_matrix(bl, sigma).shape()[0]); + for (auto sigma : range(n_sig)) + Sigma_DC[bl + sigma * n_blocks] = std::get<0>(dc_formulas(method, N_total, N_per_sigma(sigma), n_orb, U_int, J_hund)) + * nda::eye(density_matrix(bl, sigma).shape()[0]); return Sigma_DC; } + //------------------------------------------------------------------------------------ + double dc_solver::dc_energy(nda::array, 2> const &density_matrix) { + auto [N_total, N_per_sigma] = get_total_density(_spin_kind, density_matrix); + auto [n_blocks, n_sig] = density_matrix.shape(); - nda::matrix double_counting_energy_dc(nda::array, 2> const &density_matrix, double U_int, double J_hund, - std::string const method) { + long n_orb = 0; + for (auto bl : range(n_blocks)) n_orb += density_matrix(bl, 0).shape()[0]; + if (_spin_kind == spin_kind_e::NonColinear) n_orb /= 2; - auto [N_total, N_per_sigma] = get_total_density(density_matrix); + auto E_DC = std::vector(n_sig); + for (auto sigma : range(n_sig)) E_DC[sigma] = std::get<1>(dc_formulas(method, N_total, N_per_sigma(sigma), n_orb, U_int, J_hund)); - auto [n_blocks, n_sigma] = density_matrix.shape(); - auto n_orb = - stdr::fold_left(density_matrix(r_all, 0) | stdv::transform([](auto x) { return x.shape()[0]; }) | tl::to(), 0, std::plus<>()); - n_orb = (n_sigma == 2) ? n_orb : static_cast(n_orb / 2); + // E_DC must be independent of sigma (Mag appears only as Mag^2 in all current formulas). + // Guard against future formulas where this assumption may not hold. + if (n_sig > 1) + for (auto sigma : range(1, n_sig)) + if (std::abs(E_DC[sigma] - E_DC[0]) > 1.e-10) + throw std::runtime_error("dc_energy: E_DC depends on spin channel"); - auto E_DC = nda::matrix(n_blocks, n_sigma); - for (auto bl : range(n_blocks)) - for (auto sigma : range(n_sigma)) E_DC(bl, sigma) = std::get<1>(dc_formulas(method, N_total, N_per_sigma(sigma), n_orb, U_int, J_hund)); - return E_DC; + return E_DC[0]; } //------------------------------------------------------------------------------------ - dc_solver::dc_solver(long n_sigma, std::string method, double U_int, double J_hund) - : n_sigma{n_sigma}, method{std::move(method)}, U_int{U_int}, J_hund{J_hund} {}; - - nda::array, 2> dc_solver::get_density_matrix_from_gf(block_gf const &gimp) { - auto n_blocks = gimp.size() / n_sigma; - auto density_matrix = nda::array, 2>(n_blocks, n_sigma); - for (auto bl : range(n_blocks)) - for (auto sigma : range(n_sigma)) density_matrix(bl, sigma) = real(density(gimp[bl + sigma * n_blocks])); - return density_matrix; - } - //------------------------------------------------------------------------------------ - std::vector> dc_solver::dc_self_energy(block_gf const &gimp) { - auto Sigma_DC = std::vector>(gimp.size()); - auto density_matrix = get_density_matrix_from_gf(gimp); - auto n_blocks = density_matrix.extent(0); - auto Sigma_dc_mat = double_counting_sigma_dc(density_matrix, U_int, J_hund, method); - for (auto bl : range(n_blocks)) - for (auto sigma : range(n_sigma)) Sigma_DC[bl + sigma * n_blocks] = Sigma_dc_mat(bl, sigma); - return Sigma_DC; + return dc_self_energy(get_density_matrix_from_gf(gimp)); } //------------------------------------------------------------------------------------ - nda::matrix dc_solver::dc_energy(block_gf const &gimp) { - return double_counting_energy_dc(get_density_matrix_from_gf(gimp), U_int, J_hund, method); + double dc_solver::dc_energy(block_gf const &gimp) { + return dc_energy(get_density_matrix_from_gf(gimp)); } } // namespace triqs::modest diff --git a/c++/triqs_modest/double_counting.hpp b/c++/triqs_modest/double_counting.hpp index b9673b1..9709ff3 100644 --- a/c++/triqs_modest/double_counting.hpp +++ b/c++/triqs_modest/double_counting.hpp @@ -7,6 +7,7 @@ #include #include "utils/defs.hpp" #include "triqs/gfs.hpp" +#include "local_space.hpp" using namespace triqs::gfs; @@ -26,17 +27,6 @@ namespace triqs::modest { */ std::pair dc_formulas(std::string const method, double const N_tot, double const N_sigma, long const n_orb, double const U, double const J); - /** - * @brief Compute double counting correction for a DC type (method) from the density matrix of a Green's function. - * - * @param density_matrix Density matrix [orbital, spin] indices. - * @param U_int Coulomb interaction parameter. - * @param J_hund Hund's coupling interaction parameter. - * @param method DC formula (`sFLL`, `cFLL`, `sAMF`, `cAMF`, `cHeld`). - * @return A pair of \f$ \Sigma_{DC} \f$ and \f$ E_{DC} \f$. - */ - std::pair, 2>, nda::matrix> double_counting(nda::array, 2> const &density_matrix, - double U_int, double J_hund, std::string const method); /** * @ingroup double_counting @@ -49,7 +39,9 @@ namespace triqs::modest { class dc_solver { private: - // number of spins + // spin kind + spin_kind_e _spin_kind; + // number of spins (derived from spin_kind) long n_sigma; // double counting method to use std::string method; @@ -64,21 +56,21 @@ namespace triqs::modest { * @param gimp The impurity Green's function which is used to calculate the orbital-resolved density matrices to evaluate the double counting formula. * @return Density matrix of Green's function in block_matrix format. */ - nda::array, 2> get_density_matrix_from_gf(block_gf const &gimp); + nda::array, 2> get_density_matrix_from_gf(block_gf const &gimp); public: /** * @brief Construct a double counting "solver". * - * @param n_sigma Dimension of the \f$ \sigma \f$ index. + * @param spin_kind Spin kind of the correlated space (Polarized, NonPolarized, NonColinear). * @param method Double counting formula (method) to call (options: `cFLL`, `sFLL`, `cAMF`, `sAMF`, `cHeld`). * @param U_int Hubbard \f$ U \f$ to use in the DC formula. * @param J_hund Hund's coupling \f$ J \f$ to use in the DC formula. */ - dc_solver(long n_sigma, std::string method, double U_int, double J_hund); + dc_solver(spin_kind_e spin_kind, std::string method, double U_int, double J_hund); /** - * @brief Compute the double-counting self-energy. + * @brief Compute the double-counting self-energy from a Green's function. * * @param gimp The impurity Green's function which is used to calculate the orbital-resolved density matrices to * evaluate the double counting formula. @@ -87,13 +79,29 @@ namespace triqs::modest { std::vector> dc_self_energy(block_gf const &gimp); /** - * @brief Compute the double counting correction to the energy. + * @brief Compute the double counting correction to the energy from a Green's function. * * @param gimp The impurity Green's function which is used to calculate the orbital-resolved density matrices to * evaluate the double counting formula. - * @return Double counting energy term. + * @return Double counting energy \f$ E_{DC} \f$. + */ + double dc_energy(block_gf const &gimp); + + /** + * @brief Compute the double-counting self-energy from a density matrix. + * + * @param density_matrix Density matrix with shape [n_blocks, n_sigma]. + * @return Double counting self-energy term. + */ + std::vector> dc_self_energy(nda::array, 2> const &density_matrix); + + /** + * @brief Compute the double counting correction to the energy from a density matrix. + * + * @param density_matrix Density matrix with shape [n_blocks, n_sigma]. + * @return Double counting energy \f$ E_{DC} \f$. */ - nda::matrix dc_energy(block_gf const &gimp); + double dc_energy(nda::array, 2> const &density_matrix); }; } // namespace triqs::modest diff --git a/c++/triqs_modest/local_space.hpp b/c++/triqs_modest/local_space.hpp index 54a22c5..5da4350 100644 --- a/c++/triqs_modest/local_space.hpp +++ b/c++/triqs_modest/local_space.hpp @@ -21,6 +21,16 @@ namespace triqs::modest { NonColinear // σ = 0. There is no index, the spin index is grouped with other non-diagonal indices. e.g. Nambu, spin-orbit }; + /// Number of σ channels for a given spin_kind_e. + inline long n_sigma_from_spin_kind(spin_kind_e sk) { + switch (sk) { + case spin_kind_e::Polarized: return 2; + case spin_kind_e::NonPolarized: return 2; + case spin_kind_e::NonColinear: return 1; + } + __builtin_unreachable(); + } + // -------------------------------------------------------------------- /// Info on an atomic shell. struct atomic_orbs { @@ -140,7 +150,7 @@ namespace triqs::modest { [[nodiscard]] spin_kind_e spin_kind() const { return _spin_kind; }; /// Dimension of the \f$ \sigma \f$ index. - [[nodiscard]] long n_sigma() const { return _spin_kind == spin_kind_e::NonColinear ? 1 : 2; } + [[nodiscard]] long n_sigma() const { return n_sigma_from_spin_kind(_spin_kind); } /// Names of spin indices for naming blocks in block GFs. [[nodiscard]] std::vector sigma_names() const { diff --git a/python/triqs_modest/utils/dc.cpp b/python/triqs_modest/utils/dc.cpp index bcba58b..c99c632 100644 --- a/python/triqs_modest/utils/dc.cpp +++ b/python/triqs_modest/utils/dc.cpp @@ -1,5 +1,6 @@ #include #include "triqs_modest/double_counting.hpp" +#include "triqs_modest/local_space.hpp" #include "dc.wrap.cxx" diff --git a/python/triqs_modest/utils/dc.toml b/python/triqs_modest/utils/dc.toml index 66ba77f..2b1c552 100644 --- a/python/triqs_modest/utils/dc.toml +++ b/python/triqs_modest/utils/dc.toml @@ -23,7 +23,7 @@ namespaces = "triqs triqs::modest triqs::lattice" # if reject_names and reject_names match Qname : skip X # else keep X. -match_names = "triqs::modest::(dc_solver|dc_formulas)" +match_names = "triqs::modest::(dc_solver|dc_formulas|spin_kind_e)" # Only keep elements declared in files matching the pattern # match_files = "" # String must be a regex diff --git a/test/c++/dc_solver_test.cpp b/test/c++/dc_solver_test.cpp index c78b938..550de34 100644 --- a/test/c++/dc_solver_test.cpp +++ b/test/c++/dc_solver_test.cpp @@ -13,12 +13,12 @@ auto run_dc_solver_test_case(std::string filename, double threshold, std::string auto [_, obe] = one_body_elements_from_dft_converter(filename, threshold); auto E = make_embedding(obe.C_space); auto Gimp = E.extract(gloc(gloc_ref[0][0].mesh(), obe, 0.0)); - auto double_counting = dc_solver(E.n_sigma(), method, 2.0, 1.0); + auto double_counting = dc_solver(obe.C_space.spin_kind(), method, 2.0, 1.0); auto Sigma_DC = Gimp | stdv::transform([&double_counting](auto &x) { return double_counting.dc_self_energy(x); }) | tl::to(); auto E_DC = Gimp | stdv::transform([&double_counting](auto &x) { return double_counting.dc_energy(x); }) | tl::to(); // HL : need to loop over all of the data but the reference data has a different shape. TODO: update the reference data EXPECT_NEAR(real(Sigma_DC[0][0](0, 0)), dc_ref[0][0](0, 0), 1.e-5); - EXPECT_NEAR(E_DC[0](0, 0), Edc_ref[0], 1.e-5); + EXPECT_NEAR(E_DC[0], Edc_ref[0], 1.e-5); } //-------------------------------- @@ -89,6 +89,136 @@ TEST(dc_solver_tests, dc_wien2k_large_nu_cHeld) { //NOLINT // } //-------------------------------- +//-------------------------------- +// Hardcoded-value tests for dc_formulas +//-------------------------------- +TEST(dc_solver_tests, dc_formulas_cFLL) { // NOLINT + // N_tot=4, N_sigma=2, n_orb=5, U=4, J=1 + // DC_val = U*(N_tot-0.5) - 0.5*J*(N_tot-1) = 4*3.5 - 0.5*1*3 = 12.5 + // E_val = 0.5*U*N_tot*(N_tot-1) - 0.5*J*N_tot*(0.5*N_tot-1) = 24 - 2 = 22 + auto [DC_val, E_val] = dc_formulas("cFLL", 4.0, 2.0, 5, 4.0, 1.0); + EXPECT_NEAR(DC_val, 12.5, 1.e-12); + EXPECT_NEAR(E_val, 22.0, 1.e-12); +} + +TEST(dc_solver_tests, dc_formulas_sFLL) { // NOLINT + // N_tot=4, N_sigma=2.5, n_orb=5, U=4, J=1 + // Mag = 2*2.5 - 4 = 1 + // DC_val = U*(N_tot-0.5) - J*(N_sigma-0.5) = 14 - 2 = 12 + // E_val = 24 - 2 - 0.25*J*Mag^2 = 21.75 + auto [DC_val, E_val] = dc_formulas("sFLL", 4.0, 2.5, 5, 4.0, 1.0); + EXPECT_NEAR(DC_val, 12.0, 1.e-12); + EXPECT_NEAR(E_val, 21.75, 1.e-12); +} + +TEST(dc_solver_tests, dc_formulas_cAMF) { // NOLINT + // N_tot=4, N_sigma=2, n_orb=5, U=4, J=1 + // L_orbit = 0.5*(5-1) = 2 + // C = (U + 2*J*L_orbit) / (2*L_orbit + 1) = 8/5 = 1.6 + // DC_val = U*N_tot - 0.5*N_tot*C = 16 - 3.2 = 12.8 + // E_val = (0.5*U - 0.25*C)*N_tot^2 = 1.6*16 = 25.6 + auto [DC_val, E_val] = dc_formulas("cAMF", 4.0, 2.0, 5, 4.0, 1.0); + EXPECT_NEAR(DC_val, 12.8, 1.e-12); + EXPECT_NEAR(E_val, 25.6, 1.e-12); +} + +TEST(dc_solver_tests, dc_formulas_sAMF) { // NOLINT + // N_tot=4, N_sigma=2.5, n_orb=5, U=4, J=1 + // L_orbit = 2, C_factor = 8/5 = 1.6, Mag = 2*2.5 - 4 = 1 + // DC_val = U*N_tot - C_factor*N_sigma = 16 - 4 = 12 + // E_val = 0.5*U*N_tot^2 - 0.25*C_factor*N_tot^2 - 0.25*C_factor*Mag^2 = 32 - 6.4 - 0.4 = 25.2 + auto [DC_val, E_val] = dc_formulas("sAMF", 4.0, 2.5, 5, 4.0, 1.0); + EXPECT_NEAR(DC_val, 12.0, 1.e-12); + EXPECT_NEAR(E_val, 25.2, 1.e-12); +} + +TEST(dc_solver_tests, dc_formulas_cHeld) { // NOLINT + // N_tot=4, N_sigma=2, n_orb=5, U=4, J=1 + // U_mean = (4 + 4*(4-2) + 4*(4-3)) / 9 = (4+8+4)/9 = 16/9 + // DC_val = 16/9 * 3.5 = 56/9 + // E_val = 0.5 * 16/9 * 4 * 3 = 96/9 + auto [DC_val, E_val] = dc_formulas("cHeld", 4.0, 2.0, 5, 4.0, 1.0); + EXPECT_NEAR(DC_val, 56.0 / 9.0, 1.e-12); + EXPECT_NEAR(E_val, 96.0 / 9.0, 1.e-12); +} + +//-------------------------------- +// Density matrix overload tests +//-------------------------------- +TEST(dc_solver_tests, density_matrix_overload) { // NOLINT + // Create a density matrix with known traces: 1 block, 2 spins, 5 orbitals each + // N_up = 2.0, N_down = 2.0 => N_tot = 4.0 + auto density_matrix = nda::array, 2>(1, 2); + density_matrix(0, 0) = 0.4 * nda::eye(5); // trace = 2.0 + density_matrix(0, 1) = 0.4 * nda::eye(5); // trace = 2.0 + + auto solver = dc_solver(spin_kind_e::NonPolarized, "cFLL", 4.0, 1.0); + auto Sigma_DC = solver.dc_self_energy(density_matrix); + auto E_DC = solver.dc_energy(density_matrix); + + // NonPolarized: N_per_sigma = N_tot/2 = 2.0 for both spins + // cFLL: DC_val = 12.5, E_val = 22.0 + // vector layout: [bl + sigma * n_blocks], n_blocks=1 => [0]=sigma0, [1]=sigma1 + EXPECT_NEAR(real(Sigma_DC[0](0, 0)), 12.5, 1.e-12); + EXPECT_NEAR(real(Sigma_DC[1](0, 0)), 12.5, 1.e-12); + EXPECT_NEAR(E_DC, 22.0, 1.e-12); +} + +//-------------------------------- +// Polarized vs NonPolarized sFLL test (bug fix validation) +//-------------------------------- +TEST(dc_solver_tests, polarized_vs_nonpolarized_sFLL) { // NOLINT + // N_up=2.5, N_down=1.5 => N_tot=4.0 + auto density_matrix = nda::array, 2>(1, 2); + density_matrix(0, 0) = 0.5 * nda::eye(5); // trace = 2.5 + density_matrix(0, 1) = 0.3 * nda::eye(5); // trace = 1.5 + + // NonPolarized: overrides N_per_sigma to N_tot/2 = 2.0 for both spins + // => sFLL becomes identical to cFLL: DC_val = 12.5 for both spins + auto solver_np = dc_solver(spin_kind_e::NonPolarized, "sFLL", 4.0, 1.0); + auto Sigma_np = solver_np.dc_self_energy(density_matrix); + auto E_np = solver_np.dc_energy(density_matrix); + EXPECT_NEAR(real(Sigma_np[0](0, 0)), 12.5, 1.e-12); + EXPECT_NEAR(real(Sigma_np[1](0, 0)), 12.5, 1.e-12); + EXPECT_NEAR(E_np, 22.0, 1.e-12); + + // Polarized: preserves per-spin densities + // Spin up: N_sigma=2.5 => DC_val = U*(N_tot-0.5) - J*(N_sigma-0.5) = 14 - 2 = 12 + // Spin down: N_sigma=1.5 => DC_val = U*(N_tot-0.5) - J*(N_sigma-0.5) = 14 - 1 = 13 + auto solver_p = dc_solver(spin_kind_e::Polarized, "sFLL", 4.0, 1.0); + auto Sigma_p = solver_p.dc_self_energy(density_matrix); + auto E_p = solver_p.dc_energy(density_matrix); + EXPECT_NEAR(real(Sigma_p[0](0, 0)), 12.0, 1.e-12); + EXPECT_NEAR(real(Sigma_p[1](0, 0)), 13.0, 1.e-12); + + // Polarized sFLL energy: Mag = 2*N_up - N_tot = 1, E = 24 - 2 - 0.25*1 = 21.75 + EXPECT_NEAR(E_p, 21.75, 1.e-12); + + // Verify that Polarized and NonPolarized give DIFFERENT DC self-energies for sFLL + EXPECT_GT(std::abs(Sigma_p[0](0, 0) - Sigma_np[0](0, 0)), 0.1); +} + +//-------------------------------- +// NonColinear spin kind test +//-------------------------------- +TEST(dc_solver_tests, noncolinear_cFLL) { // NOLINT + // NonColinear: 1 spin channel, orbitals include spin => n_orb_matrix = 10 for 5 orbitals + // density_matrix has shape [1, 1] with a 10x10 matrix (spin folded into orbital index) + // trace = 4.0 => N_total = 4.0, N_per_sigma = N_total/2 = 2.0, n_orb = 10/2 = 5 + auto density_matrix = nda::array, 2>(1, 1); + density_matrix(0, 0) = 0.4 * nda::eye(10); // trace = 4.0 + + auto solver = dc_solver(spin_kind_e::NonColinear, "cFLL", 4.0, 1.0); + auto Sigma_DC = solver.dc_self_energy(density_matrix); + auto E_DC = solver.dc_energy(density_matrix); + + // cFLL with N_tot=4, N_sigma=2, n_orb=5, U=4, J=1: DC_val=12.5, E_val=22.0 + EXPECT_NEAR(real(Sigma_DC[0](0, 0)), 12.5, 1.e-12); + EXPECT_NEAR(E_DC, 22.0, 1.e-12); + // Vector has only 1 entry (1 block * 1 spin channel) + EXPECT_EQ(Sigma_DC.size(), 1u); +} + #ifdef LFS //-------------------------------- // LiV2O4 (Wien2k) diff --git a/test/c++/embed_desc_test.cpp b/test/c++/embed_desc_test.cpp index 6d9b9e3..6b6edfb 100644 --- a/test/c++/embed_desc_test.cpp +++ b/test/c++/embed_desc_test.cpp @@ -83,7 +83,7 @@ TEST(embed_desc_tests, embed_matrix) { double beta = 40.0; auto mesh = mesh::imfreq{beta, statistic_enum::Fermion, 251}; auto G = E.extract(gloc(mesh, obe, find_chemical_potential(target_density, obe, beta, "bisection", 1e-3, false))); - auto double_counting = dc_solver(E.n_sigma(), "sFLL", 2.0, 1.0); + auto double_counting = dc_solver(obe.C_space.spin_kind(), "sFLL", 2.0, 1.0); auto Sigma_static = G | stdv::transform([&double_counting](auto &x) { return double_counting.dc_self_energy(x); }) | tl::to(); auto Sigma_static_C = E.embed(Sigma_static);