SIRIUS 7.5.0
Electronic structure library and applications
gvec.cpp
Go to the documentation of this file.
1// Copyright (c) 2013-2017 Anton Kozhevnikov, Thomas Schulthess
2// All rights reserved.
3//
4// Redistribution and use in source and binary forms, with or without modification, are permitted provided that
5// the following conditions are met:
6//
7// 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the
8// following disclaimer.
9// 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions
10// and the following disclaimer in the documentation and/or other materials provided with the distribution.
11//
12// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
13// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
14// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
15// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
16// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
17// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
18// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
19
20/** \file gvec.cpp
21 *
22 * \brief Contains the implementation of Gvec class.
23 *
24 */
25
26#include "symmetry/lattice.hpp"
27#include "gvec.hpp"
28#include "core/serializer.hpp"
29
30namespace sirius {
31
32namespace fft {
33
35{
36 /* index of the z coordinate of G-vector: first 12 bits */
37 uint32_t j = idx__ & 0xFFF;
38 /* index of z-column: last 20 bits */
39 uint32_t i = idx__ >> 12;
40 RTE_ASSERT(i < (uint32_t)z_columns_.size());
41 RTE_ASSERT(j < (uint32_t)z_columns_[i].z.size());
42 int x = z_columns_[i].x;
43 int y = z_columns_[i].y;
44 int z = z_columns_[i].z[j];
45 return r3::vector<int>(x, y, z);
46}
47
48void Gvec::find_z_columns(double Gmax__, fft::Grid const& fft_box__)
49{
50 PROFILE("fft::Gvec::find_z_columns");
51
52 sddk::mdarray<int, 2> non_zero_columns(fft_box__.limits(0), fft_box__.limits(1));
53 non_zero_columns.zero();
54
55 num_gvec_ = 0;
56
57 auto add_new_column = [&](int i, int j) {
58
59 if (non_zero_columns(i, j)) {
60 return;
61 }
62
63 std::vector<int> zcol;
64
65 /* in general case take z in [0, Nz) */
66 int zmax = fft_box__[2] - 1;
67 /* in case of G-vector reduction take z in [0, Nz/2] for {x=0,y=0} stick */
68 if (reduce_gvec_ && !i && !j) {
69 zmax = fft_box__.limits(2).second;
70 }
71 /* loop over z-coordinates of FFT grid */
72 for (int iz = 0; iz <= zmax; iz++) {
73 /* get z-coordinate of G-vector */
74 int k = fft_box__.freq_by_coord<2>(iz);
75 /* take G+k */
76 auto vgk = r3::dot(lattice_vectors_, (r3::vector<double>(i, j, k) + vk_));
77 /* add z-coordinate of G-vector to the list */
78 if (vgk.length() <= Gmax__) {
79 zcol.push_back(k);
80 }
81 }
82
83 /* add column to the list */
84 if (zcol.size()) {
85 z_columns_.push_back(z_column_descriptor(i, j, zcol));
86 num_gvec_ += static_cast<int>(zcol.size());
87
88 non_zero_columns(i, j) = 1;
89 if (reduce_gvec_) {
90 int mi = -i;
91 int mj = -j;
92 if (mi >= fft_box__.limits(0).first && mi <= fft_box__.limits(0).second &&
93 mj >= fft_box__.limits(1).first && mj <= fft_box__.limits(1).second) {
94 non_zero_columns(mi, mj) = 1;
95 }
96 }
97 }
98 };
99
100 PROFILE_START("fft::Gvec::find_z_columns|add");
101 /* copy column order from previous G-vector set */
102 if (gvec_base_) {
103 for (int icol = 0; icol < gvec_base_->num_zcol(); icol++) {
104 int i = gvec_base_->zcol(icol).x;
105 int j = gvec_base_->zcol(icol).y;
106 add_new_column(i, j);
107 }
108 }
109
110 /* check all z-columns and add if within sphere. Only allow non-negative x-indices for reduced case */
111 for (int i = reduce_gvec_? 0 : fft_box__.limits(0).first; i <= fft_box__.limits(0).second; i++) {
112 for (int j = fft_box__.limits(1).first; j <= fft_box__.limits(1).second; j++) {
113 add_new_column(i, j);
114 }
115 }
116 PROFILE_STOP("fft::Gvec::find_z_columns|add");
117
118 if (!gvec_base_) {
119 /* put column with {x, y} = {0, 0} to the beginning */
120 for (size_t i = 0; i < z_columns_.size(); i++) {
121 if (z_columns_[i].x == 0 && z_columns_[i].y == 0) {
122 std::swap(z_columns_[i], z_columns_[0]);
123 break;
124 }
125 }
126 }
127
128 PROFILE_START("fft::Gvec::find_z_columns|sym");
129 /* now we have to remove edge G-vectors that don't form a complete shell */
130 if (bare_gvec_) {
131 auto lat_sym = sirius::find_lat_sym(this->unit_cell_lattice_vectors(), this->sym_tol_);
132
133 std::fill(non_zero_columns.at(sddk::memory_t::host), non_zero_columns.at(sddk::memory_t::host) + non_zero_columns.size(), -1);
134 for (int i = 0; i < static_cast<int>(z_columns_.size()); i++) {
135 non_zero_columns(z_columns_[i].x, z_columns_[i].y) = i;
136 }
137
138 std::vector<z_column_descriptor> z_columns_tmp;
139
140 auto remove_incomplete = [this, &z_columns_tmp, &lat_sym, &non_zero_columns](int i) // i - index of column
141 {
142 int z_min = (2 << 20);
143 int z_max = -(2 << 20);
144 for (int iz = 0; iz < static_cast<int>(z_columns_[i].z.size()); iz++) {
145 z_min = std::min(z_min, z_columns_[i].z[iz]);
146 z_max = std::max(z_max, z_columns_[i].z[iz]);
147 }
148 std::vector<int> z;
149 for (int iz = 0; iz < static_cast<int>(z_columns_[i].z.size()); iz++) {
150 bool found_for_all_sym{true};
151 /* check only first or last z coordinate inside z column */
152 if (z_columns_[i].z[iz] == z_min || z_columns_[i].z[iz] == z_max) {
153 r3::vector<int> G(z_columns_[i].x, z_columns_[i].y, z_columns_[i].z[iz]);
154 for (auto& R: lat_sym) {
155 /* apply lattice symmeetry operation to a G-vector */
156 auto G1 = r3::dot(G, R);
157 if (reduce_gvec_) {
158 if (G1[0] == 0 && G1[1] == 0) {
159 G1[2] = std::abs(G1[2]);
160 }
161 }
162 int i1 = non_zero_columns(G1[0], G1[1]);
163 if (i1 == -1) {
164 G1 = G1 * (-1);
165 i1 = non_zero_columns(G1[0], G1[1]);
166 if (i1 == -1) {
167 std::stringstream s;
168 s << "index of z-column is not found" << std::endl
169 << " G : " << G << std::endl
170 << " G1 : " << G1;
171 RTE_THROW(s);
172 }
173 }
174
175 bool found = (std::find(z_columns_[i1].z.begin(), z_columns_[i1].z.end(), G1[2]) != std::end(z_columns_[i1].z));
176 found_for_all_sym = found_for_all_sym && found;
177 } // R
178 }
179 if (found_for_all_sym) {
180 z.push_back(z_columns_[i].z[iz]);
181 }
182 } // iz
183 if (z.size()) {
184 z_columns_tmp.push_back(z_column_descriptor(z_columns_[i].x, z_columns_[i].y, z));
185 }
186 };
187
188 for (int i = 0; i < static_cast<int>(z_columns_.size()); i++) {
189 remove_incomplete(i);
190 }
191 z_columns_ = z_columns_tmp;
192
193 num_gvec_ = 0;
194 for (auto& zc : z_columns_) {
195 num_gvec_ += static_cast<int>(zc.z.size());
196 }
197 }
198 PROFILE_STOP("fft::Gvec::find_z_columns|sym");
199
200 PROFILE_START("fft::Gvec::find_z_columns|sort");
201 /* sort z-columns starting from the second or skip num_zcol of base distribution */
202 int n = (gvec_base_) ? gvec_base_->num_zcol() : 1;
203 std::sort(z_columns_.begin() + n, z_columns_.end(),
204 [](z_column_descriptor const& a, z_column_descriptor const& b) { return a.z.size() > b.z.size(); });
205 PROFILE_STOP("fft::Gvec::find_z_columns|sort");
206}
207
209{
212 /* local number of z-columns for each rank */
213 std::vector<std::vector<z_column_descriptor>> zcols_local(comm().size());
214
215 /* use already existing distribution of base G-vector set */
216 if (gvec_base_) {
217 for (int rank = 0; rank < comm().size(); rank++) {
218 for (int i = 0; i < gvec_base_->zcol_count(rank); i++) {
219 int icol = gvec_base_->zcol_offset(rank) + i;
220 /* assign column to the found rank */
221 zcols_local[rank].push_back(z_columns_[icol]);
222 /* count local number of z-columns */
223 zcol_distr_.counts[rank] += 1;
224 /* count local number of G-vectors */
225 gvec_distr_.counts[rank] += static_cast<int>(z_columns_[icol].z.size());
226 }
227 }
228 }
229
230 int n = (gvec_base_) ? gvec_base_->num_zcol() : 0;
231
232 std::vector<int> ranks;
233 for (int i = n; i < static_cast<int>(z_columns_.size()); i++) {
234 /* initialize the list of ranks to 0,1,2,... */
235 if (ranks.empty()) {
236 ranks.resize(comm().size());
237 std::iota(ranks.begin(), ranks.end(), 0);
238 }
239 /* find rank with minimum number of G-vectors */
240 auto rank_with_min_gvec = std::min_element(ranks.begin(), ranks.end(), [this](const int& a, const int& b) {
241 return gvec_distr_.counts[a] < gvec_distr_.counts[b];
242 });
243
244 /* assign column to the found rank */
245 zcols_local[*rank_with_min_gvec].push_back(z_columns_[i]);
246 /* count local number of z-columns */
247 zcol_distr_.counts[*rank_with_min_gvec] += 1;
248 /* count local number of G-vectors */
249 gvec_distr_.counts[*rank_with_min_gvec] += static_cast<int>(z_columns_[i].z.size());
250 /* exclude this rank from the search */
251 ranks.erase(rank_with_min_gvec);
252 }
253 gvec_distr_.calc_offsets();
254 zcol_distr_.calc_offsets();
255
256 /* save new ordering of z-columns */
257 z_columns_.clear();
258 for (int rank = 0; rank < comm().size(); rank++) {
259 z_columns_.insert(z_columns_.end(), zcols_local[rank].begin(), zcols_local[rank].end());
260 }
261
262 /* sanity check */
263 int ng{0};
264 for (int rank = 0; rank < comm().size(); rank++) {
265 ng += gvec_distr_.counts[rank];
266 }
267 if (ng != num_gvec_) {
268 RTE_THROW("wrong number of G-vectors");
269 }
270 this->offset_ = this->gvec_offset(this->comm().rank());
271 this->count_ = this->gvec_count(this->comm().rank());
272 this->num_zcol_local_ = this->zcol_distr_.counts[this->comm().rank()];
273}
274
276{
277 if (!bare_gvec_) {
278 return;
279 }
280
281 PROFILE("fft::Gvec::find_gvec_shells");
282
283 auto lat_sym = sirius::find_lat_sym(this->unit_cell_lattice_vectors(), this->sym_tol_);
284
286 gvec_shell_ = sddk::mdarray<int, 1>(num_gvec_, sddk::memory_t::host, "gvec_shell_");
287
288 std::fill(&gvec_shell_[0], &gvec_shell_[0] + num_gvec_, -1);
289
290 /* find G-vector shells using symmetry consideration */
291 for (int ig = 0; ig < num_gvec_; ig++) {
292 /* if the shell for this vector is not yet found */
293 if (gvec_shell_[ig] == -1) {
294 auto G = gvec<index_domain_t::global>(ig);
295 for (auto& R: lat_sym) {
296 auto G1 = r3::dot(G, R);
297 auto ig1 = index_by_gvec(G1);
298 if (ig1 == -1) {
299 G1 = G1 * (-1);
300 ig1 = index_by_gvec(G1);
301 if (ig1 == -1) {
302 RTE_THROW("symmetry-related G-vector is not found");
303 }
304 }
305 if (gvec_shell_[ig1] == -1) {
307 } else {
308 if (gvec_shell_[ig1] != num_gvec_shells_) {
309 auto gc = r3::dot(lattice_vectors_, G);
310 auto gc1 = r3::dot(lattice_vectors_, G1);
311 std::stringstream s;
312 s << "Error in G-vector shell index" << std::endl
313 << " G : " << G << std::endl
314 << " rotated G : " << G1 << std::endl
315 << " current shell index : " << num_gvec_shells_ << std::endl
316 << " rotated G shell index : " << gvec_shell_[ig1] << std::endl
317 << " length of G : " << gc.length() << std::endl
318 << " length of rotated G : " << gc1.length() << std::endl
319 << " length difference : " << std::abs(gc.length() - gc1.length());
320 RTE_THROW(s);
321 }
322 }
323 }
325 }
326 }
327 for (int ig = 0; ig < num_gvec_; ig++) {
328 if (gvec_shell_[ig] == -1) {
329 RTE_THROW("wrong G-vector shell");
330 }
331 }
332
333 gvec_shell_len_ = sddk::mdarray<double, 1>(num_gvec_shells_, sddk::memory_t::host, "gvec_shell_len_");
334 std::fill(&gvec_shell_len_[0], &gvec_shell_len_[0] + num_gvec_shells_, 0);
335
336 std::vector<int> ngv_sh(num_gvec_shells_, 0);
337
338 for (int ig = 0; ig < num_gvec_; ig++) {
339 auto g = gvec_cart<index_domain_t::global>(ig).length();
340 int igsh = gvec_shell_[ig];
341 gvec_shell_len_[igsh] += g;
342 ngv_sh[igsh]++;
343 }
344 for (int i = 0; i < num_gvec_shells_; i++) {
345 gvec_shell_len_[i] /= ngv_sh[i];
346 }
347
348 /* map from global index of G-shell to a list of local G-vectors */
349 std::map<int, std::vector<int>> gshmap;
350 for (int igloc = 0; igloc < this->count(); igloc++) {
351 int igsh = this->shell(this->offset() + igloc);
352 if (gshmap.count(igsh) == 0) {
353 gshmap[igsh] = std::vector<int>();
354 }
355 gshmap[igsh].push_back(igloc);
356 }
357
359 gvec_shell_idx_local_.resize(this->count());
360 gvec_shell_len_local_.clear();
361 for (auto it = gshmap.begin(); it != gshmap.end(); ++it) {
362 int igsh = it->first;
363 gvec_shell_len_local_.push_back(this->shell_len(igsh));
364 for (auto igloc: it->second) {
366 }
368 }
369}
370
371void
373{
374 gvec_ = sddk::mdarray<int, 2>(3, count(), sddk::memory_t::host, "gvec_");
375 gkvec_ = sddk::mdarray<double, 2>(3, count(), sddk::memory_t::host, "gkvec_");
376
377 for (int igloc = 0; igloc < count(); igloc++) {
378 int ig = offset() + igloc;
380 for (int x : {0, 1, 2}) {
381 gvec_(x, igloc) = G[x];
382 gkvec_(x, igloc) = G[x] + vk_[x];
383 }
384 }
385}
386
387void
389{
390 gvec_cart_ = sddk::mdarray<double, 2>(3, count(), sddk::memory_t::host, "gvec_cart_");
391 gkvec_cart_ = sddk::mdarray<double, 2>(3, count(), sddk::memory_t::host, "gkvec_cart_");
392 /* this arrays are allocated with GPU- friendly data layout */
393 gvec_tp_ = sddk::mdarray<double, 2>(count(), 2, sddk::memory_t::host, "gvec_tp_");
394 gkvec_tp_ = sddk::mdarray<double, 2>(count(), 2, sddk::memory_t::host, "gvec_tp_");
395 if (bare_gvec_) {
396 gvec_len_ = sddk::mdarray<double, 1>(count(), sddk::memory_t::host, "gvec_len_");
397 }
398
399 for (int igloc = 0; igloc < count(); igloc++) {
400 auto gc = r3::dot(lattice_vectors_, r3::vector<int>(&gvec_(0, igloc)));
401 auto gkc = r3::dot(lattice_vectors_, r3::vector<double>(&gkvec_(0, igloc)));
402 for (int x : {0, 1, 2}) {
403 gvec_cart_(x, igloc) = gc[x];
404 gkvec_cart_(x, igloc) = gkc[x];
405 }
406 if (bare_gvec_) {
407 gvec_len_(igloc) = gvec_shell_len_(gvec_shell_(this->offset() + igloc));
408 }
409 auto gs = r3::spherical_coordinates(gc);
410 gvec_tp_(igloc, 0) = gs[1];
411 gvec_tp_(igloc, 1) = gs[2];
412
413 auto gks = r3::spherical_coordinates(gkc);
414 gkvec_tp_(igloc, 0) = gks[1];
415 gkvec_tp_(igloc, 1) = gks[2];
416 }
417}
418
419void
420Gvec::init(fft::Grid const& fft_grid)
421{
422 PROFILE("fft::Gvec::init");
423
424 find_z_columns(Gmax_, fft_grid);
425
427
428 gvec_index_by_xy_ =
429 sddk::mdarray<int, 3>(2, fft_grid.limits(0), fft_grid.limits(1), sddk::memory_t::host, "Gvec.gvec_index_by_xy_");
430 std::fill(gvec_index_by_xy_.at(sddk::memory_t::host), gvec_index_by_xy_.at(sddk::memory_t::host) + gvec_index_by_xy_.size(),
431 -1);
432
433 /* build the full G-vector index and reverse mapping */
435 int ig{0};
436 for (size_t i = 0; i < z_columns_.size(); i++) {
437 /* starting G-vector index for a z-stick */
438 gvec_index_by_xy_(0, z_columns_[i].x, z_columns_[i].y) = ig;
439 /* pack size of a z-stick and column index in one number */
440 gvec_index_by_xy_(1, z_columns_[i].x, z_columns_[i].y) = static_cast<int>((z_columns_[i].z.size() << 20) + i);
441 for (size_t j = 0; j < z_columns_[i].z.size(); j++) {
442 gvec_full_index_[ig++] = static_cast<uint32_t>((i << 12) + j);
443 }
444 }
445 if (ig != num_gvec_) {
446 RTE_THROW("wrong G-vector count");
447 }
448
449 for (int ig = 0; ig < num_gvec_; ig++) {
450 auto gv = gvec<index_domain_t::global>(ig);
451 if (index_by_gvec(gv) != ig) {
452 std::stringstream s;
453 s << "wrong G-vector index: ig=" << ig << " gv=" << gv << " index_by_gvec(gv)=" << index_by_gvec(gv);
454 RTE_THROW(s);
455 }
456 }
457
458 /* first G-vector must be (0, 0, 0); never remove this check!!! */
460 if (g0[0] || g0[1] || g0[2]) {
461 RTE_THROW("first G-vector is not zero");
462 }
463
465
466 if (gvec_base_) {
467 /* the size of the mapping is equal to the local number of G-vectors in the base set */
468 gvec_base_mapping_ = sddk::mdarray<int, 1>(gvec_base_->count(), sddk::memory_t::host, "gvec_base_mapping_");
469 /* loop over local G-vectors of a base set */
470 for (int igloc = 0; igloc < gvec_base_->count(); igloc++) {
471 /* G-vector in lattice coordinates */
472 auto G = gvec_base_->gvec<index_domain_t::local>(igloc);
473 /* global index of G-vector in the current set */
474 int ig = index_by_gvec(G);
475 /* the same MPI rank must store this G-vector */
476 ig -= offset();
477 if (ig >= 0 && ig < count()) {
478 gvec_base_mapping_(igloc) = ig;
479 } else {
480 std::stringstream s;
481 s << "local G-vector index is not found" << std::endl
482 << " G-vector: " << G << std::endl
483 << " G-vector index in base distribution : " << gvec_base_->offset() + igloc << std::endl
484 << " G-vector index in base distribution (by G-vector): " << gvec_base_->index_by_gvec(G) << std::endl
485 << " G-vector index in new distribution : " << index_by_gvec(G) << std::endl
486 << " offset in G-vector index for this rank: " << offset() << std::endl
487 << " local number of G-vectors for this rank: " << count();
488 RTE_THROW(s);
489 }
490 }
491 }
492 // TODO: add a check for gvec_base (there is already a test for this).
495}
496
497std::pair<int, bool> Gvec::index_g12_safe(r3::vector<int> const& g1__, r3::vector<int> const& g2__) const
498{
499 auto v = g1__ - g2__;
500 int idx = index_by_gvec(v);
501 bool conj{false};
502 if (idx < 0) {
503 idx = index_by_gvec(v * (-1));
504 conj = true;
505 }
506 if (idx < 0 || idx >= num_gvec()) {
507 std::stringstream s;
508 s << "wrong index of G-G' vector" << std::endl
509 << " G: " << g1__ << std::endl
510 << " G': " << g2__ << std::endl
511 << " G - G': " << v << std::endl
512 << " idx: " << idx;
513 RTE_THROW(s);
514 }
515 return std::make_pair(idx, conj);
516}
517
519{
520 /* reduced G-vector set does not have negative z for x=y=0 */
521 if (reduced() && G__[0] == 0 && G__[1] == 0 && G__[2] < 0) {
522 return -1;
523 }
524 int ig0 = gvec_index_by_xy_(0, G__[0], G__[1]);
525 if (ig0 == -1) {
526 return -1;
527 }
528 /* index of the column */
529 int icol = gvec_index_by_xy_(1, G__[0], G__[1]) & 0xFFFFF;
530 /* quick exit if z is out of bounds */
531 if (G__[2] < z_columns_[icol].z_min || G__[2] > z_columns_[icol].z_max) {
532 return -1;
533 }
534 /* size of the column */
535 int col_size = gvec_index_by_xy_(1, G__[0], G__[1]) >> 20;
536
537 /* three possible options for the z-column location
538
539 frequency ... -4, -3, -2, -1, 0, 1, 2, 3, 4 ...
540 -----------------------------------------------------------------------------
541 G-vec ordering
542 #1 (all negative) ___ 0 1 2 3 __________________
543 #2 (negative and positive) ____________ 3 4 0 1 2 _________
544 #3 (all positive) _____________________ 0 1 2 3 ___
545
546 Remember how FFT frequencies are stored: first positive frequences, then negative in the reverse order
547
548 subtract first z-coordinate in column from the current z-coordinate of G-vector: in case #1 or #3 this
549 already gives a proper offset, in case #2 storage of FFT frequencies must be taken into account
550 */
551 int z0 = G__[2] - z_columns_[icol].z[0];
552 /* calculate proper offset */
553 int offs = (z0 >= 0) ? z0 : z0 + col_size;
554 /* full index */
555 int ig = ig0 + offs;
556 RTE_ASSERT(ig >= 0 && ig < num_gvec());
557 return ig;
558}
559
560Gvec send_recv(mpi::Communicator const& comm__, Gvec const& gv_src__, int source__, int dest__)
561{
562 serializer s;
563
564 if (comm__.rank() == source__) {
565 ::sirius::fft::serialize(s, gv_src__);
566 }
567
568 s.send_recv(comm__, source__, dest__);
569
570 Gvec gv(gv_src__.comm());
571
572 if (comm__.rank() == dest__) {
573 ::sirius::fft::deserialize(s, gv);
574 }
575 return gv;
576}
577
578void Gvec_fft::build_fft_distr()
579{
580 /* calculate distribution of G-vectors and z-columns for the FFT communicator */
581 gvec_distr_fft_ = mpi::block_data_descriptor(comm_fft().size());
582
583 for (int rank = 0; rank < comm_fft().size(); rank++) {
584 for (int i = 0; i < comm_ortho_fft().size(); i++) {
585 /* fine-grained rank */
586 int r = rank_map_(rank, i);
587 gvec_distr_fft_.counts[rank] += gvec().gvec_count(r);
588 }
589 }
590 for (int i = 0; i < comm_ortho_fft().size(); i++) {
591 /* fine-grained rank */
592 int r = rank_map_(comm_fft().rank(), i);
594 }
595 /* get offsets of G-vectors */
596 gvec_distr_fft_.calc_offsets();
597}
598
600{
601 /* build a table of {offset, count} values for G-vectors in the swapped distribution;
602 * we are preparing to swap plane-wave coefficients from a default slab distribution to a FFT-friendly
603 * distribution
604 * +--------------+ +----+----+----+
605 * | : : | | | | |
606 * +--------------+ |....|....|....|
607 * | : : | -> | | | |
608 * +--------------+ |....|....|....|
609 * | : : | | | | |
610 * +--------------+ +----+----+----+
611 *
612 * i.e. we will make G-vector slabs more fat (pile-of-slabs) and at the same time reshulffle wave-functions
613 * between columns of the 2D MPI grid */
615 for (int i = 0; i < comm_ortho_fft_.size(); i++) {
617 }
618 gvec_fft_slab_.calc_offsets();
619
620 RTE_ASSERT(gvec_fft_slab_.offsets.back() + gvec_fft_slab_.counts.back() == gvec_distr_fft_.counts[comm_fft().rank()]);
621
624 for (int i = 0; i < comm_ortho_fft_.size(); i++) {
625 for (int j = 0; j < comm_fft_.size(); j++) {
626 int r = rank_map_(j, i);
627 /* get array of G-vectors of rank r */
628 auto gv = this->gvec_.gvec_local(r);
629
630 if (j == comm_fft_.rank()) {
631 for (int ig = 0; ig < gvec_fft_slab_.counts[i]; ig++) {
632 for (int x : {0, 1, 2}) {
633 gvec_array_(x, gvec_fft_slab_.offsets[i] + ig) = gv(x, ig);
634 }
635 }
636 }
637 }
638 }
640}
641
642Gvec_fft::Gvec_fft(Gvec const& gvec__, mpi::Communicator const& comm_fft__, mpi::Communicator const& comm_ortho_fft__)
643 : gvec_(gvec__)
644 , comm_fft_(comm_fft__)
645 , comm_ortho_fft_(comm_ortho_fft__)
646{
647 if (comm_fft_.size() * comm_ortho_fft_.size() != gvec_.comm().size()) {
648 std::stringstream s;
649 s << "wrong size of communicators" << std::endl
650 << " comm_fft_.size() = " << comm_fft_.size() << std::endl
651 << " comm_ortho_fft_.size() = " << comm_ortho_fft_.size() << std::endl
652 << " gvec_.comm().size() = " << gvec_.comm().size();
653 RTE_THROW(s);
654 }
656 rank_map_.zero();
657 /* get a global rank */
659 gvec_.comm().allreduce(&rank_map_(0, 0), gvec_.comm().size());
660
661 build_fft_distr();
662
663 pile_gvec();
664}
665
666Gvec_shells::Gvec_shells(Gvec const& gvec__)
667 : comm_(gvec__.comm())
668 , gvec_(gvec__)
669{
670 PROFILE("fft::Gvec_shells");
671
672 a2a_send_ = mpi::block_data_descriptor(comm_.size());
673 a2a_recv_ = mpi::block_data_descriptor(comm_.size());
674
675 /* split G-vector shells between ranks in cyclic order */
676 spl_num_gsh_ = splindex_block_cyclic<>(gvec_.num_shells(), n_blocks(comm_.size()), block_id(comm_.rank()), 1);
677
678 /* each rank sends a fraction of its local G-vectors to other ranks */
679 /* count this fraction */
680 for (int igloc = 0; igloc < gvec_.count(); igloc++) {
681 int ig = gvec_.offset() + igloc;
682 int igsh = gvec_.shell(ig);
683 a2a_send_.counts[spl_num_gsh_.location(igsh).ib]++;
684 }
685 a2a_send_.calc_offsets();
686 /* sanity check: total number of elements to send is equal to the local number of G-vector */
687 if (a2a_send_.size() != gvec_.count()) {
688 RTE_THROW("wrong number of G-vectors");
689 }
690 /* count the number of elements to receive */
691 for (int r = 0; r < comm_.size(); r++) {
692 for (int igloc = 0; igloc < gvec_.gvec_count(r); igloc++) {
693 /* index of G-vector in the original distribution */
694 int ig = gvec_.gvec_offset(r) + igloc;
695 /* index of the G-vector shell */
696 int igsh = gvec_.shell(ig);
697 if (spl_num_gsh_.location(igsh).ib == comm_.rank()) {
698 a2a_recv_.counts[r]++;
699 }
700 }
701 }
702 a2a_recv_.calc_offsets();
703 /* sanity check: sum of local sizes in the remapped order is equal to the total number of G-vectors */
704 int ng = gvec_count_remapped();
705 comm_.allreduce(&ng, 1);
706 if (ng != gvec_.num_gvec()) {
707 RTE_THROW("wrong number of G-vectors");
708 }
709
710 /* local set of G-vectors in the remapped order */
711 gvec_remapped_ = sddk::mdarray<int, 2>(3, gvec_count_remapped(), sddk::memory_t::host, "gvec_remapped_");
712 gvec_shell_remapped_ = sddk::mdarray<int, 1>(gvec_count_remapped(), sddk::memory_t::host, "gvec_shell_remapped_");
713 std::vector<int> counts(comm_.size(), 0);
714 for (int r = 0; r < comm_.size(); r++) {
715 for (int igloc = 0; igloc < gvec_.gvec_count(r); igloc++) {
716 int ig = gvec_.gvec_offset(r) + igloc;
717 int igsh = gvec_.shell(ig);
718 auto G = gvec_.gvec<index_domain_t::global>(ig);
719 if (spl_num_gsh_.location(igsh).ib == comm_.rank()) {
720 for (int x : {0, 1, 2}) {
721 gvec_remapped_(x, a2a_recv_.offsets[r] + counts[r]) = G[x];
722 }
723 gvec_shell_remapped_(a2a_recv_.offsets[r] + counts[r]) = igsh;
724 counts[r]++;
725 }
726 }
727 }
728 for (int ig = 0; ig < gvec_count_remapped(); ig++) {
729 idx_gvec_[gvec_remapped(ig)] = ig;
730 }
731 /* sanity check */
732 for (int igloc = 0; igloc < this->gvec_count_remapped(); igloc++) {
733 auto G = this->gvec_remapped(igloc);
734 if (this->index_by_gvec(G) != igloc) {
735 RTE_THROW("Wrong remapped index of G-vector");
736 }
737 int igsh = this->gvec_shell_remapped(igloc);
738 if (igsh != this->gvec().shell(G)) {
739 RTE_THROW("Wrong remapped shell of G-vector");
740 }
741 }
742}
743
744void serialize(serializer& s__, Gvec const& gv__)
745{
746 serialize(s__, gv__.vk_);
747 serialize(s__, gv__.Gmax_);
748 serialize(s__, gv__.lattice_vectors_);
749 serialize(s__, gv__.reduce_gvec_);
750 serialize(s__, gv__.bare_gvec_);
751 serialize(s__, gv__.num_gvec_);
752 serialize(s__, gv__.num_gvec_shells_);
753 serialize(s__, gv__.gvec_full_index_);
754 serialize(s__, gv__.gvec_shell_);
755 serialize(s__, gv__.gvec_shell_len_);
756 serialize(s__, gv__.gvec_index_by_xy_);
757 serialize(s__, gv__.z_columns_);
758 serialize(s__, gv__.gvec_distr_);
759 serialize(s__, gv__.zcol_distr_);
760 serialize(s__, gv__.gvec_base_mapping_);
761 serialize(s__, gv__.offset_);
762 serialize(s__, gv__.count_);
763}
764
765void deserialize(serializer& s__, Gvec& gv__)
766{
767 deserialize(s__, gv__.vk_);
768 deserialize(s__, gv__.Gmax_);
769 deserialize(s__, gv__.lattice_vectors_);
770 deserialize(s__, gv__.reduce_gvec_);
771 deserialize(s__, gv__.bare_gvec_);
772 deserialize(s__, gv__.num_gvec_);
773 deserialize(s__, gv__.num_gvec_shells_);
774 deserialize(s__, gv__.gvec_full_index_);
775 deserialize(s__, gv__.gvec_shell_);
776 deserialize(s__, gv__.gvec_shell_len_);
777 deserialize(s__, gv__.gvec_index_by_xy_);
778 deserialize(s__, gv__.z_columns_);
779 deserialize(s__, gv__.gvec_distr_);
780 deserialize(s__, gv__.zcol_distr_);
781 deserialize(s__, gv__.gvec_base_mapping_);
782 deserialize(s__, gv__.offset_);
783 deserialize(s__, gv__.count_);
784}
785
786} // namespace sddk
787
788} // namespace sirius
Helper class to create FFT grids of given sizes and compute indices in space- and frequency domains.
Definition: fft3d_grid.hpp:38
const std::pair< int, int > & limits(int idim__) const
Limits of a given dimension.
Definition: fft3d_grid.hpp:102
int freq_by_coord(int x__) const
Get frequency by coordinate.
Definition: fft3d_grid.hpp:132
mpi::Communicator const & comm_fft() const
Return FFT communicator.
Definition: gvec.hpp:811
int count() const
Local number of G-vectors for FFT-friendly distribution for this rank.
Definition: gvec.hpp:829
mpi::block_data_descriptor gvec_fft_slab_
Distribution of G-vectors inside FFT-friendly "fat" slab.
Definition: gvec.hpp:790
mpi::Communicator const & comm_ortho_fft() const
Return a communicator that is orthogonal to the FFT communicator.
Definition: gvec.hpp:817
mpi::Communicator const & comm_ortho_fft_
Communicator which is orthogonal to FFT communicator.
Definition: gvec.hpp:781
sddk::mdarray< int, 2 > rank_map_
Mapping of MPI ranks used to split G-vectors to a 2D grid.
Definition: gvec.hpp:793
sddk::mdarray< double, 2 > gkvec_cart_array_
Cartesian coordinaes of a local set of G+k-vectors.
Definition: gvec.hpp:800
mpi::Communicator const & comm_fft_
Communicator for the FFT.
Definition: gvec.hpp:778
void update_gkvec_cart()
Update Cartesian coordinates after a change in lattice vectors.
Definition: gvec.hpp:898
sddk::mdarray< int, 2 > gvec_array_
Lattice coordinates of a local set of G-vectors.
Definition: gvec.hpp:797
mpi::block_data_descriptor gvec_distr_fft_
Distribution of G-vectors for FFT.
Definition: gvec.hpp:784
Gvec const & gvec_
Reference to the G-vector instance.
Definition: gvec.hpp:775
Gvec const & gvec() const
Return the original (not reshuffled) G-vector class.
Definition: gvec.hpp:847
void pile_gvec()
Stack together the G-vector slabs to make a larger ("fat") slab for a FFT driver.
Definition: gvec.cpp:599
int num_zcol_local_
Local number of z-columns.
Definition: gvec.hpp:787
A set of G-vectors for FFTs and G+k basis functions.
Definition: gvec.hpp:130
int offset_
Offset in the global index for the local part of G-vectors.
Definition: gvec.hpp:236
int gvec_count(int rank__) const
Number of G-vectors for a fine-grained distribution.
Definition: gvec.hpp:498
int num_zcol() const
Return global number of z-columns.
Definition: gvec.hpp:688
sddk::mdarray< double, 1 > gvec_shell_len_
Radii (or lengths) of G-vector shells in a.u.^-1.
Definition: gvec.hpp:170
mpi::block_data_descriptor zcol_distr_
Fine-grained distribution of z-columns.
Definition: gvec.hpp:193
sddk::mdarray< uint32_t, 1 > gvec_full_index_
Mapping between G-vector index [0:num_gvec_) and a full index.
Definition: gvec.hpp:161
sddk::mdarray< double, 2 > gvec_cart_
Cartiesian coordinaes of a local set of G-vectors.
Definition: gvec.hpp:221
int num_gvec() const
Return the total number of G-vectors within the cutoff.
Definition: gvec.hpp:478
std::vector< double > gvec_shell_len_local_
Radii of G-vector shells in the local index counting [0, num_gvec_shells_local)
Definition: gvec.hpp:179
int offset() const
Offset (in the global index) of G-vectors for a fine-grained distribution for a current MPI rank.
Definition: gvec.hpp:520
sddk::mdarray< int, 2 > gvec_
Lattice coordinates of a local set of G-vectors.
Definition: gvec.hpp:215
std::vector< z_column_descriptor > z_columns_
Global list of non-zero z-columns.
Definition: gvec.hpp:187
sddk::mdarray< int, 1 > gvec_shell_
Index of the shell to which the given G-vector belongs.
Definition: gvec.hpp:164
sddk::mdarray< int, 1 > gvec_base_mapping_
Mapping between current and base G-vector sets.
Definition: gvec.hpp:211
void distribute_z_columns()
Distribute z-columns between MPI ranks.
Definition: gvec.cpp:208
sddk::mdarray< double, 1 > gvec_len_
Length of the local fraction of G-vectors.
Definition: gvec.hpp:227
r3::vector< double > vk_
k-vector of G+k.
Definition: gvec.hpp:133
int zcol_offset(int rank__) const
Offset in the global index of z-columns for a given rank.
Definition: gvec.hpp:491
int count() const
Number of G-vectors for a fine-grained distribution for the current MPI rank.
Definition: gvec.hpp:506
void find_gvec_shells()
Find a list of G-vector shells.
Definition: gvec.cpp:275
int gvec_offset(int rank__) const
Offset (in the global index) of G-vectors for a fine-grained distribution.
Definition: gvec.hpp:512
bool bare_gvec_
True if this a list of G-vectors without k-point shift.
Definition: gvec.hpp:148
int count_
Local number of G-vectors.
Definition: gvec.hpp:239
int num_shells() const
Return number of G-vector shells.
Definition: gvec.hpp:532
double sym_tol_
Symmetry tolerance of the real-space lattice.
Definition: gvec.hpp:245
void init_gvec_cart_local()
Initialize Cartesian coordinates of the local fraction of G-vectors.
Definition: gvec.cpp:388
int num_gvec_shells_local_
Local number of G-vector shells for the local number of G-vectors.
Definition: gvec.hpp:176
int zcol_count(int rank__) const
Number of z-columns for a fine-grained distribution.
Definition: gvec.hpp:484
mpi::block_data_descriptor gvec_distr_
Fine-grained distribution of G-vectors.
Definition: gvec.hpp:190
r3::vector< int > gvec(int ig__) const
Return G vector in fractional coordinates.
Definition: gvec.hpp:540
sddk::mdarray< double, 2 > gkvec_cart_
Cartesian coordinaes of a local set of G+k-vectors.
Definition: gvec.hpp:224
bool reduce_gvec_
Indicates that G-vectors are reduced by inversion symmetry.
Definition: gvec.hpp:145
int shell(int ig__) const
Return index of the G-vector shell by the G-vector index.
Definition: gvec.hpp:608
int index_by_gvec(r3::vector< int > const &G__) const
Return a global G-vector index in the range [0, num_gvec) by the G-vector.
Definition: gvec.cpp:518
int num_zcol_local_
Local number of z-columns.
Definition: gvec.hpp:242
int num_gvec_shells_
Number of G-vector shells (groups of G-vectors with the same length).
Definition: gvec.hpp:167
Gvec const * gvec_base_
Set of G-vectors on which the current G-vector distribution can be based.
Definition: gvec.hpp:198
double Gmax_
Cutoff for |G+k| vectors.
Definition: gvec.hpp:136
double shell_len(int igs__) const
Return length of the G-vector shell.
Definition: gvec.hpp:619
sddk::mdarray< double, 2 > gkvec_
Lattice coordinates of a local set of G+k-vectors.
Definition: gvec.hpp:218
auto const & gvec_local() const
Return local list of G-vectors.
Definition: gvec.hpp:726
r3::vector< int > gvec_by_full_index(uint32_t idx__) const
Return corresponding G-vector for an index in the range [0, num_gvec).
Definition: gvec.cpp:34
std::vector< int > gvec_shell_idx_local_
Mapping between local index of G-vector and local G-shell index.
Definition: gvec.hpp:182
void init(fft::Grid const &fft_grid)
Initialize everything.
Definition: gvec.cpp:420
r3::matrix< double > lattice_vectors_
Reciprocal lattice vectors.
Definition: gvec.hpp:139
int num_gvec_
Total number of G-vectors.
Definition: gvec.hpp:151
void init_gvec_local()
Initialize lattice coordinates of the local fraction of G-vectors.
Definition: gvec.cpp:372
void find_z_columns(double Gmax__, fft::Grid const &fft_box__)
Find z-columns of G-vectors inside a sphere with Gmax radius.
Definition: gvec.cpp:48
MPI communicator wrapper.
int size() const
Size of the communicator (number of ranks).
void allreduce(T *buffer__, int count__) const
Perform the in-place (the output buffer is used as the input buffer) all-to-all reduction.
int rank() const
Rank of MPI process inside communicator.
Simple implementation of 3d vector.
Definition: r3.hpp:45
void zero(memory_t mem__, size_t idx0__, size_t n__)
Zero n elements starting from idx0.
Definition: memory.hpp:1316
size_t size() const
Return total size (number of elements) of the array.
Definition: memory.hpp:1207
Serialize and deserialize objects.
Definition: serializer.hpp:36
Declaration and implementation of Gvec class.
Crystal lattice functions.
auto spherical_coordinates(vector< double > vc)
Transform Cartesian coordinates [x,y,z] to spherical coordinates [r,theta,phi].
Definition: r3.hpp:590
Namespace of the SIRIUS library.
Definition: sirius.f90:5
@ global
Global index.
strong_type< int, struct __block_id_tag > block_id
ID of the block.
Definition: splindex.hpp:108
strong_type< int, struct __n_blocks_tag > n_blocks
Number of blocks to which the global index is split.
Definition: splindex.hpp:105
auto conj(double x__)
Return complex conjugate of a number. For a real value this is the number itself.
Definition: math_tools.hpp:165
Serializer for simple data structures.
Descriptor of the z-column (x,y fixed, z varying) of the G-vectors.
Definition: gvec.hpp:49
int y
Y-coordinate (can be negative and positive).
Definition: gvec.hpp:53
int x
X-coordinate (can be negative and positive).
Definition: gvec.hpp:51