From e9975c834ce28e79c42e3ce08998808f578ce756 Mon Sep 17 00:00:00 2001 From: Qianruipku Date: Sun, 27 Aug 2023 17:27:46 +0800 Subject: [PATCH 1/4] Fix: bug of NAN in SDFT --- .../hamilt_stodft/sto_func.cpp | 39 ++++++++++++------- 1 file changed, 24 insertions(+), 15 deletions(-) diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_func.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_func.cpp index f1c200a6ca..d9a0a33429 100644 --- a/source/module_hamilt_pw/hamilt_stodft/sto_func.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/sto_func.cpp @@ -68,32 +68,38 @@ REAL Sto_Func:: nxfd(REAL rawe) template REAL Sto_Func:: fdlnfd(REAL e) { - REAL e_mu = (e - mu) / this->tem ; - if(e_mu > 36) + REAL e_mu = (e - mu) / this->tem; + if (e_mu > 36) return 0; - else if(e_mu < -36) + else if (e_mu < -36) return 0; else { REAL f = 1 / (1 + exp(e_mu)); - return (f * log(f) + (1.0-f) * log(1.0-f)); + if (f == 0 || f == 1) + return 0; + else + return (f * log(f) + (1.0 - f) * log(1.0 - f)); } } template REAL Sto_Func:: nfdlnfd(REAL rawe) { - REAL Ebar = (Emin + Emax)/2; - REAL DeltaE = (Emax - Emin)/2; - REAL ne_mu = (rawe * DeltaE + Ebar - mu) / this->tem ; - if(ne_mu > 36) + REAL Ebar = (Emin + Emax) / 2; + REAL DeltaE = (Emax - Emin) / 2; + REAL ne_mu = (rawe * DeltaE + Ebar - mu) / this->tem; + if (ne_mu > 36) return 0; - else if(ne_mu < -36) + else if (ne_mu < -36) return 0; else { REAL f = 1 / (1 + exp(ne_mu)); - return f * log(f) + (1-f) * log(1-f); + if (f == 0 || f == 1) + return 0; + else + return f * log(f) + (1 - f) * log(1 - f); } } @@ -103,14 +109,17 @@ REAL Sto_Func:: n_root_fdlnfd(REAL rawe) REAL Ebar = (Emin + Emax)/2; REAL DeltaE = (Emax - Emin)/2; REAL ne_mu = (rawe * DeltaE + Ebar - mu) / this->tem ; - if(ne_mu > 72) + if(ne_mu > 36) return 0; - else if(ne_mu < -72) + else if(ne_mu < -36) return 0; else { REAL f = 1 / (1 + exp(ne_mu)); - return sqrt(-f * log(f) - (1-f) * log(1-f)); + if (f == 0 || f == 1) + return 0; + else + return sqrt(-f * log(f) - (1-f) * log(1-f)); } } @@ -170,7 +179,7 @@ REAL Sto_Func::ngauss(REAL rawe) REAL DeltaE = (Emax - Emin)/2; REAL e = rawe * DeltaE + Ebar; REAL a = pow((targ_e-e),2)/2.0/pow(sigma,2); - if(a > 32) + if(a > 72) return 0; else return exp(-a) /sqrt(TWOPI) / sigma ; @@ -183,7 +192,7 @@ REAL Sto_Func::nroot_gauss(REAL rawe) REAL DeltaE = (Emax - Emin)/2; REAL e = rawe * DeltaE + Ebar; REAL a = pow((targ_e-e),2)/4.0/pow(sigma,2); - if(a > 32) + if(a > 72) return 0; else return exp(-a) /sqrt(sqrt(TWOPI) * sigma) ; From 57bca757f2cdf1cd11ba5d3692681375847e97a3 Mon Sep 17 00:00:00 2001 From: Qianruipku Date: Thu, 19 Oct 2023 12:59:25 +0800 Subject: [PATCH 2/4] fix bug when fixed smearing and nspin=2 --- source/module_elecstate/occupy.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/source/module_elecstate/occupy.cpp b/source/module_elecstate/occupy.cpp index 5fa40ac1fb..afd9a7d9a6 100644 --- a/source/module_elecstate/occupy.cpp +++ b/source/module_elecstate/occupy.cpp @@ -131,8 +131,12 @@ void Occupy::iweights(const int nks, const std::vector& isk) { assert(is < 2); + double degspin = 2.0; + if (GlobalV::NSPIN == 4) + degspin = 1.0; + if (is != -1) + degspin = 1.0; - double degspin = (GlobalV::NSPIN == 1) ? 2.0 : 1.0; double ib_mind = nelec / degspin; int ib_min = std::ceil(ib_mind); if (ib_min != int(ib_mind)) From 7274099beb67c3e1d8fb8c578706fd469bcf1e66 Mon Sep 17 00:00:00 2001 From: Qianruipku Date: Thu, 19 Oct 2023 13:23:22 +0800 Subject: [PATCH 3/4] update input --- docs/advanced/input_files/input-main.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index c2446eb52d..b8721155ab 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -881,7 +881,7 @@ calculations. - **Type**: String - **Description**: It indicates which occupation and smearing method is used in the calculation. - - **fixed**: fixed occupations. + - **fixed**: fixed occupations (available for non-coductors only) - **gauss** or **gaussian**: Gaussian smearing method. - **mp**: methfessel-paxton smearing method; recommended for metals. - **fd**: Fermi-Dirac smearing method: $f=1/\{1+\exp[(E-\mu)/kT]\}$ and smearing_sigma below is the temperature $T$ (in Ry). From dca582a0d1c42307d2872cd22a5bb16c6f3aea92 Mon Sep 17 00:00:00 2001 From: Qianruipku Date: Mon, 23 Oct 2023 16:21:07 +0800 Subject: [PATCH 4/4] fix UT of elecstate_occupy --- source/module_elecstate/test/elecstate_occupy_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/module_elecstate/test/elecstate_occupy_test.cpp b/source/module_elecstate/test/elecstate_occupy_test.cpp index 9b58996ddb..e95470cd98 100644 --- a/source/module_elecstate/test/elecstate_occupy_test.cpp +++ b/source/module_elecstate/test/elecstate_occupy_test.cpp @@ -225,7 +225,7 @@ TEST_F(OccupyTest, IweightsWarning) ekb(0, 0) = 0.1; testing::internal::CaptureStdout(); - EXPECT_EXIT(occupy.iweights(1, wk, 1, 1.0, ekb, ef, wg, 0, isk);, ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(occupy.iweights(1, wk, 1, 1.0, ekb, ef, wg, -1, isk);, ::testing::ExitedWithCode(0), ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("It is not a semiconductor or insulator. Change 'smearing_method'.")); }