Skip to content

Commit 8d4c15e

Browse files
jinzx10jinzx10
committed
becke partition initial commit
1 parent ba0cb21 commit 8d4c15e

File tree

5 files changed

+417
-10
lines changed

5 files changed

+417
-10
lines changed

‎source/module_base/grid/partition.cpp‎

Lines changed: 124 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,124 @@
1+
#include "module_base/grid/partition.h"
2+
#include "module_base/constants.h"
3+
4+
#include <cmath>
5+
#include <numeric>
6+
#include <vector>
7+
8+
namespace Grid {
9+
namespace Partition {
10+
11+
double w_becke(
12+
int nR0,
13+
double* drR,
14+
double* dRR,
15+
int nR,
16+
int* iR,
17+
int c
18+
) {
19+
std::vector<double> P(nR, 1.0);
20+
for (int i = 0; i < nR; ++i) {
21+
int I = iR[i];
22+
for (int j = i + 1; j < nR; ++j) {
23+
int J = iR[j];
24+
double mu = (drR[I] - drR[J]) / dRR[I*nR0 + J];
25+
double s = s_becke(mu);
26+
P[I] *= s;
27+
P[J] *= (1.0 - s); // s_becke(-mu) = 1 - s_becke(mu)
28+
}
29+
}
30+
return P[c] / std::accumulate(P.begin(), P.end(), 0.0);
31+
}
32+
33+
double s_becke(double mu) {
34+
/*
35+
* Becke's iterated polynomials (3rd order)
36+
*
37+
* s(mu) = 0.5 * (1 - p(p(p(mu))))
38+
*
39+
* p(x) = 0.5 * x * (3 - x^2)
40+
*
41+
*/
42+
double p = 0.5 * mu * (3.0 - mu*mu);
43+
p = 0.5 * p * (3.0 - p*p);
44+
p = 0.5 * p * (3.0 - p*p);
45+
return 0.5 * (1.0 - p);
46+
}
47+
48+
49+
double s_stratmann(double mu, double a) {
50+
/*
51+
* Stratmann's piecewise cell function
52+
*
53+
* s(mu) = 0.5 * (1 - g(mu/a))
54+
*
55+
* / -1 x <= -1
56+
* |
57+
* g(x) = | (35x - 35x^3 + 21x^5 - 5x^7) / 16 |x| < 1
58+
* |
59+
* \ +1 x >= +1
60+
*
61+
*/
62+
double x = mu / a;
63+
double x2 = x * x;
64+
double h = 0.625 * x * (35 + x2 * (-35 + x2 * (21 - 5 * x2)));
65+
66+
bool mid = std::abs(x) < 1;
67+
double g = !mid * (1 - 2 * std::signbit(x)) + mid * h;
68+
return 0.5 * (1 - g);
69+
}
70+
71+
72+
double u_stratmann_mod(double y, double b) {
73+
using ModuleBase::PI;
74+
bool core = y <= b;
75+
bool edge = !core && y < 1.0;
76+
return core + edge * 0.5 * (std::cos(PI * (y - b) / (1.0 - b)) + 1.0);
77+
}
78+
79+
80+
double s_stratmann_mod(double mu, double y, double a, double b) {
81+
// Modified Stratmann's cell function by Knuth et al.
82+
// y = |r-R(J)| / Rcut(J)
83+
return 1.0 + u_stratmann_mod(y, b) * (s_stratmann(mu, a) - 1.0);
84+
}
85+
86+
87+
//double weight(
88+
// int nRtot,
89+
// double* drR,
90+
// double* dRR,
91+
// int nR,
92+
// int* iR,
93+
// int c,
94+
// CellFuncType type,
95+
// double* Rcut
96+
//) {
97+
// // unnormalized weight
98+
// std::vector<double> P(nR, 1.0);
99+
//
100+
// // confocal ellipsoidal coordinates
101+
// std::vector<double> mu(nR*nR, 0.0);
102+
// for (int i = 0; i < nR; ++i) {
103+
// int I = iR[i];
104+
// for (int j = i + 1; j < nR; ++j) {
105+
// int J = iR[j];
106+
// mu[I*nR + J] = (drR[I] - drR[J]) / dRR[I*nRtot + J];
107+
// mu[J*nR + I] = -mu[I*nR + J];
108+
// }
109+
// }
110+
//
111+
// for (int i = 0; i < nR; ++i) {
112+
// if (i == c) {
113+
// continue;
114+
// }
115+
// int J = iR[i];
116+
// double mu = (drR[I] - drR[J]) / dRR[I*nRtot + J];
117+
// P[i] *= s_becke(mu);
118+
// }
119+
//
120+
// return P[c] / std::accumulate(P.begin(), P.end(), 0.0);
121+
//}
122+
123+
} // end of namespace Partition
124+
} // end of namespace Grid

‎source/module_base/grid/partition.h‎

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
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+
StratmannMod
11+
};
12+
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+
double* drR,
41+
double* dRR,
42+
int nR,
43+
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 dRcRnn Distance between the grid center and its nearest neighbor.
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+
double* drR,
71+
double* dRR,
72+
int nR,
73+
int* iR,
74+
int c,
75+
double dRcRnn,
76+
double a = 0.64
77+
);
78+
79+
// Stratmann's piecewise cell function
80+
double s_stratmann(double mu, double a);
81+
82+
83+
/**
84+
* Knuth, F., Carbogno, C., Atalla, V., Blum, V., & Scheffler, M. (2015).
85+
* All-electron formalism for total energy strain derivatives and stress
86+
* tensor components for numeric atom-centered orbitals.
87+
* Computer Physics Communications, 190, 33-50.
88+
*/
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
91+
92+
93+
} // end of namespace Partition
94+
} // end of namespace Grid
95+
96+
#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/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)