Skip to content

Refactor dc_solver to use spin_kind_e and fix spin-polarized double counting #6

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jaemolihm wants to merge 1 commit into
base: unstable
Choose a base branch
from
Open

Refactor dc_solver to use spin_kind_e and fix spin-polarized double counting #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
168 changes: 68 additions & 100 deletions c++/triqs_modest/double_counting.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -44,129 +42,99 @@ namespace triqs::modest {
throw std::runtime_error("Not implemented! Complain to the developers");
}
}
//------------------------------------------------------------------------------------

std::pair<nda::array<nda::matrix<double>, 2>, nda::matrix<double>>
//double_counting_from_gf(block2_gf<imfreq, matrix_valued> const &GC, double U_int, double J_hund,
double_counting(nda::array<nda::matrix<dcomplex>, 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<nda::matrix<double>, 2>(n_atoms, n_sigma);
auto E_DC = nda::matrix<double>(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<double>(n_orb);
}

for (auto atom : range(n_atoms)) {
auto Nsigma = nda::zeros<double>(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<long>(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<double, nda::vector<double>> get_total_density(nda::array<nda::matrix<dcomplex>, 2> const &density_matrix) {
auto [n_blocks, n_sigma] = density_matrix.shape();
auto N_per_sigma = nda::zeros<double>(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<double, nda::vector<double>> get_total_density(spin_kind_e spin_kind,
nda::array<nda::matrix<dcomplex>, 2> const &density_matrix) {
auto [n_blocks, n_sig] = density_matrix.shape();
auto N_per_sigma = nda::zeros<double>(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<nda::matrix<double>, 2> double_counting_sigma_dc(nda::array<nda::matrix<dcomplex>, 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<nda::matrix<dcomplex>, 2> dc_solver::get_density_matrix_from_gf(block_gf<imfreq, matrix_valued> const &gimp) {
auto n_blocks = gimp.size() / n_sigma;
auto density_matrix = nda::array<nda::matrix<dcomplex>, 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<nda::matrix<dcomplex>> dc_solver::dc_self_energy(nda::array<nda::matrix<dcomplex>, 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<std::vector>(), 0, std::plus<>());
n_orb = (n_sigma == 2) ? n_orb : static_cast<long>(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<nda::matrix<double>, 2>(n_blocks, n_sigma);
auto Sigma_DC = std::vector<nda::matrix<dcomplex>>(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<double>(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<dcomplex>(density_matrix(bl, sigma).shape()[0]);
return Sigma_DC;
}

//------------------------------------------------------------------------------------
double dc_solver::dc_energy(nda::array<nda::matrix<dcomplex>, 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> double_counting_energy_dc(nda::array<nda::matrix<dcomplex>, 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<double>(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<std::vector>(), 0, std::plus<>());
n_orb = (n_sigma == 2) ? n_orb : static_cast<long>(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<double>(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<nda::matrix<double>, 2> dc_solver::get_density_matrix_from_gf(block_gf<imfreq, matrix_valued> const &gimp) {
auto n_blocks = gimp.size() / n_sigma;
auto density_matrix = nda::array<nda::matrix<double>, 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<nda::matrix<dcomplex>> dc_solver::dc_self_energy(block_gf<imfreq, matrix_valued> const &gimp) {
auto Sigma_DC = std::vector<nda::matrix<dcomplex>>(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<double> dc_solver::dc_energy(block_gf<imfreq, matrix_valued> 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<imfreq, matrix_valued> const &gimp) {
return dc_energy(get_density_matrix_from_gf(gimp));
}

} // namespace triqs::modest
46 changes: 27 additions & 19 deletions c++/triqs_modest/double_counting.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
#include <nda/nda.hpp>
#include "utils/defs.hpp"
#include "triqs/gfs.hpp"
#include "local_space.hpp"

using namespace triqs::gfs;

Expand All @@ -26,17 +27,6 @@ namespace triqs::modest {
*/
std::pair<double, double> 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<nda::array<nda::matrix<double>, 2>, nda::matrix<double>> double_counting(nda::array<nda::matrix<dcomplex>, 2> const &density_matrix,
double U_int, double J_hund, std::string const method);

/**
* @ingroup double_counting
Expand All @@ -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;
Expand All @@ -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<nda::matrix<double>, 2> get_density_matrix_from_gf(block_gf<imfreq, matrix_valued> const &gimp);
nda::array<nda::matrix<dcomplex>, 2> get_density_matrix_from_gf(block_gf<imfreq, matrix_valued> 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.
Expand All @@ -87,13 +79,29 @@ namespace triqs::modest {
std::vector<nda::matrix<dcomplex>> dc_self_energy(block_gf<imfreq, matrix_valued> 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<imfreq, matrix_valued> 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<nda::matrix<dcomplex>> dc_self_energy(nda::array<nda::matrix<dcomplex>, 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<double> dc_energy(block_gf<imfreq, matrix_valued> const &gimp);
double dc_energy(nda::array<nda::matrix<dcomplex>, 2> const &density_matrix);
};

} // namespace triqs::modest
12 changes: 11 additions & 1 deletion c++/triqs_modest/local_space.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand Down Expand Up @@ -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<std::string> sigma_names() const {
Expand Down
1 change: 1 addition & 0 deletions python/triqs_modest/utils/dc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#include <c2py/c2py.hpp>

#include "triqs_modest/double_counting.hpp"
#include "triqs_modest/local_space.hpp"

#include "dc.wrap.cxx"
2 changes: 1 addition & 1 deletion python/triqs_modest/utils/dc.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading