Skip to content

Commit 7cb7664

Browse files
jinzx10jinzx10
authored
Becke's partition for multi-center grid integration (#5292)
* becke partition initial commit * stratmann * stratmann revised * temporarily remove modified stratmann * use std::vector instead of builtin array
1 parent 1a11dac commit 7cb7664

File tree

6 files changed

+534
-10
lines changed

6 files changed

+534
-10
lines changed

‎source/module_base/grid/partition.cpp‎

Lines changed: 142 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
#include "module_base/grid/partition.h"
2+
#include "module_base/constants.h"
3+
4+
#include <cmath>
5+
#include <functional>
6+
#include <numeric>
7+
#include <algorithm>
8+
#include <vector>
9+
#include <cassert>
10+
11+
namespace Grid {
12+
namespace Partition {
13+
14+
const double stratmann_a = 0.64;
15+
16+
double w_becke(
17+
int nR0,
18+
const double* drR,
19+
const double* dRR,
20+
int nR,
21+
const int* iR,
22+
int c
23+
) {
24+
assert(nR > 0 && nR0 >= nR);
25+
std::vector<double> P(nR, 1.0);
26+
for (int i = 0; i < nR; ++i) {
27+
int I = iR[i];
28+
for (int j = i + 1; j < nR; ++j) {
29+
int J = iR[j];
30+
double mu = (drR[I] - drR[J]) / dRR[I*nR0 + J];
31+
double s = s_becke(mu);
32+
P[I] *= s;
33+
P[J] *= (1.0 - s); // s(-mu) = 1 - s(mu)
34+
}
35+
}
36+
return P[c] / std::accumulate(P.begin(), P.end(), 0.0);
37+
}
38+
39+
40+
double s_becke(double mu) {
41+
/*
42+
* Becke's iterated polynomials (3rd order)
43+
*
44+
* s(mu) = 0.5 * (1 - p(p(p(mu))))
45+
*
46+
* p(x) = 0.5 * x * (3 - x^2)
47+
*
48+
*/
49+
double p = 0.5 * mu * (3.0 - mu*mu);
50+
p = 0.5 * p * (3.0 - p*p);
51+
p = 0.5 * p * (3.0 - p*p);
52+
return 0.5 * (1.0 - p);
53+
}
54+
55+
56+
double w_stratmann(
57+
int nR0,
58+
const double* drR,
59+
const double* dRR,
60+
const double* drR_thr,
61+
int nR,
62+
int* iR,
63+
int c
64+
) {
65+
assert(nR > 0 && nR0 >= nR);
66+
int I = iR[c], J = 0;
67+
68+
// If r falls within the exclusive zone of a center, return immediately.
69+
for (int j = 0; j < nR; ++j) {
70+
J = iR[j];
71+
if (drR[J] <= drR_thr[J]) {
72+
return static_cast<double>(I == J);
73+
}
74+
}
75+
76+
// Even if the grid point does not fall within the exclusive zone of any
77+
// center, the normalized weight could still be 0 or 1, and this can be
78+
// figured out by examining the unnormalized weight alone.
79+
80+
// Swap the grid center to the first position in iteration for convenience.
81+
// Restore the original order before return.
82+
std::swap(iR[0], iR[c]);
83+
84+
std::vector<double> P(nR);
85+
for (int j = 1; j < nR; ++j) {
86+
J = iR[j];
87+
double mu = (drR[I] - drR[J]) / dRR[I*nR0 + J];
88+
P[j] = s_stratmann(mu);
89+
}
90+
P[0] = std::accumulate(P.begin() + 1, P.end(), 1.0,
91+
std::multiplies<double>());
92+
93+
if (P[0] == 0.0 || P[0] == 1.0) {
94+
std::swap(iR[0], iR[c]); // restore the original order
95+
return P[0];
96+
}
97+
98+
// If it passes all the screening above, all unnormalized weights
99+
// have to be calculated in order to get the normalized weight.
100+
101+
std::for_each(P.begin() + 1, P.end(), [](double& s) { s = 1.0 - s; });
102+
for (int i = 1; i < nR; ++i) {
103+
I = iR[i];
104+
for (int j = i + 1; j < nR; ++j) {
105+
J = iR[j];
106+
double mu = (drR[I] - drR[J]) / dRR[I*nR0 + J];
107+
double s = s_stratmann(mu);
108+
P[i] *= s;
109+
P[j] *= (1.0 - s); // s(-mu) = 1 - s(mu)
110+
}
111+
}
112+
113+
std::swap(iR[0], iR[c]); // restore the original order
114+
return P[0] / std::accumulate(P.begin(), P.end(), 0.0);
115+
}
116+
117+
118+
double s_stratmann(double mu) {
119+
/*
120+
* Stratmann's piecewise cell function
121+
*
122+
* s(mu) = 0.5 * (1 - g(mu/a))
123+
*
124+
* / -1 x <= -1
125+
* |
126+
* g(x) = | (35x - 35x^3 + 21x^5 - 5x^7) / 16 |x| < 1
127+
* |
128+
* \ +1 x >= +1
129+
*
130+
*/
131+
double x = mu / stratmann_a;
132+
double x2 = x * x;
133+
double h = 0.0625 * x * (35 + x2 * (-35 + x2 * (21 - 5 * x2)));
134+
135+
bool mid = std::abs(x) < 1;
136+
double g = !mid * (1 - 2 * std::signbit(x)) + mid * h;
137+
return 0.5 * (1.0 - g);
138+
}
139+
140+
141+
} // end of namespace Partition
142+
} // end of namespace Grid

