18 static constexpr float radii[] = {
21 1.82f, 0.63f, 1.75f, 1.775f, 1.50f,
23 2.27f, 1.73f, 1.50f, 2.10f, 1.90f,
25 2.75f, 1.95f, 1.32f, 1.95f, 1.06f,
26 1.13f, 1.19f, 1.26f, 1.13f, 1.63f,
27 1.40f, 1.39f, 1.87f, 1.48f, 0.83f,
29 2.65f, 2.02f, 1.61f, 1.42f, 1.33f,
30 1.75f, 2.00f, 1.20f, 1.22f, 1.63f,
31 1.72f, 1.58f, 1.93f, 2.17f, 1.12f,
33 3.01f, 2.41f, 1.83f, 1.86f, 1.62f,
34 1.79f, 1.76f, 1.74f, 1.96f, 1.69f,
35 1.66f, 1.63f, 1.61f, 1.59f, 1.57f,
36 1.54f, 1.53f, 1.40f, 1.22f, 1.26f,
37 1.30f, 1.58f, 1.22f, 1.72f, 1.66f,
38 1.55f, 1.96f, 2.02f, 1.73f, 1.21f,
40 3.24f, 2.57f, 2.12f, 1.84f, 1.60f,
41 1.75f, 1.71f, 1.67f, 1.66f, 1.65f,
42 1.64f, 1.63f, 1.62f, 1.61f, 1.60f,
43 1.59f, 1.58f, 1.00f, 1.00f, 1.00f,
44 1.00f, 1.00f, 1.00f, 1.00f, 1.00f,
45 1.00f, 1.00f, 1.00f, 1.00f, 1.00f,
49 static_assert(radii[
static_cast<int>(
El::D)] == 1.2f,
"Hmm");
50 static_assert(
sizeof(radii) /
sizeof(radii[0]) ==
51 static_cast<int>(
El::END) + 1,
"Hmm");
52 return radii[
static_cast<int>(el)];
59 static constexpr float radii[] = {
62 0.53f, 0.21f, 0.05f, 1.90f, 1.12f,
64 0.93f, 0.51f, 0.33f, 0.20f, 0.39f,
66 1.31f, 0.94f, 0.69f, 0.36f, 0.48f,
67 0.33f, 0.26f, 0.48f, 0.34f, 0.43f,
68 0.51f, 0.54f, 0.41f, 0.20f, 0.28f,
70 1.28f, 1.12f, 0.84f, 0.53f, 0.42f,
71 0.35f, 0.31f, 0.32f, 0.49f, 0.58f,
72 0.61f, 0.72f, 0.56f, 0.49f, 0.70f,
74 1.61f, 1.29f, 0.97f, 0.81f, 0.79f,
75 0.92f, 0.91f, 0.90f, 0.89f, 0.88f,
76 0.70f, 0.85f, 0.84f, 0.83f, 0.82f,
77 0.81f, 0.80f, 0.52f, 0.58f, 0.36f,
78 0.32f, 0.33f, 0.51f, 0.51f, 0.51f,
79 0.90f, 0.69f, 0.59f, 0.70f, 0.61f,
81 1.74f, 1.42f, 1.06f, 0.88f, 0.72f,
82 0.46f, 0.65f, 0.65f, 0.79f, 0.79f,
83 0.77f, 0.76f, 1.00f, 1.00f, 1.00f,
84 1.00f, 1.00f, 1.00f, 1.00f, 1.00f,
85 1.00f, 1.00f, 1.00f, 1.00f, 1.00f,
86 1.00f, 1.00f, 1.00f, 1.00f, 1.00f,
90 static_assert(radii[
static_cast<int>(
El::D)] == 1.40f,
"Hmm");
91 return radii[
static_cast<int>(el)];
95 case El::H:
return 1.4f;
96 case El::D:
return 1.4f;
97 case El::O:
return 1.08f;
98 case El::C:
return 2.0f;
99 case El::N:
return 1.12f;
100 default:
return 1.6f;
148 int du = (
int) std::floor(
r / mask.spacing[0]);
149 int dv = (
int) std::floor(
r / mask.spacing[1]);
150 int dw = (
int) std::floor(
r / mask.spacing[2]);
151 double max_spacing2 =
sq(std::max(std::max(mask.unit_cell.a / mask.nu,
152 mask.unit_cell.b / mask.nv),
153 mask.unit_cell.c / mask.nw)) + 1
e-6;
154 if (2 *
du >= mask.nu || 2 *
dv >= mask.nv || 2 *
dw >= mask.nw)
155 fail(
"grid operation failed: radius bigger than half the unit cell?");
156 std::vector<std::array<int,3>>
stencil1;
157 std::vector<std::array<int,3>>
stencil2;
158 for (
int w = -
dw; w <=
dw; ++w)
159 for (
int v = -
dv; v <=
dv; ++v)
160 for (
int u = -
du; u <=
du; ++u) {
162 double r2 = mask.unit_cell.orthogonalize_difference(
fdelta).length_sq();
163 if (
r2 <=
r *
r &&
r2 != 0.) {
164 std::array<int,3>
wvu{{w <= 0 ? w : w - mask.nw,
165 v <= 0 ? v : v - mask.nv,
166 u <= 0 ? u : u - mask.nu}};
177 size_t idx = mask.index_near_zero(
p.u +
wvu[2],
p.v +
wvu[1],
p.w +
wvu[0]);
178 if (mask.data[idx] == value) {
189 size_t idx = mask.index_near_zero(
p.u +
wvu[2],
p.v +
wvu[1],
p.w +
wvu[0]);
190 if (mask.data[idx] != value) {
197 size_t idx = mask.index_near_zero(
p.u +
wvu[2],
p.v +
wvu[1],
p.w +
wvu[0]);
198 if (mask.data[idx] != value)