Skip to content

Commit 3f8d425

Browse files
jinzx10jinzx10
committed
stratmann
1 parent 8d4c15e commit 3f8d425

File tree

4 files changed

+172
-37
lines changed

4 files changed

+172
-37
lines changed

‎source/module_base/grid/partition.cpp‎

Lines changed: 67 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,9 @@
88
namespace Grid {
99
namespace Partition {
1010

11+
const double stratmann_a = 0.64;
12+
const double stratmann_mod_b = 0.8;
13+
1114
double w_becke(
1215
int nR0,
1316
double* drR,
@@ -46,7 +49,62 @@ double s_becke(double mu) {
4649
}
4750

4851

49-
double s_stratmann(double mu, double a) {
52+
double w_stratmann(
53+
int nR0,
54+
double* drR,
55+
double* dRR,
56+
double* drR_thr,
57+
int nR,
58+
int* iR,
59+
int c
60+
) {
61+
62+
// if r falls within the exclusive zone of the center
63+
// whom this grid point belongs to
64+
int I = iR[c];
65+
if (drR[I] <= drR_thr[I]) {
66+
return 1.0;
67+
}
68+
69+
// if r falls within the exclusive zone of a center
70+
// other than the one whom this grid point belongs to
71+
for (int J = 0; J < nR; ++J) {
72+
// no need to exclude J == c because it was checked before
73+
if (drR[iR[J]] <= drR_thr[iR[J]]) {
74+
return 0.0;
75+
}
76+
}
77+
78+
double Pc = 1.0;
79+
for (int i = 0; i < c; ++i) {
80+
int J = iR[i];
81+
double mu = (drR[I] - drR[J]) / dRR[I*nR0 + J];
82+
Pc *= s_stratmann(mu);
83+
}
84+
for (int i = c + 1; i < nR; ++i) {
85+
int J = iR[i];
86+
double mu = (drR[I] - drR[J]) / dRR[I*nR0 + J];
87+
Pc *= s_stratmann(mu);
88+
}
89+
if (Pc == 0.0 || Pc == 1.0) {
90+
return Pc;
91+
}
92+
93+
std::vector<double> P(nR, 1.0);
94+
for (int i = 0; i < nR; ++i) {
95+
int I = iR[i];
96+
for (int j = i + 1; j < nR; ++j) {
97+
int J = iR[j];
98+
double mu = (drR[I] - drR[J]) / dRR[I*nR0 + J];
99+
double s = s_stratmann(mu);
100+
P[I] *= s;
101+
P[J] *= (1.0 - s); // s(-mu) = 1 - s(mu)
102+
}
103+
}
104+
return P[c] / std::accumulate(P.begin(), P.end(), 0.0);
105+
}
106+
107+
double s_stratmann(double mu) {
50108
/*
51109
* Stratmann's piecewise cell function
52110
*
@@ -59,28 +117,29 @@ double s_stratmann(double mu, double a) {
59117
* \ +1 x >= +1
60118
*
61119
*/
62-
double x = mu / a;
120+
double x = mu / stratmann_a;
63121
double x2 = x * x;
64-
double h = 0.625 * x * (35 + x2 * (-35 + x2 * (21 - 5 * x2)));
122+
double h = 0.0625 * x * (35 + x2 * (-35 + x2 * (21 - 5 * x2)));
65123

66124
bool mid = std::abs(x) < 1;
67125
double g = !mid * (1 - 2 * std::signbit(x)) + mid * h;
68126
return 0.5 * (1 - g);
69127
}
70128

71129

72-
double u_stratmann_mod(double y, double b) {
130+
double u_stratmann_mod(double y) {
73131
using ModuleBase::PI;
74-
bool core = y <= b;
132+
bool core = y <= stratmann_mod_b;
75133
bool edge = !core && y < 1.0;
76-
return core + edge * 0.5 * (std::cos(PI * (y - b) / (1.0 - b)) + 1.0);
134+
return core + edge * 0.5 * (1.0 +
135+
std::cos(PI * (y - stratmann_mod_b) / (1.0 - stratmann_mod_b)));
77136
}
78137

79138

80-
double s_stratmann_mod(double mu, double y, double a, double b) {
139+
double s_stratmann_mod(double mu, double y) {
81140
// Modified Stratmann's cell function by Knuth et al.
82141
// y = |r-R(J)| / Rcut(J)
83-
return 1.0 + u_stratmann_mod(y, b) * (s_stratmann(mu, a) - 1.0);
142+
return 1.0 + u_stratmann_mod(y) * (s_stratmann(mu) - 1.0);
84143
}
85144

86145

‎source/module_base/grid/partition.h‎

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,8 @@ enum class Type {
1010
StratmannMod
1111
};
1212

13+
extern const double stratmann_a;
14+
extern const double stratmann_mod_b;
1315

1416
/**
1517
* @brief Becke's partition weight.
@@ -56,7 +58,7 @@ double s_becke(double mu);
5658
*
5759
* @see w_becke
5860
*
59-
* @param dRcRnn Distance between the grid center and its nearest neighbor.
61+
* @param drR_thr Radius of exclusive zone of each center.
6062
*
6163
* Reference:
6264
* Stratmann, R. E., Scuseria, G. E., & Frisch, M. J. (1996).
@@ -69,15 +71,14 @@ double w_stratmann(
6971
int nR0,
7072
double* drR,
7173
double* dRR,
74+
double* drR_thr,
7275
int nR,
7376
int* iR,
74-
int c,
75-
double dRcRnn,
76-
double a = 0.64
77+
int c
7778
);
7879

7980
// Stratmann's piecewise cell function
80-
double s_stratmann(double mu, double a);
81+
double s_stratmann(double mu);
8182

8283

8384
/**
@@ -86,8 +87,8 @@ double s_stratmann(double mu, double a);
8687
* tensor components for numeric atom-centered orbitals.
8788
* Computer Physics Communications, 190, 33-50.
8889
*/
89-
double s_stratmann_mod(double mu, double y, double a = 0.64, double b = 0.8);
90-
double u_stratmann_mod(double y, double b); // used by s_stratmann_mod
90+
double s_stratmann_mod(double mu, double y);
91+
double u_stratmann_mod(double y); // used by s_stratmann_mod
9192

9293

9394
} // end of namespace Partition

‎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_partition.cpp‎

Lines changed: 89 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -98,8 +98,32 @@ std::vector<Param> test_params = {
9898
},
9999
};
100100

101+
std::vector<double> dist_R_R(const std::vector<Vec3>& R) {
102+
// tabulate dRR[I,J] = || R[I] - R[J] ||
103+
size_t nR = R.size();
104+
std::vector<double> dRR(nR*nR, 0.0);
105+
for (size_t I = 0; I < nR; I++) {
106+
for (size_t J = I + 1; J < nR; J++) {
107+
double d = norm(R[I] - R[J]);
108+
dRR[I*nR + J] = d;
109+
dRR[J*nR + I] = d;
110+
}
111+
}
112+
return dRR;
113+
}
101114

102-
TEST(PartitionTest, Becke) {
115+
class PartitionTest: public ::testing::Test {
116+
protected:
117+
PartitionTest();
118+
119+
// grid & weight for one-center integration
120+
std::vector<double> r;
121+
std::vector<double> w;
122+
123+
const double tol = 1e-5;
124+
};
125+
126+
PartitionTest::PartitionTest() {
103127
// angular grid & weight
104128
std::vector<double> r_ang, w_ang;
105129
int lmax = 25;
@@ -114,7 +138,8 @@ TEST(PartitionTest, Becke) {
114138

115139
// complete grid & weight for one-center integration
116140
size_t ngrid = w_rad.size() * w_ang.size();
117-
std::vector<double> r(3*ngrid), w(ngrid);
141+
r.resize(3*ngrid);
142+
w.resize(ngrid);
118143

119144
size_t ir = 0;
120145
for (size_t i = 0; i < w_rad.size(); i++) {
@@ -126,53 +151,95 @@ TEST(PartitionTest, Becke) {
126151
++ir;
127152
}
128153
}
154+
}
155+
129156

130-
// each Param defines a test function
157+
TEST_F(PartitionTest, Becke) {
131158
for (const Param& param : test_params) {
132159
double val = 0.0;
133160
double val_ref = ref(param.a, param.n);
134161

135-
// tabulate dRR[I,J] = |R[I] - R[J]|
162+
// tabulate || R[I] - R[J] ||
163+
std::vector<double> dRR(dist_R_R(param.R));
164+
165+
// all centers are involved
136166
size_t nR = param.R.size();
137-
std::vector<double> dRR(nR*nR, 0.0);
138-
for (size_t I = 0; I < nR; I++) {
139-
for (size_t J = I + 1; J < nR; J++) {
140-
double d = norm(param.R[I] - param.R[J]);
141-
dRR[I*nR + J] = d;
142-
dRR[J*nR + I] = d;
167+
std::vector<int> iR(nR);
168+
std::iota(iR.begin(), iR.end(), 0);
169+
170+
for (size_t I = 0; I < nR; ++I) { // for each center
171+
for (size_t i = 0; i < w.size(); i++) {
172+
Vec3 ri = Vec3{r[3*i], r[3*i+1], r[3*i+2]} + param.R[I];
173+
174+
// tabulate || r - R[J] ||
175+
std::vector<double> drR(nR);
176+
for (size_t J = 0; J < nR; ++J) {
177+
drR[J] = norm(ri - param.R[J]);
178+
}
179+
180+
// partition weight for this grid point
181+
double w_part = Grid::Partition::w_becke(
182+
drR.size(), drR.data(), dRR.data(),
183+
iR.size(), iR.data(), I
184+
);
185+
186+
val += w_part * w[i] * func(ri, param.R, param.a, param.n);
143187
}
144188
}
145189

190+
EXPECT_NEAR(val, val_ref, tol);
191+
}
192+
}
193+
194+
195+
TEST_F(PartitionTest, Stratmann) {
196+
for (const Param& param : test_params) {
197+
double val = 0.0;
198+
double val_ref = ref(param.a, param.n);
199+
200+
// tabulate || R[I] - R[J] ||
201+
std::vector<double> dRR(dist_R_R(param.R));
202+
146203
// all centers are involved
204+
size_t nR = param.R.size();
147205
std::vector<int> iR(nR);
148206
std::iota(iR.begin(), iR.end(), 0);
149207

208+
// radii of exclusive zone
209+
using Grid::Partition::stratmann_a;
210+
std::vector<double> drR_thr(nR);
211+
for (size_t I = 0; I < nR; ++I) {
212+
double dRRmin = 1e100;
213+
for (size_t J = 0; J < nR; ++J) {
214+
if (J != I) {
215+
dRRmin = std::min(dRRmin, dRR[I*nR + J]);
216+
}
217+
}
218+
drR_thr[I] = 0.5 * (1-stratmann_a) * dRRmin;
219+
}
220+
150221
for (size_t I = 0; I < nR; ++I) { // for each center
151-
for (size_t i = 0; i < ngrid; i++) {
152-
Vec3 ri = {
153-
r[3*i] + param.R[I][0],
154-
r[3*i+1] + param.R[I][1],
155-
r[3*i+2] + param.R[I][2]
156-
};
157-
158-
// tabulate drR[J] = |r - R[J]|
222+
for (size_t i = 0; i < w.size(); i++) {
223+
Vec3 ri = Vec3{r[3*i], r[3*i+1], r[3*i+2]} + param.R[I];
224+
225+
// tabulate || r - R[J] ||
159226
std::vector<double> drR(nR);
160227
for (size_t J = 0; J < nR; ++J) {
161228
drR[J] = norm(ri - param.R[J]);
162229
}
163230

164231
// partition weight for this grid point
165-
double w_part = Grid::Partition::w_becke(
166-
drR.size(), drR.data(), dRR.data(), iR.size(), iR.data(), I
232+
double w_part = Grid::Partition::w_stratmann(
233+
drR.size(), drR.data(), dRR.data(), drR_thr.data(),
234+
iR.size(), iR.data(), I
167235
);
168236

169237
val += w_part * w[i] * func(ri, param.R, param.a, param.n);
170238
}
171239
}
172240

173-
EXPECT_NEAR(val, val_ref, 1e-5);
241+
EXPECT_NEAR(val, val_ref, tol);
174242
}
175-
176243
}
177244

178245
int main(int argc, char** argv)

0 commit comments

Comments
 (0)