occ bug fix...EJB

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

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 

Nwpw/pspw/lib/kinetic/Kinetic.cpp

@@ -139,7 +139,7 @@ double Kinetic_Operator::ke_ave(double *psi, double *occ)
    for (auto ms=0; ms<mypneb->ispin; ++ms)
    for (auto q=0; q<mypneb->neq[ms]; ++q) 
    {
-     double wght = occ[mypneb->msntoindex(ms,q)];
+     double wght = occ ? occ[mypneb->msntoindex(ms,q)] : 1.0;
      for (k=0; k<ksize1; ++k) 
      {
         ave += tg[k] * (psi[k1] * psi[k1] + psi[k2] * psi[k2])*wght;

Nwpw/pspw/lib/psi/psi.cpp

@@ -326,17 +326,16 @@ void psi_read0(Pneb *mypneb, int *version, int nfft[], double unita[],
    {
       if (myparall->is_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 */

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