diff --git a/Nwpw/nwpwlib/utilities/util.cpp b/Nwpw/nwpwlib/utilities/util.cpp index 18d640e4..7bb5099f 100644 --- a/Nwpw/nwpwlib/utilities/util.cpp +++ b/Nwpw/nwpwlib/utilities/util.cpp @@ -1028,9 +1028,22 @@ double util_occupation_distribution(const int smeartype, const double e) } else if (smeartype == 2) { // Gaussian f = 0.5 * std::erfc(e); } else if (smeartype == 3) { // Hermite smearing - double exp_term = std::exp(-e * e); - double hermite_correction = (2.0 * e * e - 1.0) * exp_term; - f = exp_term + hermite_correction / sqrt_pi; + double exp_term = std::exp(-e * e); // Gaussian term + double hermite_correction = (2.0 * e * e - 1.0) * exp_term; // Hermite correction term + + // Compute the original Hermite-Gaussian function + double f_original = exp_term + hermite_correction / sqrt_pi; + + // Find fmax by evaluating the original function at e=0 + double fmax = exp_term + (2.0 * 0 * 0 - 1.0) * exp_term / sqrt_pi; // f_original at e=0 + + if (e <= 0.0) { + // Set value to 1.0 for e <= e_max + f = 1.0; + } else { + // Scale the Hermite-Gaussian function for e > e_max + f = f_original / fmax; + } } else if (smeartype == 4) { // Marzari-Vanderbilt double factor = std::sqrt(0.125 / std::atan(1.0)); // atan(1.0) = pi/4 f = std::exp(-(e + sqrt_half) * (e + sqrt_half)) * factor + 0.5 * std::erfc(e + sqrt_half); @@ -1130,8 +1143,14 @@ double util_smearcorrection(const int smeartype, const double smearkT, const dou smearcorrection -= smearkT * std::exp(-x * x) / (4.0 * sqrt_pi); } else if (smeartype == 3) { // Hermite smearing double exp_term = std::exp(-x * x); - double hermite_poly = 2.0 * x * x - 1.0; - smearcorrection += smearkT * exp_term * hermite_poly / sqrt_pi; + double hermite_correction = (2.0 * x * x - 1.0) * exp_term; + double f_original = exp_term + hermite_correction / sqrt_pi; + double fmax = std::exp(-0.0) + (2.0 * 0.0 * 0.0 - 1.0) * std::exp(-0.0) / sqrt_pi; + + if (x >0.0) { // Smearing correction - No correction needed for x <= 0, as f(x) = 1.0 + // Scale the correction for x > 0 + smearcorrection -= smearkT * f_original / fmax; + } } else if (smeartype == 4) { // Marzari-Vanderbilt correction smearcorrection -= smearkT * exp(-(x +sqrt_half)*(x+sqrt_half))*(1.0+sqrt_two*x) / (2.0*sqrt_pi); } else if (smeartype == 5) { // Methfessel-Paxton diff --git a/Nwpw/pspw/cpsd/cpsd.cpp b/Nwpw/pspw/cpsd/cpsd.cpp index dd460dea..8599a700 100644 --- a/Nwpw/pspw/cpsd/cpsd.cpp +++ b/Nwpw/pspw/cpsd/cpsd.cpp @@ -436,7 +436,7 @@ int cpsd(MPI_Comm comm_world0, std::string &rtdbstring) if (smeartype==5) std::cout << "Methfessel-Paxton" << std::endl; if (smeartype==6) std::cout << "Cold smearing" << std::endl; if (smeartype==7) std::cout << "Lorentzian" << std::endl; - if (smeartype==8) std::cout << "No correction" << std::endl; + if (smeartype==8) std::cout << "step" << std::endl; if (smeartype>=0) { std::cout << " smearing parameter = " << Ffmt(9,3) << smearkT diff --git a/Nwpw/pspw/lib/kinetic/Kinetic.cpp b/Nwpw/pspw/lib/kinetic/Kinetic.cpp index df58d04e..7485d14b 100644 --- a/Nwpw/pspw/lib/kinetic/Kinetic.cpp +++ b/Nwpw/pspw/lib/kinetic/Kinetic.cpp @@ -139,7 +139,7 @@ double Kinetic_Operator::ke_ave(double *psi, double *occ) for (auto ms=0; msispin; ++ms) for (auto q=0; qneq[ms]; ++q) { - double wght = occ[mypneb->msntoindex(ms,q)]; + double wght = occ ? occ[mypneb->msntoindex(ms,q)] : 1.0; for (k=0; kis_master()) { - dread(4,occ, ne[0]+ne[1]); - - if (reverse) - { - std::reverse(occ, occ + ne[0]); // Reverse the first section if descending - - // Check and reverse the second part - if (*ispin>1) - std::reverse(occ + ne[0], occ + ne[0] + ne[1]); // Reverse the second section if descending and *ispin > 1 - } - + dread(4,occ, ne[0]+ne[1]); + + if (reverse) + { + std::reverse(occ, occ + ne[0]); // Reverse the first section if descending + + // Check and reverse the second part + if (*ispin>1) + std::reverse(occ + ne[0], occ + ne[0] + ne[1]); // Reverse the second section if descending and *ispin > 1 + } } myparall->Brdcst_Values(0, 0, ne[0]+ne[1], occ); } @@ -415,12 +414,13 @@ bool psi_read(Pneb *mypneb, char *filename, bool wvfnc_initialize, double *psi2, psi_read0(mypneb, &version, nfft, unita, &ispin, ne, psi2, occupation, occ2, filename,false); - if (isDescending(occ2, ne[0])) - { - if (myparall->base_stdio_print) - coutput << " - reversing order of psi and occupation" << std::endl; - psi_read0(mypneb, &version, nfft, unita, &ispin, ne, psi2, occupation, occ2, filename,true); - } + if (occ2) + if (isDescending(occ2, ne[0])) + { + if (myparall->base_stdio_print) + coutput << " - reversing order of psi and occupation" << std::endl; + psi_read0(mypneb, &version, nfft, unita, &ispin, ne, psi2, occupation, occ2, filename,true); + } } /* generate new psi */ diff --git a/Nwpw/pspw/minimizer/pspw_minimizer.cpp b/Nwpw/pspw/minimizer/pspw_minimizer.cpp index 7fa5d36b..9066817c 100644 --- a/Nwpw/pspw/minimizer/pspw_minimizer.cpp +++ b/Nwpw/pspw/minimizer/pspw_minimizer.cpp @@ -422,7 +422,7 @@ int pspw_minimizer(MPI_Comm comm_world0, std::string &rtdbstring, std::ostream & if (mymolecule.smeartype==5) coutput << "Methfessel-Paxton" << std::endl; if (mymolecule.smeartype==6) coutput << "Cold smearing" << std::endl; if (mymolecule.smeartype==7) coutput << "Lorentzian" << std::endl; - if (mymolecule.smeartype==8) coutput << "No correction" << std::endl; + if (mymolecule.smeartype==8) coutput << "step" << std::endl; if (mymolecule.smeartype>=0) { coutput << " smearing parameter = " << Ffmt(9,3) << mymolecule.smearkT