Skip to content

source_estate: give Occupy's stateless smearing helpers a real interface, removing two access hacks #7983

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

Merged
mohanchen merged 1 commit into from
Sep 17, 2026
+90 −69
Merged

source_estate: give Occupy's stateless smearing helpers a real interface, removing two access hacks #7983

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
9 changes: 6 additions & 3 deletions source/source_estate/elecstate_tools.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,15 +59,18 @@ void calculate_weights(const ModuleBase::matrix& ekb,
const int band_offset = band_output.local_offset();
if (PARAM.globalv.two_fermi)
{
Occupy::iweights(nks, klist->wk, nbands, band_offset, nelec_spin[0], ekb, eferm.ef_up, wg, 0, klist->isk);
Occupy::iweights(nks, klist->wk, nbands, band_offset, nelec_spin[1], ekb, eferm.ef_dw, wg, 1, klist->isk);
Occupy::iweights(nks, klist->wk, nbands, band_offset, nelec_spin[0], ekb, eferm.ef_up, wg, 0, klist->isk,
PARAM.inp.nspin);
Occupy::iweights(nks, klist->wk, nbands, band_offset, nelec_spin[1], ekb, eferm.ef_dw, wg, 1, klist->isk,
PARAM.inp.nspin);
// ef = ( ef_up + ef_dw ) / 2.0_dp need??? mohan add 2012-04-16
// Keep independent Fermi levels for the two spin channels.
}
else
{
// A spin selector of -1 requests the combined-spin occupation path.
Occupy::iweights(nks, klist->wk, nbands, band_offset, PARAM.inp.nelec, ekb, eferm.ef, wg, -1, klist->isk);
Occupy::iweights(nks, klist->wk, nbands, band_offset, PARAM.inp.nelec, ekb, eferm.ef, wg, -1, klist->isk,
PARAM.inp.nspin);
}
}
else if (Occupy::use_gaussian_broadening)
Expand Down
35 changes: 19 additions & 16 deletions source/source_estate/occupy.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ void Occupy::decision(const std::string& name, const std::string& smearing_metho
}
else if (smearing_method == "mp3")
{
// acually any order Methfessel-Paxton method can be supported in Occupy::w1gauss()
// acually any order Methfessel-Paxton method can be supported in occupy_smearing::w1gauss()
// however the parameter is string instead of int
ModuleBase::WARNING_QUIT(
"occupy",
Expand Down Expand Up @@ -128,6 +128,8 @@ void Occupy::decision(const std::string& name, const std::string& smearing_metho
* @param wg output: weight for each k, each band.
* @param is the spin index now.
* @param isk distinguish k point belong to which spin.
* @param nspin number of spin channels; 2 selects a single spin through isk when
* is != -1, and 4 halves the spin degeneracy.
*/
void Occupy::iweights(
const int nks,
Expand All @@ -139,11 +141,12 @@ void Occupy::iweights(
double& ef,
ModuleBase::matrix& wg,
const int& is, //<- is should be -1, 0, or 1. -1 means set all spins, and 0 means spin up, 1 means spin down.
const std::vector<int>& isk)
const std::vector<int>& isk,
const int nspin)
{
assert(is < 2);
double degspin = 2.0;
if (PARAM.inp.nspin == 4) {
if (nspin == 4) {
degspin = 1.0;
}
if (is != -1) {
Expand All @@ -163,7 +166,7 @@ void Occupy::iweights(
for (int ik = 0; ik < nks; ++ik)
{
// when NSPIN=2, only calculate spin up or spin down with TWO_FERMI mode(nupdown != 0)
if (PARAM.inp.nspin == 2 && isk[ik] != is && is != -1)
if (nspin == 2 && isk[ik] != is && is != -1)
{
continue;
}
Expand Down Expand Up @@ -222,7 +225,7 @@ void Occupy::gweights(const int nks,
// Calculate the Fermi energy ef
//===============================
// call efermig
Occupy::efermig(ekb, nband, nks, nelec, wk, smearing_sigma, ngauss, ef, is, isk);
occupy_smearing::efermig(ekb, nband, nks, nelec, wk, smearing_sigma, ngauss, ef, is, isk);
demet = 0.0;

for (int ik = 0; ik < nks; ik++)
Expand All @@ -238,15 +241,15 @@ void Occupy::gweights(const int nks,
// Calculate the gaussian weights
//================================
// call wgauss
wg(ik, ib) = wk[ik] * Occupy::wgauss((ef - ekb(ik, ib)) / smearing_sigma, ngauss);
wg(ik, ib) = wk[ik] * occupy_smearing::wgauss((ef - ekb(ik, ib)) / smearing_sigma, ngauss);

//====================================================================
// The correct form of the band energy is \int e n(e) de for e<ef
// demet is the correction to add to the sum of eigenvalues
//====================================================================
// Mohan fix bug 2010-1-9
// call w1gauss
demet += wk[ik] * smearing_sigma * Occupy::w1gauss((ef - ekb(ik, ib)) / smearing_sigma, ngauss);
demet += wk[ik] * smearing_sigma * occupy_smearing::w1gauss((ef - ekb(ik, ib)) / smearing_sigma, ngauss);
}
}

Expand All @@ -267,7 +270,7 @@ void Occupy::gweights(const int nks,
* @param is spin
* @param isk array to point out each k belong to which spin
*/
void Occupy::efermig(const ModuleBase::matrix& ekb,
void occupy_smearing::efermig(const ModuleBase::matrix& ekb,
const int nband,
const int nks,
const double& nelec,
Expand Down Expand Up @@ -320,8 +323,8 @@ void Occupy::efermig(const ModuleBase::matrix& ekb,
int changetime = 0;
while (true)
{
const double sumkup = Occupy::sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, eup, is, isk);
const double sumklw = Occupy::sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, elw, is, isk);
const double sumkup = occupy_smearing::sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, eup, is, isk);
const double sumklw = occupy_smearing::sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, elw, is, isk);

if (changetime > 1000)
{
Expand All @@ -335,7 +338,7 @@ void Occupy::efermig(const ModuleBase::matrix& ekb,
std::cout << " sumklw = " << sumklw << std::endl;
std::cout << " sumkup - nelec = " << sumkup - nelec << std::endl;
std::cout << " sumklw - nelec = " << sumklw - nelec << std::endl;
ModuleBase::WARNING_QUIT("Occupy::efermig", "ERROS in SMEARING");
ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "ERROS in SMEARING");
// no need to break; quit directly
}
else if ((sumkup - nelec) < -eps)
Expand All @@ -360,7 +363,7 @@ void Occupy::efermig(const ModuleBase::matrix& ekb,
// change ef value
//======================
ef = (eup + elw) / 2.0;
const double sumkmid = sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, ef, is, isk);
const double sumkmid = occupy_smearing::sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, ef, is, isk);

if (std::abs(sumkmid - nelec) < eps)
{
Expand Down Expand Up @@ -392,7 +395,7 @@ void Occupy::efermig(const ModuleBase::matrix& ekb,
* @param isk array to point out each k belong to which spin
* @return (double) the number of states
*/
double Occupy::sumkg(const ModuleBase::matrix& ekb,
double occupy_smearing::sumkg(const ModuleBase::matrix& ekb,
const int nband,
const int nks,
const std::vector<double>& wk,
Expand All @@ -416,7 +419,7 @@ double Occupy::sumkg(const ModuleBase::matrix& ekb,
//===========================
// call wgauss
//===========================
sum1 += Occupy::wgauss((e - ekb(ik, ib)) / smearing_sigma, ngauss);
sum1 += occupy_smearing::wgauss((e - ekb(ik, ib)) / smearing_sigma, ngauss);
}
sum2 += wk[ik] * sum1;
}
Expand All @@ -433,7 +436,7 @@ double Occupy::sumkg(const ModuleBase::matrix& ekb,
return sum2;
}

double Occupy::wgauss(const double& x, const int n)
double occupy_smearing::wgauss(const double& x, const int n)
{
// ModuleBase::TITLE("Occupy","wgauss");
//=====================================================================
Expand Down Expand Up @@ -515,7 +518,7 @@ double Occupy::wgauss(const double& x, const int n)
return wga;
} // end function wgauss

double Occupy::w1gauss(const double& x, const int n)
double occupy_smearing::w1gauss(const double& x, const int n)
{
//========================================================================
// w1gauss(x,n) = \int_{-\infty}^x y delta(y) dy
Expand Down
74 changes: 48 additions & 26 deletions source/source_estate/occupy.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,50 @@
#include "source_base/matrix.h"
#include "source_base/vector3.h"

/**
* @brief Smearing weight functions and the Fermi-level bisection built on them.
*
* Each of these depends only on its arguments -- none of them reads any of
* Occupy's state -- so they are free functions rather than private static
* members of that class.
*/
namespace occupy_smearing
{
/**
* @brief smeared occupation of a state, as a function of (ef - e) / sigma
* @param n smearing type: -99 Fermi-Dirac, -1 cold, 0 Gaussian,
* n > 0 n-th order Methfessel-Paxton
*/
double wgauss(const double& x, const int n);

/// @brief first moment of the smearing function, used for the -TS term
/// @param n smearing type, as in wgauss()
double w1gauss(const double& x, const int n);

/// @brief number of electrons below a trial Fermi level e
double sumkg(const ModuleBase::matrix& ekb,
const int nband,
const int nks,
const std::vector<double>& wk,
const double& smearing_sigma,
const int ngauss,
const double& e,
const int& is,
const std::vector<int>& isk);

/// @brief find the Fermi level by bisection on sumkg()
void efermig(const ModuleBase::matrix& ekb,
const int nbnd,
const int nks,
const double& nelec,
const std::vector<double>& wk,
const double& smearing_sigma,
const int ngauss,
double& ef,
const int& is,
const std::vector<int>& isk);
} // namespace occupy_smearing

class Occupy
{

Expand Down Expand Up @@ -35,6 +79,8 @@ class Occupy
// fixed occupations
static bool fixed_occupations;

/// @param nspin number of spin channels; 2 selects a single spin through
/// isk when is != -1, and 4 halves the spin degeneracy
static void iweights(const int nks,
const std::vector<double>& wk,
const int nband,
Expand All @@ -44,7 +90,8 @@ class Occupy
double& ef,
ModuleBase::matrix& wg,
const int& is,
const std::vector<int>& isk);
const std::vector<int>& isk,
const int nspin);

static void gweights(const int nks,
const std::vector<double>& wk,
Expand All @@ -66,31 +113,6 @@ class Occupy
static double wsweight(const ModuleBase::Vector3<double> &r, ModuleBase::Vector3<double> *rws,const int nrws);

private:
static void efermig(const ModuleBase::matrix& ekb,
const int nbnd,
const int nks,
const double& nelec,
const std::vector<double>& wk,
const double& smearing_sigma,
const int ngauss,
double& ef,
const int& is,
const std::vector<int>& isk);

static double sumkg(const ModuleBase::matrix& ekb,
const int nband,
const int nks,
const std::vector<double>& wk,
const double& smearing_sigma,
const int ngauss,
const double& e,
const int& is,
const std::vector<int>& isk);

static double wgauss(const double& x, const int n);

static double w1gauss(const double& x, const int n);

//============================
// Needed in tweights
//============================
Expand Down
41 changes: 17 additions & 24 deletions source/source_estate/test/elecstate_occupy_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,9 +1,7 @@
#include <string>
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#define private public
#include "source_io/module_parameter/parameter.h"
#undef private

/***************************************************************
* unit test of class Occupy
Expand All @@ -14,9 +12,7 @@
* - Occupy::Occupy()
* - Occupy::decision()
*/
#define private public
#include "source_estate/occupy.h"
#undef private
class OccupyTest : public ::testing::Test
{
protected:
Expand Down Expand Up @@ -182,21 +178,19 @@ TEST_F(OccupyTest, DecisionArbitrary)

TEST_F(OccupyTest, IweightsNOSPIN)
{
PARAM.input.nspin = 1;
double ef = 0.0;
ModuleBase::matrix wg(1, 1);
std::vector<double> wk(1, 2.0);
ModuleBase::matrix ekb(1, 1);
std::vector<int> isk(1);
ekb(0, 0) = 0.1;
occupy.iweights(1, wk, 1, 0, 2.0, ekb, ef, wg, 0, isk);
occupy.iweights(1, wk, 1, 0, 2.0, ekb, ef, wg, 0, isk, /*nspin=*/1);
EXPECT_DOUBLE_EQ(wg(0, 0), 2.0);
EXPECT_DOUBLE_EQ(ef, 0.1);
}

TEST_F(OccupyTest, IweightsSPIN)
{
PARAM.input.nspin = 2;
double ef_up = 0.0;
double ef_dw = 0.0;
ModuleBase::matrix wg(2, 1);
Expand All @@ -207,8 +201,8 @@ TEST_F(OccupyTest, IweightsSPIN)
isk[1] = 1;
ekb(0, 0) = 0.1;
ekb(1, 0) = 0.2;
occupy.iweights(2, wk, 1, 0, 1.0, ekb, ef_up, wg, 0, isk);
occupy.iweights(2, wk, 1, 0, 1.0, ekb, ef_dw, wg, 1, isk);
occupy.iweights(2, wk, 1, 0, 1.0, ekb, ef_up, wg, 0, isk, /*nspin=*/2);
occupy.iweights(2, wk, 1, 0, 1.0, ekb, ef_dw, wg, 1, isk, /*nspin=*/2);
EXPECT_DOUBLE_EQ(wg(0, 0), 1.0);
EXPECT_DOUBLE_EQ(wg(1, 0), 1.0);
EXPECT_DOUBLE_EQ(ef_up, 0.1);
Expand All @@ -217,7 +211,6 @@ TEST_F(OccupyTest, IweightsSPIN)

TEST_F(OccupyTest, IweightsWarning)
{
PARAM.input.nspin = 1;
double ef = 0.0;
ModuleBase::matrix wg(1, 1);
std::vector<double> wk(1, 2.0);
Expand All @@ -226,28 +219,28 @@ TEST_F(OccupyTest, IweightsWarning)
ekb(0, 0) = 0.1;

testing::internal::CaptureStdout();
EXPECT_EXIT(occupy.iweights(1, wk, 1, 0, 1.0, ekb, ef, wg, -1, isk);, ::testing::ExitedWithCode(1), "");
EXPECT_EXIT(occupy.iweights(1, wk, 1, 0, 1.0, ekb, ef, wg, -1, isk, /*nspin=*/1);, ::testing::ExitedWithCode(1), "");
output = testing::internal::GetCapturedStdout();
EXPECT_THAT(output, testing::HasSubstr("It is not a semiconductor or insulator. Please do not set 'smearing_method=fixed', and try other options."));
}

TEST_F(OccupyTest, Wgauss)
{
EXPECT_DOUBLE_EQ(occupy.wgauss(0.0, 0), 0.5);
EXPECT_DOUBLE_EQ(occupy.wgauss(0.0, -1), 0.4006259784506005);
EXPECT_DOUBLE_EQ(occupy.wgauss(0.0, -99), 0.5);
EXPECT_DOUBLE_EQ(occupy.wgauss(0.0, 1), 0.5);
EXPECT_DOUBLE_EQ(occupy.wgauss(0.0, 2), 0.5);
EXPECT_DOUBLE_EQ(occupy.wgauss(10, 0), 1.0);
EXPECT_DOUBLE_EQ(occupy_smearing::wgauss(0.0, 0), 0.5);
EXPECT_DOUBLE_EQ(occupy_smearing::wgauss(0.0, -1), 0.4006259784506005);
EXPECT_DOUBLE_EQ(occupy_smearing::wgauss(0.0, -99), 0.5);
EXPECT_DOUBLE_EQ(occupy_smearing::wgauss(0.0, 1), 0.5);
EXPECT_DOUBLE_EQ(occupy_smearing::wgauss(0.0, 2), 0.5);
EXPECT_DOUBLE_EQ(occupy_smearing::wgauss(10, 0), 1.0);
}

TEST_F(OccupyTest, W1gauss)
{
EXPECT_DOUBLE_EQ(occupy.w1gauss(0.0, 0), -0.28209479177387814);
EXPECT_DOUBLE_EQ(occupy.w1gauss(0.0, -1), -0.1710991401561083);
EXPECT_DOUBLE_EQ(occupy.w1gauss(0.0, -99), -0.69314718055994529);
EXPECT_DOUBLE_EQ(occupy.w1gauss(0.0, 1), -0.14104739588693907);
EXPECT_DOUBLE_EQ(occupy.w1gauss(0.0, 2), -0.10578554691520431);
EXPECT_DOUBLE_EQ(occupy_smearing::w1gauss(0.0, 0), -0.28209479177387814);
EXPECT_DOUBLE_EQ(occupy_smearing::w1gauss(0.0, -1), -0.1710991401561083);
EXPECT_DOUBLE_EQ(occupy_smearing::w1gauss(0.0, -99), -0.69314718055994529);
EXPECT_DOUBLE_EQ(occupy_smearing::w1gauss(0.0, 1), -0.14104739588693907);
EXPECT_DOUBLE_EQ(occupy_smearing::w1gauss(0.0, 2), -0.10578554691520431);
}

TEST_F(OccupyTest, Sumkg)
Expand All @@ -260,7 +253,7 @@ TEST_F(OccupyTest, Sumkg)
double e = 0.0;
int is = 0;
std::vector<int> isk = {0, 0};
EXPECT_DOUBLE_EQ(occupy.sumkg(ekb, 1, 1, wk, smearing_sigma, ngauss, e, is, isk), 1.0);
EXPECT_DOUBLE_EQ(occupy_smearing::sumkg(ekb, 1, 1, wk, smearing_sigma, ngauss, e, is, isk), 1.0);
}

TEST_F(OccupyTest, Efermig)
Expand All @@ -274,7 +267,7 @@ TEST_F(OccupyTest, Efermig)
int is = 0;
std::vector<int> isk = {0, 0};
double ef = 0.0;
occupy.efermig(ekb, 1, 1, 1.0, wk, smearing_sigma, ngauss, ef, is, isk);
occupy_smearing::efermig(ekb, 1, 1, 1.0, wk, smearing_sigma, ngauss, ef, is, isk);
EXPECT_NEAR(ef, -0.5, 1e-13);
}

Expand Down
Loading