temporarily remove modified stratmann

Author: jinzx10

source/module_base/grid/partition.cpp

@@ -10,8 +10,6 @@ namespace Grid {
 namespace Partition {
 
 const double stratmann_a = 0.64;
-const double stratmann_mod_b = 0.8;
-
 
 double w_becke(
     int nR0,
@@ -135,76 +133,5 @@ double s_stratmann(double mu) {
 }
 
 
-double w_stratmann_mod(
-    int nR0,
-    const double* drR,
-    const double* dRR,
-    const double* drR_thr,
-    const double* Rcut,
-    int nR,
-    int* iR,
-    int c
-) {
-    int I = iR[c], J = 0;
-
-    // If r falls within the exclusive zone of a center, return immediately.
-    for (int j = 0; j < nR; ++j) {
-        J = iR[j];
-        if (drR[J] <= drR_thr[J]) {
-            return static_cast<double>(I == J);
-        }
-    }
-
-    // Even if the grid point does not fall within the exclusive zone of any
-    // center, the normalized weight could still be 0 or 1, and this can be
-    // figured out by examining the unnormalized weight alone.
-
-    // Move the center to the first position for convenience. Swap back later.
-    std::swap(iR[0], iR[c]);
-
-    double P[nR];
-    for (int j = 1; j < nR; ++j) {
-        J = iR[j];
-        double mu = (drR[I] - drR[J]) / dRR[I*nR0 + J];
-        P[j] = s_stratmann(mu);
-    }
-    P[0] = std::accumulate(P + 1, P + nR, 1.0, std::multiplies<double>());
-
-    if (P[0] == 0.0 || P[0] == 1.0) {
-        std::swap(iR[0], iR[c]);
-        return P[0];
-    }
-
-    // If it passes all the screening, all unnormalized weights have to be
-    // calculated in order to get the normalized weight.
-
-    std::for_each(P + 1, P + nR, [](double& s) { s = 1.0 - s; });
-    for (int i = 1; i < nR; ++i) {
-        I = iR[i];
-        for (int j = i + 1; j < nR; ++j) {
-            J = iR[j];
-            double mu = (drR[I] - drR[J]) / dRR[I*nR0 + J];
-            double s = s_stratmann(mu);
-            P[i] *= s;
-            P[j] *= (1.0 - s); // s(-mu) = 1 - s(mu)
-        }
-    }
-
-    std::swap(iR[0], iR[c]);
-    return P[0] / std::accumulate(P, P + nR, 0.0);
-}
-
-
-double s_stratmann_mod(double mu, double y) {
-    using ModuleBase::PI;
-    bool core = y <= stratmann_mod_b;
-    bool edge = !core && y < 1.0;
-    double u = core + edge * 0.5 * (1.0 +
-            std::cos(PI * (y - stratmann_mod_b) / (1.0 - stratmann_mod_b)));
-    return 1.0 + u * (s_stratmann(mu) - 1.0);
-}
-
-
-
 } // end of namespace Partition
 } // end of namespace Grid

source/module_base/grid/partition.h

@@ -7,11 +7,9 @@ namespace Partition {
 enum class Type {
     Becke,
     Stratmann,
-    StratmannMod
 };
 
 extern const double stratmann_a;
-extern const double stratmann_mod_b;
 
 /**
  * @brief Becke's partition weight.
@@ -80,32 +78,6 @@ double w_stratmann(
 // Stratmann's piecewise cell function
 double s_stratmann(double mu);
 
-
-/**
- * @brief Becke's partition weight with a modified Stratmann's scheme
- * by Knuth et al.
- *
- * Reference:
- * Knuth, F., Carbogno, C., Atalla, V., Blum, V., & Scheffler, M. (2015).
- * All-electron formalism for total energy strain derivatives and stress
- * tensor components for numeric atom-centered orbitals.
- * Computer Physics Communications, 190, 33-50.
- *
- */
-double w_stratmann_mod(
-    int nR0,
-    const double* drR,
-    const double* dRR,
-    const double* drR_thr,
-    const double* Rcut,
-    int nR,
-    int* iR,
-    int c
-);
-
-double s_stratmann_mod(double mu, double y); // y = || r - R[J] || / Rcut[J]
-
-
 } // end of namespace Partition
 } // end of namespace Grid
 

source/module_base/grid/test/test_partition.cpp

@@ -8,6 +8,8 @@
 #include <array>
 #include <numeric>
 #include <chrono>
+#include <random>
+#include <algorithm>
 
 #ifdef __MPI
 #include <mpi.h>
@@ -154,7 +156,7 @@ PartitionTest::PartitionTest() {
     std::vector<double> r_rad, w_rad;
     int nrad = 60;
     int Rcut = 7.0;
-    int mult = 4;
+    int mult = 2;
     Grid::Radial::baker(nrad, Rcut, r_rad, w_rad, mult);
 
     // complete grid & weight for one-center integration
@@ -189,6 +191,10 @@ TEST_F(PartitionTest, Becke) {
         std::vector<int> iR(nR);
         std::iota(iR.begin(), iR.end(), 0);
 
+        std::random_device rd;
+        std::mt19937 g(rd());
+        std::shuffle(iR.begin(), iR.end(), g);
+
         for (size_t I = 0; I < nR; ++I) { // for each center
             for (size_t i = 0; i < w.size(); i++) {
                 Vec3 ri = Vec3{r[3*i], r[3*i+1], r[3*i+2]} + param.R[I];
@@ -271,6 +277,7 @@ TEST_F(PartitionTest, Stratmann) {
     printf("time elapsed = %8.3e seconds\n", dur.count());
 }
 
+
 int main(int argc, char** argv)
 {
 #ifdef __MPI