‎source/module_base/grid/partition.h‎

Lines changed: 84 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,84 @@
1+
#ifndef GRID_PARTITION_H
2+
#define GRID_PARTITION_H
3+
4+
namespace Grid {
5+
namespace Partition {
6+
7+
enum class Type {
8+
Becke,
9+
Stratmann,
10+
};
11+
12+
extern const double stratmann_a;
13+
14+
/**
15+
* @brief Becke's partition weight.
16+
*
17+
* This function computes the normalized Becke's partition weight
18+
* for a grid point associated with a selected set of centers, given
19+
* the grid point's distance to centers and inter-center distance.
20+
*
21+
* @param nR0 Total number of centers given by drR & dRR.
22+
* @param drR Distance between the grid point and centers.
23+
* @param dRR Distance between centers. dRR[I*nR0 + J] is the
24+
* distance between center I and J.
25+
* @param nR Number of centers involved in the weight calculation.
26+
* nR <= nR0. Length of iR.
27+
* @param iR Indices of centers involved.
28+
* Each element is a distinctive index in [0, nR0).
29+
* @param c iR[c] is the index of the center whom this grid point
30+
* belongs to.
31+
*
32+
* Reference:
33+
* Becke, A. D. (1988).
34+
* A multicenter numerical integration scheme for polyatomic molecules.
35+
* The Journal of chemical physics, 88(4), 2547-2553.
36+
*
37+
*/
38+
double w_becke(
39+
int nR0,
40+
const double* drR,
41+
const double* dRR,
42+
int nR,
43+
const int* iR,
44+
int c
45+
);
46+
47+
// Becke's cell function (iterated polynomial)
48+
double s_becke(double mu);
49+
50+
51+
/**
52+
* @brief Becke's partition weight with Stratmann's scheme.
53+
*
54+
* This function is similar to w_becke, but the cell function adopts
55+
* the one proposed by Stratmann et al, which enables some screening.
56+
*
57+
* @see w_becke
58+
*
59+
* @param drR_thr Radius of exclusive zone of each center.
60+
*
61+
* Reference:
62+
* Stratmann, R. E., Scuseria, G. E., & Frisch, M. J. (1996).
63+
* Achieving linear scaling in exchange-correlation density functional
64+
* quadratures.
65+
* Chemical physics letters, 257(3-4), 213-223.
66+
*
67+
*/
68+
double w_stratmann(
69+
int nR0,
70+
const double* drR,
71+
const double* dRR,
72+
const double* drR_thr,
73+
int nR,
74+
int* iR,
75+
int c
76+
);
77+
78+
// Stratmann's piecewise cell function
79+
double s_stratmann(double mu);
80+
81+
} // end of namespace Partition
82+
} // end of namespace Grid
83+
84+
#endif

‎source/module_base/grid/radial.cpp‎

Lines changed: 3 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -2,14 +2,6 @@
22

33
#include <cmath>
44

5-
namespace {
6-
7-
const double pi = std::acos(-1.0);
8-
const double inv_ln2 = 1.0 / std::log(2.0);
9-
10-
} // end of anonymous namespace
11-
12-
135
namespace Grid {
146
namespace Radial {
157

@@ -43,6 +35,9 @@ void murray(int n, double R, double* r, double* w) {
4335

4436

4537
void treutler_m4(int n, double R, double* r, double* w, double alpha) {
38+
const double pi = std::acos(-1.0);
39+
const double inv_ln2 = 1.0 / std::log(2.0);
40+
4641
for (int i = 1; i <= n; ++i) {
4742
double x = std::cos(i * pi / (n + 1));
4843
double beta = std::sqrt((1.0 + x) / (1.0 - x));

‎source/module_base/grid/test/CMakeLists.txt‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,3 +12,11 @@ AddTest(
1212
SOURCES test_radial.cpp
1313
../radial.cpp
1414
)
15+
16+
AddTest(
17+
TARGET test_partition
18+
SOURCES test_partition.cpp
19+
../partition.cpp
20+
../radial.cpp
21+
../delley.cpp
22+
)

‎source/module_base/grid/test/test_delley.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -101,8 +101,8 @@ TEST_F(DelleyTest, Accuracy) {
101101
ModuleBase::Ylm::sph_harm(func_lmax,
102102
grid[3*i], grid[3*i+1], grid[3*i+2], ylm_real);
103103
double tmp = 0.0;
104-
for (size_t i = 0; i < coef.size(); ++i) {
105-
tmp += coef[i] * ylm_real[i];
104+
for (size_t j = 0; j < coef.size(); ++j) {
105+
tmp += coef[j] * ylm_real[j];
106106
}
107107
val += weight[i] * tmp * tmp;
108108
}

0 commit comments

Comments
 (0)