PSCF v1.3
Basis.tpp
1#ifndef PRDC_BASIS_TPP
2#define PRDC_BASIS_TPP
3
4/*
5* PSCF - Polymer Self-Consistent Field
6*
7* Copyright 2015 - 2025, The Regents of the University of Minnesota
8* Distributed under the terms of the GNU General Public License.
9*/
10
11#include "Basis.h"
12#include "TWave.h"
13#include "groupFile.h"
14#include <prdc/crystal/UnitCell.h>
15#include <prdc/crystal/SpaceGroup.h>
16#include <prdc/crystal/shiftToMinimum.h>
17#include <pscf/mesh/Mesh.h>
18#include <pscf/mesh/MeshIterator.h>
19#include <util/signal/Signal.h>
20#include <util/format/Dbl.h>
21#include <util/math/Constants.h>
22
23#include <algorithm>
24#include <vector>
25#include <set>
26#include <fstream>
27
28namespace Pscf {
29namespace Prdc {
30
31 /*
32 * Constructor.
33 */
34 template <int D>
36 : waves_(),
37 stars_(),
38 waveIds_(),
39 starIds_(),
40 nWave_(0),
41 nBasisWave_(0),
42 nStar_(0),
43 nBasis_(0),
44 signalPtr_(nullptr),
45 unitCellPtr_(0),
46 meshPtr_(0),
47 isInitialized_(false)
48 {
49 signalPtr_ = new Signal<void>();
50 }
51
52 /*
53 * Destructor.
54 */
55 template <int D>
57 {
58 delete signalPtr_;
59 }
60
61 /*
62 * Construct basis for pseudo-spectral scft.
63 */
64 template <int D>
65 void Basis<D>::makeBasis(Mesh<D> const & mesh,
66 UnitCell<D> const & unitCell,
67 std::string groupName)
68 {
69 SpaceGroup<D> group;
70 readGroup(groupName, group);
71 makeBasis(mesh, unitCell, group);
72 }
73
74 /*
75 * Construct a symmetry-adapted basis for pseudo-spectral scft.
76 */
77 template <int D>
78 void Basis<D>::makeBasis(Mesh<D> const & mesh,
79 UnitCell<D> const & unitCell,
80 SpaceGroup<D> const & group)
81 {
82 // Precondition: Check compatibility of mesh with space group
83 group.checkMeshDimensions(mesh.dimensions());
84
85 // Save pointers to mesh and unit cell
86 meshPtr_ = &mesh;
87 unitCellPtr_ = &unitCell;
88
89 // Allocate arrays
90 nWave_ = mesh.size();
91 waves_.allocate(nWave_);
92 waveIds_.allocate(nWave_);
93
94 // Make sorted array of waves
95 makeWaves();
96
97 // Identify stars of waves that are related by symmetry
98 makeStars(group);
99
100 // Apply validity test suite
101 bool valid = isValid();
102 if (!valid) {
103 UTIL_THROW("Basis failed validity test suite");
104 }
105
106 // Mark as initialized
107 isInitialized_ = true;
108
109 // Notify any observers of successful basis initialization
110 signal().notify();
111 }
112
113 /*
114 * Construct ordered list of waves.
115 *
116 * On exit:
117 * - Array waves_ contains list of waves ordered by sqNorm.
118 * - Each wave has indicesDft, indicesBz and sqNorm set.
119 * - Array stars_ is still empty.
120 */
121 template <int D>
122 void Basis<D>::makeWaves()
123 {
124 IntVec<D> meshDimensions = mesh().dimensions();
125 std::vector< TWave<D> > twaves;
126 twaves.reserve(nWave_);
127
128 // Loop over dft mesh to generate all waves, add to twaves
129 TWave<D> w;
130 IntVec<D> v;
131 MeshIterator<D> itr(mesh().dimensions());
132 for (itr.begin(); !itr.atEnd(); ++itr) {
133 w.indicesDft = itr.position();
134 v = shiftToMinimum(w.indicesDft, meshDimensions, *unitCellPtr_);
135 w.indicesBz = v;
136 w.sqNorm = unitCell().ksq(v);
137 twaves.push_back(w);
138 }
139
140 // Sort temporary array twaves
141 TWaveNormComp<D> comp;
142 std::sort(twaves.begin(), twaves.end(), comp);
143
144 // Copy temporary array twaves into member variable waves_
145 for (int i = 0; i < nWave_; ++i) {
146 waves_[i].sqNorm = twaves[i].sqNorm;
147 waves_[i].indicesDft = twaves[i].indicesDft;
148 waves_[i].indicesBz = twaves[i].indicesBz;
149 }
150
151 }
152
153 /*
154 * Complete construction of a basis by grouping presorted waves into
155 * stars and completing initialization of all Wave and Star objects.
156 */
157 template <int D>
158 void Basis<D>::makeStars(SpaceGroup<D> const & group)
159 {
160
161 /*
162 * Conceptual definitions:
163 *
164 * - A "list" is a set of wavevectors of equal magnitude.
165 * - A "star" is a set of wavevectors that are related by symmetry.
166 *
167 * Each list may contain one or more complete stars. Complete lists
168 * are identified as an intermediate step in identification of stars.
169 *
170 * During initial processing, wavevectors are temporarily stored in
171 * TWave<D> objects. The following local containers of TWave<D>
172 * objects are used:
173 *
174 * list - a std::set of waves of equal norm (a "list")
175 * star - a std::set of symmetry-related waves (a "star")
176 * tempStar - a sorted star, sorted by descending indicesBz
177 * tempList - a sorted list, with contiguous sorted stars
178 *
179 * Reasons for choice of some C++ standard lib container types:
180 * std::set for "list", "star" prevents duplicates, allows removal
181 * std::vector for tempStar allows use of std::sort
182 */
183
184 // Local TWave<D> containers and associated iterators
185 std::set< TWave<D>, TWaveDftComp<D> > list;
186 std::set< TWave<D>, TWaveDftComp<D> > star;
187 std::vector< TWave<D> > tempStar;
188 GArray< TWave<D> > tempList;
189 typename std::set< TWave<D>, TWaveDftComp<D> >::iterator rootItr;
190 typename std::set< TWave<D>, TWaveDftComp<D> >::iterator setItr;
191
192 // Local variables
193 TWave<D> wave;
194 Basis<D>::Star newStar;
195 std::complex<double> coeff;
196 double Gsq;
197 double Gsq_max;
198 double phase_diff;
199 const double twoPi = 2.0*Constants::Pi;
200 const double epsilon = 1.0E-8;
201 IntVec<D> meshDimensions = mesh().dimensions();
202 IntVec<D> rootVecBz; // BZ indices for root of this star
203 IntVec<D> rootVecDft; // DFT indices for root of this star
204 IntVec<D> vec; // Indices of temporary wavevector
205 IntVec<D> nVec; // Indices of inverse of a wavevector
206 int listBegin = 0; // id of first wave in this list
207 int listEnd = 0; // (id of last wave in this list) + 1
208 int listSize; // listEnd - listBegin
209 int starBegin = 0; // id of first wave in this star
210 int i, j, k;
211 bool cancel;
212
213 /*
214 * Overview of algorithm:
215 *
216 * Precondition: Wavevectors in the array waves_ are sorted in
217 * nondecreasing order by wavevector norm.
218 *
219 * Loop over index i of array waves_ {
220 *
221 * Search for end of a "list" (i.e., contiguous block of waves
222 * of equal magnitude) by identifying changes in magnitude.
223 * The resulting list has indices [listBegin,listEnd-1].
224 * Set newList true.
225 *
226 * // Each list may contain one or more stars.
227 * // Process the newly identified list to identify stars
228 * If (newList) {
229 *
230 * Copy all waves in the range into std::set list
231 *
232 * Set rootItr to the first wave in the list
233 *
234 * // Loop over stars within the list
235 * while (list.size() > 0) {
236 *
237 * // To generate a star from a root wave rootItr,
238 * // loop over symmetry operations of space group.
239 * For each symmetry operation group[j] {
240 * Compute vec = (rootItr->indicesBz)*group[j]
241 * Set phase = rootItr->indicesBz .dot. (group[j].t)
242 * Check for cancellation of the star, set "cancel" flag
243 * Add wave to std::set<TWave> star if not added before
244 * // Here, use of a std::set simplifies test of uniqueness
245 * }
246 *
247 * Copy all waves from star to std::vector<TWave> tempStar
248 * Sort tempStar by indicesBz, in descending order
249 * // Here, use of a std::vector for tempStar allows sorting
250 *
251 * // Add waves in star to tempList and remove from list
252 * For each wave in tempStar {
253 * Append the wave to GArray<TWave> tempList
254 * Erase the wave from std::set<TWave> list
255 * // Here, use of a std::set for list simplifies erasure
256 * }
257 *
258 * Initialize a Star object named newStar
259 * Assign values to members beginId, endId, size, cancel
260 *
261 * // Assign values of newStar.invertFlag, rootItr, nextInvert
262 * if (nextInvert == -1) {
263 * // This is the second star in pair
264 * newStar.invertFlag = -1;
265 * nextInvert = 1;
266 * Set rootItr to the first wave in remaining list
267 * } else {
268 * Search for inverse of rootItr in this star
269 * if inverse is in this star {
270 * newStar.invertFlag = 0
271 * nextInvert = 1;
272 * Set rootItr to the first wave in remaining list
273 * } else
274 * Search for inverse of rootItr in this remaining list
275 * if the inverse is in the remaining list {
276 * newStar.invertFlag = 1
277 * nextInvert = -1;
278 * set rootItr to inverse of current root
279 * }
280 * }
281 *
282 * Append newStar object to GArray<Star> stars_
283 *
284 * } // end loop over stars in a single list
285 *
286 * // At this point, tempList contains the contents of the
287 * // waves_ array occupying the range [beginId, endId-1],
288 * // grouped by stars, with waves within each star sorted
289 * // by indexBz.
290 *
291 * // Overwrite the block of array waves_ with indices in the
292 * // range [beginId, endId-1] with the contents of tempList.
293 * For each wave in tempList {
294 * Copy a TWave in tempList to a Basis:Wave in waves_
295 * Assign a complex coefficient of unit norm to the Wave
296 * }
297 *
298 * // At this point, coefficients of waves have unit magnitude
299 * // and correct relative phases within each star, but not the
300 * // final absolute phases or magnitude.
301 *
302 * } // finish processing of one list (waves of equal norm)
303 *
304 * } // End initial processing of all waves and stars
305 *
306 * // Set phases of wave coefficients
307 * For each star in array stars_ {
308 * if star is closed under inversion (star.invertFlag == 0) {
309 * Find the inverse of every wave in this star
310 * if star is cancelled {
311 * Set coefficients of all waves to zero
312 * } else {
313 * Set the root to the first wave in the star
314 * For each wave in star:
315 * Divide coeff by the root coefficient
316 * }
317 * if (coeffs of root & inverse are not complex conjugates){
318 * Divide all coeffs by a common phasor chosen to obtain
319 * complex conjugate coefficients for root and inverse
320 * }
321 * }
322 * } else
323 * if (star.invertFlag == 1) {
324 * Find the inverse of every wave in this star and next star
325 * Set root of this star to the 1st wave in this star
326 * Set partner to the inverse of the root of this star
327 * If this star is cancelled {
328 * Set coefficients in this star and next to zero
329 * } else {
330 * For each wave in this star:
331 * Divide coeff by the root coefficient
332 * }
333 * For each wave in the next star:
334 * Divide coeff by the partner coefficient
335 * }
336 * }
337 * }
338 * // Note: If star.invertFlag = -1, do nothing because properties
339 * // of this star were all set when processing its partner.
340 * }
341 *
342 * // For all waves, normalize coefficients and set starId
343 * For each star in array stars_ {
344 * For each wave in this star {
345 * Set Wave.starId
346 * Divide coefficient by sqrt(double(star.size))
347 * }
348 * }
349 *
350 * // For all waves, set implicit member and add to look up table
351 * For each wave in array waves_ {
352 * Set Wave::implicit attribute
353 * Set waveIds_[rank] = i
354 * }
355 */
356
357 // Loop over all waves (initial processing of waves)
358 nBasis_ = 0;
359 nBasisWave_ = 0;
360 Gsq_max = waves_[0].sqNorm;
361 for (i = 1; i <= nWave_; ++i) {
362
363 // Determine if this wave begins a new list
364 bool newList = false;
365 if (i == nWave_) {
366 listEnd = i;
367 listSize = listEnd - listBegin;
368 newList = true;
369 } else {
370 Gsq = waves_[i].sqNorm;
371 if (Gsq > Gsq_max + epsilon) {
372 Gsq_max = Gsq;
373 listEnd = i;
374 listSize = listEnd - listBegin;
375 newList = true;
376 }
377 }
378
379 // Process completed list of wavectors of equal norm
380 if (newList) {
381
382 // Copy waves of equal norm into std::set "list"
383 list.clear();
384 tempList.clear();
385 for (j = listBegin; j < listEnd; ++j) {
386 wave.indicesDft = waves_[j].indicesDft;
387 wave.indicesBz = waves_[j].indicesBz;
388 wave.sqNorm = waves_[j].sqNorm;
389 if (j > listBegin) {
390 UTIL_CHECK( std::abs(wave.sqNorm-waves_[j].sqNorm)
391 < 2.0*epsilon );
392 }
393 list.insert(wave);
394 }
395
396 // On entry to each iteration of the loop over stars,
397 // rootItr and nextInvert are known. The iterator rootItr
398 // points to the wave in the remaining list that will be
399 // used as the root of the next star. The flag nextInvert
400 // is equal to -1 iff the previous star was the first of
401 // a pair that are open under inversion, and is equal
402 // to + 1 otherwise.
403
404 // Initial values for first star in this list
405 rootItr = list.begin();
406 int nextInvert = 1;
407
408 // Loop over stars with a list of waves of equal norm,
409 // removing each star from the list as it is identified.
410 // The root of the next star must have been chosen on
411 // entry to each iteration of this loop.
412
413 while (list.size() > 0) {
414
415 rootVecBz = rootItr->indicesBz;
416 rootVecDft = rootItr->indicesDft;
417 Gsq = rootItr->sqNorm;
418 cancel = false;
419 star.clear();
420
421 // Construct a star from root vector, by applying every
422 // symmetry operation in the group to the root wavevector.
423 for (j = 0; j < group.size(); ++j) {
424
425 // Apply symmetry (i.e., multiply by rotation matrix)
426 // vec = rotated wavevector.
427 vec = rootVecBz*group[j];
428
429 // Check that rotated vector has same norm as root.
430 UTIL_CHECK(std::abs(Gsq - unitCell().ksq(vec)) < epsilon);
431
432 // Initialize TWave object associated with rotated wave
433 wave.sqNorm = Gsq;
434 wave.indicesBz = shiftToMinimum(vec, meshDimensions,
435 *unitCellPtr_);
436 wave.indicesDft = vec;
437 mesh().shift(wave.indicesDft);
438
439 // Compute phase for coeff. of wave in basis function.
440 // Convention -pi < phase <= pi.
441 wave.phase = 0.0;
442 for (k = 0; k < D; ++k) {
443 wave.phase += rootVecBz[k]*(group[j].t(k));
444 }
445 while (wave.phase > 0.5) {
446 wave.phase -= 1.0;
447 }
448 while (wave.phase <= -0.5) {
449 wave.phase += 1.0;
450 }
451 wave.phase *= twoPi;
452
453 // Check for cancellation of star: The star is
454 // cancelled if application of any symmetry operation
455 // in the group to the root vector yields a rotated
456 // vector equivalent to the root vector but with a
457 // nonzero phase, creating a contradiction.
458
459 if (wave.indicesDft == rootVecDft) {
460 if (std::abs(wave.phase) > 1.0E-6) {
461 cancel = true;
462 }
463 }
464
465 // Search for an equivalent wave already in the star.
466 // Note: Equivalent waves have equal DFT indices.
467 setItr = star.find(wave);
468
469 if (setItr == star.end()) {
470
471 // If no equivalent wave is found in the star,
472 // then add this wave to the star
473 star.insert(wave);
474
475 } else {
476
477 // If an equivalent wave is found, check if the
478 // phases are equivalent. If not, the star is
479 // cancelled.
480
481 phase_diff = setItr->phase - wave.phase;
482 while (phase_diff > 0.5) {
483 phase_diff -= 1.0;
484 }
485 while (phase_diff <= -0.5) {
486 phase_diff += 1.0;
487 }
488 if (std::abs(phase_diff) > 1.0E-6) {
489 cancel = true;
490 }
491
492 }
493
494 }
495
496 // Copy all waves from std::set star to std::vector tempStar
497 tempStar.clear();
498 setItr = star.begin();
499 for ( ; setItr != star.end(); ++setItr) {
500 tempStar.push_back(*setItr);
501 }
502
503 // Sort tempStar, in descending order by indicesBz.
504 TWaveBzComp<D> waveBzComp;
505 std::sort(tempStar.begin(), tempStar.end(), waveBzComp);
506
507 // Append contents of tempStar to tempList, erase from list
508 int tempStarSize = tempStar.size();
509 for (j = 0; j < tempStarSize; ++j) {
510 list.erase(tempStar[j]);
511 tempList.append(tempStar[j]);
512 }
513 UTIL_CHECK((int)(tempList.size()+list.size()) == listSize);
514
515 // If this star is not cancelled, increment the number of
516 // basis functions (nBasis_) & waves in basis (nBasisWave_)
517 if (!cancel) {
518 ++nBasis_;
519 nBasisWave_ += star.size();
520 }
521
522 // Initialize a Star object
523 // newStar.eigen = Gsq;
524 newStar.beginId = starBegin;
525 newStar.endId = newStar.beginId + star.size();
526 newStar.size = star.size();
527 newStar.cancel = cancel;
528 // Note: newStar.starInvert is not yet known
529
530 // Determine invertFlag, rootItr and nextInvert
531 if (nextInvert == -1) {
532
533 // If this star is 2nd of a pair related by inversion,
534 // set root for next star to 1st wave of remaining list.
535
536 newStar.invertFlag = -1;
537 rootItr = list.begin();
538 nextInvert = 1;
539
540 } else {
541
542 // If this star is not the 2nd of a pair of partners,
543 // then determine if it is closed under inversion.
544
545 // Compute inverse nVec of root vector in FBZ
546 nVec.negate(rootVecBz);
547
548 // Shift inverse nVec to the DFT mesh
549 (*meshPtr_).shift(nVec);
550
551 // Search for inverse of root vector within this star
552 bool inverseFound = false;
553 setItr = star.begin();
554 for ( ; setItr != star.end(); ++setItr) {
555 if (nVec == setItr->indicesDft) {
556 inverseFound = true;
557 break;
558 }
559 }
560
561 if (inverseFound) {
562
563 // If this star is closed under inversion, the root
564 // of next star is the 1st vector of remaining list.
565
566 newStar.invertFlag = 0;
567 rootItr = list.begin();
568 nextInvert = 1;
569
570 } else {
571
572 // This star is open under inversion, and is the
573 // first star of a pair related by inversion
574
575 newStar.invertFlag = 1;
576 nextInvert = -1;
577
578 // Find inverse of the root of this star in the
579 // remaining list, and use this inverse as the
580 // root of the next star.
581
582 setItr = list.begin();
583 for ( ; setItr != list.end(); ++setItr) {
584 if (nVec == setItr->indicesDft) {
585 inverseFound = true;
586 rootItr = setItr;
587 break;
588 }
589 }
590 // If inverseFound, then rootVecDft = nVec
591
592 // Failure to find the inverse here is an error:
593 // It must be either in this star or remaining list
594
595 if (!inverseFound) {
596 std::cout << "Inverse not found for: " << "\n";
597 std::cout << " vec (ft):"
598 << rootVecDft <<"\n";
599 std::cout << " vec (bz):"
600 << rootVecBz <<"\n";
601 std::cout << "-vec (dft):" << nVec << "\n";
602 UTIL_CHECK(inverseFound);
603 }
604
605 }
606
607 }
608
609 stars_.append(newStar);
610 starBegin = newStar.endId;
611
612 }
613 // End loop over stars within a list.
614
615 UTIL_CHECK(list.size() == 0);
616 UTIL_CHECK(tempList.size() == listEnd - listBegin);
617
618 // Copy tempList into corresponding section of waves_,
619 // overwriting the section of waves_ used to create the list.
620 // Compute a complex coefficient of unit norm for each wave.
621 for (j = 0; j < tempList.size(); ++j) {
622 k = j + listBegin;
623 waves_[k].indicesDft = tempList[j].indicesDft;
624 waves_[k].indicesBz = tempList[j].indicesBz;
625 waves_[k].sqNorm = tempList[j].sqNorm;
626 coeff = std::complex<double>(0.0, tempList[j].phase);
627 coeff = exp(coeff);
628 if (std::abs(imag(coeff)) < 1.0E-6) {
629 coeff = std::complex<double>(real(coeff), 0.0);
630 }
631 if (std::abs(real(coeff)) < 1.0E-6) {
632 coeff = std::complex<double>(0.0, imag(coeff));
633 }
634 waves_[k].coeff = coeff;
635 }
636
637 // Processing of list is now complete.
638 // Here, waves_[k].coeff has unit absolute magnitude, and
639 // correct relative phases for waves within a star, but
640 // the coeff may not be unity for the first or last wave in
641 // the star.
642
643 listBegin = listEnd;
644 }
645 // Finished processing a list of waves of equal norm
646
647 } // End loop over all waves
648 nStar_ = stars_.size();
649 // Complete initial processing of all lists and stars
650
651 /*
652 * Conventions for phases of wave coefficients (imposed below):
653 * - Coefficients of the root of each star and its inverse must
654 * be complex conjugates.
655 * - In a closed star (starInvert = 0), the coefficient of the
656 * root must have a non-negative real part. If the root
657 * coefficient is pure imaginary, the imaginary part must
658 * be negative.
659 * - In a pair of open stars that are related by inversion
660 * symmetry, the coefficients of the first wave of the first
661 * star and its inverse must have real coefficients.
662 */
663
664 // Final processing of phases of of waves in stars:
665 std::complex<double> rootCoeff; // Coefficient of root wave
666 std::complex<double> partCoeff; // Coefficient of partner of root
667 std::complex<double> d;
668 int rootId, partId;
669 for (i = 0; i < nStar_; ++i) {
670
671 // Treat open and closed stars differently
672
673 if (stars_[i].invertFlag == 0) {
674
675 // First, assign inverseId for each wave in star
676 for (j = stars_[i].beginId; j < stars_[i].endId; ++j) {
677 if (waves_[j].inverseId < 0) { // if inverseId is unassigned
678
679 // Compute nVec = inverse of root, shifted to DFT mesh
680 nVec.negate(waves_[j].indicesBz);
681 (*meshPtr_).shift(nVec);
682
683 // Find inverse
684
685 // Check in the position that the inverse is
686 // expected to be located for a typical star
687 k = stars_[i].endId - 1 - (j - stars_[i].beginId);
688 if (nVec == waves_[k].indicesDft) {
689 waves_[j].inverseId = k;
690 waves_[k].inverseId = j;
691 } else {
692 // Inverse not in expected position, search full star
693 // (this usually occurs for stars on the edge of the
694 // Brillouin zone)
695 for (k = j; k < stars_[i].endId; ++k) {
696 if (nVec == waves_[k].indicesDft) {
697 waves_[j].inverseId = k;
698 waves_[k].inverseId = j;
699 break;
700 }
701 }
702 }
703
704 // For invertFlag == 0, failure to find nVec in this
705 // star is a fatal error
706 if (waves_[j].inverseId < 0) {
707 std::cout << "\n";
708 std::cout << "Inverse not found in closed star"
709 << std::endl;
710 std::cout << "G = " << waves_[j].indicesBz
711 << ", coeff = " << waves_[j].coeff
712 << std::endl;
713 std::cout << "All waves in star " << i
714 << std::endl;
715 for (k=stars_[i].beginId; k < stars_[i].endId; ++k) {
716 std::cout << waves_[k].indicesBz << " "
717 << waves_[k].coeff << std::endl;
718 }
719 UTIL_CHECK(waves_[j].inverseId >= 0);
720 }
721
722 }
723 }
724
725 // Identify root of this star (star i)
726 // Set the root to be the first wave in the star
727
728 rootId = stars_[i].beginId;
729 stars_[i].waveBz = waves_[rootId].indicesBz;
730
731 if (stars_[i].cancel) {
732
733 // If the star is cancelled, set all coefficients to zero
734 std::complex<double> czero(0.0, 0.0);
735 for (j = stars_[i].beginId; j < stars_[i].endId; ++j) {
736 waves_[j].coeff = czero;
737 }
738
739 } else { // if not cancelled
740
741 // Set partId to index of the inverse of the root
742 partId = waves_[rootId].inverseId;
743
744 // Divide all coefficients by the root coefficient
745 rootCoeff = waves_[rootId].coeff;
746 for (j = stars_[i].beginId; j < stars_[i].endId; ++j) {
747 waves_[j].coeff /= rootCoeff;
748 }
749 rootCoeff = waves_[rootId].coeff;
750 UTIL_CHECK(std::abs(real(rootCoeff) - 1.0) < 1.0E-9);
751 UTIL_CHECK(std::abs(imag(rootCoeff)) < 1.0E-9);
752
753 // Require coefficients of root and inverse are conjugates
754 if (partId != rootId) {
755
756 partCoeff = waves_[partId].coeff;
757 UTIL_CHECK(std::abs(std::abs(partCoeff) - 1.0) < 1.0E-9);
758 if (std::abs(partCoeff - rootCoeff) > 1.0E-6) {
759 d = sqrt(partCoeff);
760 if (real(d) < -1.0E-4) {
761 d = -d;
762 } else
763 if (std::abs(real(d)) <= 1.0E-4) {
764 if (imag(d) < 0.0) {
765 d = -d;
766 }
767 }
768 for (j=stars_[i].beginId; j < stars_[i].endId; ++j){
769 waves_[j].coeff /= d;
770 }
771 }
772
773 }
774
775 } // end if (cancel) ... else ...
776
777 } // end if (stars_[i].invertFlag == 0)
778 else
779 if (stars_[i].invertFlag == 1) {
780
781 // Process a pair of open stars related by inversion.
782
783 // Preconditions:
784 UTIL_CHECK(stars_[i].size == stars_[i+1].size);
785 UTIL_CHECK(stars_[i].cancel == stars_[i+1].cancel);
786 // UTIL_CHECK(stars_[i+1].invertFlag == -1);
787
788 // First, assign inverseId for each wave in pair of stars
789 for (j = stars_[i].beginId; j < stars_[i].endId; ++j) {
790 if (waves_[j].inverseId < 0) { // if inverseId is unassigned
791
792 // Compute nVec = inverse of root, shifted to DFT mesh
793 nVec.negate(waves_[j].indicesBz);
794 (*meshPtr_).shift(nVec);
795
796 // Find inverse
797
798 // Check in the position that the inverse is expected
799 // to be located for a typical pair of stars
800 k = stars_[i+1].endId - 1 - (j - stars_[i].beginId);
801 if (nVec == waves_[k].indicesDft) {
802 waves_[j].inverseId = k;
803 waves_[k].inverseId = j;
804 } else {
805 // Inverse not in expected position, search full star
806 // (this usually occurs for stars on the edge of the
807 // Brillouin zone)
808 k = stars_[i+1].beginId;
809 for ( ; k < stars_[i+1].endId; ++k) {
810 if (nVec == waves_[k].indicesDft) {
811 waves_[j].inverseId = k;
812 waves_[k].inverseId = j;
813 break;
814 }
815 }
816 }
817
818 // For invertFlag = 1, failure to find nVec in the
819 // next star is a fatal error
820 UTIL_CHECK(waves_[j].inverseId >= 0);
821
822 }
823 }
824
825 // Check that inverseId was assigned for all waves in the
826 // next star
827 for (j = stars_[i+1].beginId; j < stars_[i+1].endId; ++j) {
828 UTIL_CHECK(waves_[j].inverseId >= 0);
829 }
830
831 // Identify root of this star (star i)
832 // Set the root to be the first wave in the star
833 rootId = stars_[i].beginId;
834 stars_[i].waveBz = waves_[rootId].indicesBz;
835
836 // Identify root of the next star (star i+1)
837 // Set the root to be the inverse of the root of star i
838 partId = waves_[rootId].inverseId;
839 stars_[i+1].waveBz = waves_[partId].indicesBz;
840
841 if (stars_[i].cancel) {
842
843 std::complex<double> czero(0.0, 0.0);
844 for (j = stars_[i].beginId; j < stars_[i].endId; ++j) {
845 waves_[j].coeff = czero;
846 }
847 for (j = stars_[i+1].beginId; j < stars_[i+1].endId; ++j) {
848 waves_[j].coeff = czero;
849 }
850
851 } else { // if star is not cancelled
852
853 // Divide all coefficients in this star by root coeff
854 rootCoeff = waves_[rootId].coeff;
855 for (j = stars_[i].beginId; j < stars_[i].endId; ++j) {
856 waves_[j].coeff /= rootCoeff;
857 }
858
859 // Divide coefficients in next star by a partner coeff.
860 partCoeff = waves_[partId].coeff;
861 for (j = stars_[i+1].beginId; j < stars_[i+1].endId; ++j) {
862 waves_[j].coeff /= partCoeff;
863 }
864
865 } // end if (cancel) ... else ...
866
867 } // end if (invertFlag==0) ... else if (invertFlag==1) ...
868
869 // Note: If invertFlag == -1, do nothing and continue
870 // Related stars with invertFlag == 1 and -1 are treated together
871
872 } // end loop over stars
873
874 // For all waves, normalize coefficients and set starId
875 for (i = 0; i < nStar_; ++i) {
876 double snorm = 1.0/sqrt(double(stars_[i].size));
877 for (j = stars_[i].beginId; j < stars_[i].endId; ++j) {
878 waves_[j].coeff *= snorm;
879 waves_[j].starId = i;
880 }
881 }
882
883 // Set tiny real and imaginary parts to zero (due to round-off)
884 for (i = 0; i < nWave_; ++i) {
885 if (std::abs(real(waves_[i].coeff)) < 1.0E-8) {
886 waves_[i].coeff
887 = std::complex<double>(0.0, imag(waves_[i].coeff));
888 }
889 if (std::abs(imag(waves_[i].coeff)) < 1.0E-8) {
890 waves_[i].coeff
891 = std::complex<double>(real(waves_[i].coeff), 0.0);
892 }
893 }
894
895 // For each wave, set implicit attribute and add to look-up table
896 for (i = 0; i < nWave_; ++i) {
897 vec = waves_[i].indicesDft;
898
899 // Validity check - check that vec is in dft mesh
900 for (j = 0; j < D; ++j) {
901 UTIL_CHECK(vec[j] >= 0);
902 UTIL_CHECK(vec[j] < meshDimensions[j]);
903 }
904
905 // Set implicit attribute
906 if ((vec[D-1] + 1) > (meshDimensions[D-1]/2 + 1)) {
907 waves_[i].implicit = true;
908 } else {
909 waves_[i].implicit = false;
910 }
911
912 // Look up table for waves
913 waveIds_[mesh().rank(vec)] = i;
914 }
915
916 // Set Star::starId and Star::basisId, add to look-up table starIds_.
917 starIds_.allocate(nBasis_);
918 j = 0;
919 for (i = 0; i < nStar_; ++i) {
920 stars_[i].starId = i;
921 if (stars_[i].cancel) {
922 stars_[i].basisId = -1;
923 } else {
924 stars_[i].basisId = j;
925 UTIL_CHECK(j < nBasis_);
926 starIds_[j] = i;
927 ++j;
928 }
929 }
930 UTIL_CHECK(j == nBasis_);
931
932 }
933
934 // Return value of nBasis
935 template <int D>
937 { return nBasis_; }
938
939
940 template <int D>
941 void Basis<D>::outputWaves(std::ostream& out, bool outputAll) const
942 {
943 out << "N_wave" << std::endl;
944 if (outputAll) {
945 out << " " << nWave_ << std::endl;
946 } else {
947 out << " " << nBasisWave_ << std::endl;
948 }
949 int i, j, k, starId;
950 k = 0;
951 for (i = 0; i < nWave_; ++i) {
952 starId = waves_[i].starId;
953 if (outputAll || (!stars_[starId].cancel)) {
954 out << Int(k, 8);
955 out << Int(i, 8);
956 for (j = 0; j < D; ++j) {
957 out << Int(waves_[i].indicesBz[j], 5);
958 }
959 out << Int(waves_[i].starId, 6);
960 out << " " << Dbl(waves_[i].coeff.real(), 15);
961 out << " " << Dbl(waves_[i].coeff.imag(), 15);
962 out << std::endl;
963 k++;
964 }
965 }
966 }
967
968
969 template <int D>
970 void Basis<D>::outputStars(std::ostream& out, bool outputAll) const
971 {
972 // Output number of stars in appropriate format
973 if (outputAll) {
974 out << "N_star" << std::endl
975 << " " << nStar_ << std::endl;
976 } else {
977 out << "N_basis" << std::endl
978 << " " << nBasis_ << std::endl;
979 }
980
981 // Loop over stars
982 int i, j;
983 for (i = 0; i < nStar_; ++i) {
984 if (outputAll || (!stars_[i].cancel)) {
985 out << Int(stars_[i].basisId, 6); // basisId
986 out << Int(i, 6); // starId
987 out << Int(stars_[i].size, 5)
988 << Int(stars_[i].beginId, 8)
989 << Int(stars_[i].endId, 8)
990 << Int(stars_[i].invertFlag, 4);
991 if (outputAll) {
992 out << Int(stars_[i].cancel, 4);
993 }
994 for (j = 0; j < D; ++j) {
995 out << Int(stars_[i].waveBz[j], 6);
996 }
997 out << std::endl;
998 }
999 }
1000 }
1001
1002 template <int D>
1004 {
1005 double Gsq;
1006 IntVec<D> v;
1007 int is, ib, iw, iwp, j;
1008
1009 // Check total number of waves == # of grid points
1010 if (nWave_ != mesh().size()) {
1011 std::cout << "nWave != size of mesh" << std::endl;
1012 return false;
1013 }
1014
1015 // Loop over dft mesh to check consistency of waveIds_ and waves_
1016 MeshIterator<D> itr(mesh().dimensions());
1017 for (itr.begin(); !itr.atEnd(); ++itr) {
1018 v = itr.position();
1019 iw = waveId(v);
1020 if (wave(iw).indicesDft != v) {
1021 std::cout << "Inconsistent waveId and Wave::indicesDft"
1022 << std::endl;
1023 return false;
1024 }
1025 }
1026
1027 // Loop over elements of waves_, check consistency of wave data.
1028 for (iw = 0; iw < nWave_; ++iw) {
1029
1030 // Check sqNorm
1031 v = waves_[iw].indicesBz;
1032 Gsq = unitCell().ksq(v);
1033 if (std::abs(Gsq - waves_[iw].sqNorm) > 1.0E-8) {
1034 std::cout << "\n";
1035 std::cout << "Incorrect sqNorm:" << "\n"
1036 << "wave.indicesBz = " << "\n"
1037 << "wave.sqNorm = " << waves_[iw].sqNorm << "\n"
1038 << "|v|^{2} = " << Gsq << "\n";
1039 return false;
1040 }
1041
1042 // Check that wave indicesBz is an image of indicesDft
1043 mesh().shift(v);
1044 if (v != waves_[iw].indicesDft) {
1045 std::cout << "\n";
1046 std::cout << "shift(indicesBz) != indicesDft" << std::endl;
1047 return false;
1048 }
1049
1050 // Compare Wave::starId to Star::beginId and Star::endId
1051 is = waves_[iw].starId;
1052 if (iw < stars_[is].beginId) {
1053 std::cout << "\n";
1054 std::cout << "Wave::starId < Star::beginId" << std::endl;
1055 return false;
1056 }
1057 if (iw >= stars_[is].endId) {
1058 std::cout << "\n";
1059 std::cout << "Wave::starId >= Star::endId" << std::endl;
1060 return false;
1061 }
1062
1063 // Check that inverseId has been assigned
1064 if (waves_[iw].inverseId < 0) {
1065 std::cout << "\n";
1066 std::cout << "Wave::inverseId not assigned\n";
1067 std::cout << "G = " << waves_[iw].indicesBz << std::endl;
1068 return false;
1069 }
1070
1071 // Check that inverseId points to the correct wave
1072 v.negate(waves_[iw].indicesBz);
1073 mesh().shift(v);
1074 iwp = waves_[iw].inverseId;
1075 if (waves_[iwp].indicesDft != v) {
1076 std::cout << "\n";
1077 std::cout << "Wave::inverseId is not inverse" << std::endl;
1078 std::cout << "G = " << waves_[iw].indicesBz << std::endl;
1079 std::cout << "-G (from inverseId) = "
1080 << waves_[iwp].indicesBz << std::endl;
1081 return false;
1082 }
1083
1084 // Check that inverseId of the inverse wave is correct
1085 if (waves_[iwp].inverseId != iw) {
1086 std::cout << "\n";
1087 std::cout << "Wave::inverseId values do not agree\n";
1088 std::cout << "+G = " << waves_[iw].indicesBz << std::endl;
1089 std::cout << "-G = " << waves_[iwp].indicesBz << std::endl;
1090 return false;
1091 }
1092
1093 // Check that either this wave or its inverse is explicit
1094 if (waves_[iw].implicit == true && waves_[iwp].implicit == true)
1095 {
1096 std::cout << "\n";
1097 std::cout << "Wave and its inverse are both implicit";
1098 std::cout << "+G = " << waves_[iw].indicesBz << std::endl;
1099 std::cout << "-G = " << waves_[iwp].indicesBz << std::endl;
1100 return false;
1101 }
1102 }
1103
1104 // Loop over all stars (elements of stars_ array)
1105 int nWave = 0;
1106 for (is = 0; is < nStar_; ++is) {
1107
1108 // Check star size
1109 nWave += stars_[is].size;
1110 if (stars_[is].size != stars_[is].endId - stars_[is].beginId) {
1111 std::cout << "\n";
1112 std::cout << "Inconsistent Star::size:" << std::endl;
1113 std::cout << "Star id " << is << std::endl;
1114 std::cout << "star size " << stars_[is].size << std::endl;
1115 std::cout << "Star begin " << stars_[is].beginId << std::endl;
1116 std::cout << "Star end " << stars_[is].endId << std::endl;
1117 return false;
1118 }
1119 if (is > 0) {
1120 if (stars_[is].beginId != stars_[is-1].endId) {
1121 std::cout << "\n";
1122 std::cout << "Star ranges not consecutive:" << std::endl;
1123 std::cout << "Star id " << is << std::endl;
1124 std::cout << "stars_[" << is << "]" << ".beginId = "
1125 << stars_[is].beginId << std::endl;
1126 std::cout << "stars_[" << is-1 << "]" << ".endId = "
1127 << stars_[is-1].endId << std::endl;
1128 return false;
1129 }
1130 }
1131
1132 // Check waveBz indices of star
1133 if (stars_[is].invertFlag == -1) {
1134 v.negate(stars_[is-1].waveBz);
1135 v = shiftToMinimum(v, mesh().dimensions(), *unitCellPtr_);
1136 if (stars_[is].waveBz != v) {
1137 std::cout << "\n";
1138 std::cout << "waveBz of star is not inverse of waveBz "
1139 << "of previous star" << std::endl;
1140 std::cout << "star id " << is << std::endl;
1141 std::cout << "waveBz " << stars_[is].waveBz << std::endl;
1142 std::cout << "waveBz (previous star) "
1143 << stars_[is-1].waveBz << std::endl;
1144 return false;
1145 }
1146 } else {
1147 v = waves_[stars_[is].beginId].indicesBz;
1148 if (stars_[is].waveBz != v) {
1149 std::cout << "\n";
1150 std::cout << "waveBz of star != first wave of star"
1151 << std::endl;
1152 std::cout << "star id " << is << std::endl;
1153 std::cout << "waveBz " << stars_[is].waveBz
1154 << std::endl;
1155 std::cout << "first wave " << v << std::endl;
1156 return false;
1157 }
1158 }
1159
1160 // Check star ids of waves in star
1161 for (iw = stars_[is].beginId; iw < stars_[is].endId; ++iw) {
1162 if (waves_[iw].starId != is) {
1163 std::cout << "\n";
1164 std::cout << "Inconsistent Wave::starId :" << std::endl;
1165 std::cout << "star id " << is << std::endl;
1166 std::cout << "star beginId " << stars_[is].beginId << "\n";
1167 std::cout << "star endId " << stars_[is].endId << "\n";
1168 std::cout << "wave id " << iw << "\n";
1169 std::cout << "wave starId " << waves_[iw].starId << "\n";
1170 return false;
1171 }
1172 }
1173
1174 // Check Star::starId is equal to array index
1175 if (stars_[is].starId != is) {
1176 std::cout << "\n";
1177 std::cout << "stars_[is].starId != is for "
1178 << "is = " << is << "\n";
1179 return false;
1180 }
1181
1182 // Check Star::basisId and starIds_ look up table
1183 ib = stars_[is].basisId;
1184 if (stars_[is].cancel) {
1185 if (ib != -1) {
1186 std::cout << "\n";
1187 std::cout << "basisId != -1 for cancelled star\n";
1188 std::cout << "star id = " << is << "\n";
1189 return false;
1190 }
1191 } else {
1192 if (starIds_[ib] != is) {
1193 std::cout << "\n";
1194 std::cout << "starIds_[stars_[is].basisId] != is for: \n";
1195 std::cout << "is = " << is << "\n";
1196 std::cout << "ib = stars_[is].basisId = " << ib << "\n";
1197 std::cout << "starIds_[ib] = " << starIds_[ib]
1198 << "\n";
1199 return false;
1200 }
1201 }
1202
1203 // Check ordering of waves in star
1204 for (iw = stars_[is].beginId + 1; iw < stars_[is].endId; ++iw) {
1205 if (waves_[iw].indicesBz > waves_[iw-1].indicesBz) {
1206 std::cout << "\n";
1207 std::cout << "Failure of ordering by indicesB within star"
1208 << std::endl;
1209 return false;
1210 }
1211 if (waves_[iw].indicesBz == waves_[iw-1].indicesBz) {
1212 std::cout << "\n";
1213 std::cout << "Equal values of indicesBz within star"
1214 << std::endl;
1215 return false;
1216 }
1217 }
1218
1219 // Check that all coefficients are zero if star is cancelled
1220 if (stars_[is].cancel) {
1221 for (iw = stars_[is].beginId + 1; iw < stars_[is].endId; ++iw)
1222 {
1223 if (std::abs(waves_[iw].coeff) > 1.0E-8) {
1224 std::cout << "\n";
1225 std::cout << "Nonzero coefficient in a cancelled star"
1226 << "\n";
1227 std::cout << "G = " << waves_[iw].indicesBz
1228 << " coeff = " << waves_[iw].coeff
1229 << "\n";
1230 return false;
1231 }
1232 }
1233 }
1234
1235 } // End do loop over all stars
1236
1237 // Check that all waves in mesh are accounted for in stars
1238 if (stars_[nStar_-1].endId != mesh().size()) {
1239 std::cout << "\n";
1240 std::cout << "Star endId of last star != mesh size" << std::endl;
1241 return false;
1242 }
1243 if (nWave != mesh().size()) {
1244 std::cout << "\n";
1245 std::cout << "Sum of star sizes != mesh size" << std::endl;
1246 return false;
1247 }
1248
1249 // Loop over closed stars and related pairs of stars.
1250 // Test closure under inversion and conjugacy of coefficients.
1251 std::complex<double> cdel;
1252 bool cancel;
1253 is = 0;
1254 while (is < nStar_) {
1255 cancel = stars_[is].cancel;
1256
1257 if (stars_[is].invertFlag == 0) {
1258
1259 // Test that star is closed under inversion and real
1260 int begin = stars_[is].beginId;
1261 int end = stars_[is].endId;
1262 for (iw = begin; iw < end; ++iw) {
1263 iwp = waves_[iw].inverseId;
1264 if (waves_[iwp].starId != is) {
1265 std::cout << "\n";
1266 std::cout << "Inverse not found in closed star"
1267 << std::endl;
1268 std::cout << "G = " << waves_[iw].indicesBz
1269 << "coeff = " << waves_[iw].coeff
1270 << std::endl;
1271 std::cout << "All waves in star " << is << "\n";
1272 for (j=begin; j < end; ++j) {
1273 std::cout << waves_[j].indicesBz << " "
1274 << waves_[j].coeff << "\n";
1275 }
1276 return false;
1277 }
1278 if (!cancel) {
1279 cdel = std::conj(waves_[iwp].coeff);
1280 cdel -= waves_[iw].coeff;
1281 if (std::abs(cdel) > 1.0E-8) {
1282 std::cout << "\n";
1283 std::cout << "Function for closed star is not real:"
1284 << "\n";
1285 std::cout << "+G = " << waves_[iw].indicesBz
1286 << " coeff = " << waves_[iw].coeff
1287 << "\n";
1288 std::cout << "-G = " << waves_[iwp].indicesBz
1289 << " coeff = " << waves_[iwp].coeff
1290 << "\n";
1291 std::cout << "Coefficients are not conjugates."
1292 << "\n";
1293 std::cout << "All waves in star " << is << "\n";
1294 for (j=begin; j < end; ++j) {
1295 std::cout << waves_[j].indicesBz << " "
1296 << waves_[j].coeff << "\n";
1297 }
1298 return false;
1299 }
1300 }
1301 }
1302
1303 // Finished processing a closed star, increment counter is
1304 ++is;
1305
1306 } else {
1307
1308 // Test pairs of open stars
1309
1310 if (stars_[is].invertFlag != 1) {
1311 std::cout << "\n";
1312 std::cout << "Expected invertFlag == 1" << std::endl;
1313 return false;
1314 }
1315 if (stars_[is+1].invertFlag != -1) {
1316 std::cout << "\n";
1317 std::cout << "Expected invertFlag == -1" << std::endl;
1318 return false;
1319 }
1320 if (stars_[is+1].size != stars_[is].size) {
1321 std::cout << "\n";
1322 std::cout << "Partner stars of different size" << std::endl;
1323 return false;
1324 }
1325 if (stars_[is+1].cancel != stars_[is].cancel) {
1326 std::cout << "\n";
1327 std::cout << "Partners stars with different cancel flags"
1328 << std::endl;
1329 return false;
1330 }
1331
1332 // Begin and end wave ids for the first and second stars
1333 int begin1 = stars_[is].beginId;
1334 int end1 = stars_[is].endId;
1335 int begin2 = stars_[is+1].beginId;
1336 int end2 = stars_[is+1].endId;
1337
1338 // Check that inverse is in next star and check for
1339 // conjugate coefficients
1340
1341 // Loop over waves in first star
1342 for (iw = begin1; iw < end1; ++iw) {
1343 iwp = waves_[iw].inverseId;
1344 if (waves_[iwp].starId != is + 1) {
1345 std::cout << "\n";
1346 std::cout << "Inverse not found for G in open star"
1347 << std::endl;
1348 std::cout << "First star id = " << is << std::endl;
1349 std::cout << "+G = " << waves_[iw].indicesBz
1350 << "coeff = " << waves_[iw].coeff
1351 << std::endl;
1352 std::cout << "Waves in star " << is
1353 << " (starInvert ==1):" << "\n";
1354 for (j = begin1; j < end1; ++j) {
1355 std::cout << waves_[j].indicesBz << " "
1356 << waves_[j].coeff << "\n";
1357 }
1358 std::cout << "Waves in star " << is+1
1359 << " (starInvert == -1):" << "\n";
1360 for (j=begin2; j < end2; ++j) {
1361 std::cout << waves_[j].indicesBz << " "
1362 << waves_[j].coeff << "\n";
1363 }
1364 return false;
1365 }
1366 if (!cancel) {
1367 cdel = std::conj(waves_[iwp].coeff);
1368 cdel -= waves_[iw].coeff;
1369 if (std::abs(cdel) > 1.0E-8) {
1370 std::cout << "\n";
1371 std::cout << "Error of coefficients in open stars:"
1372 << "\n";
1373 std::cout << "First star id = " << is << std::endl;
1374 std::cout << "+G = " << waves_[iw].indicesBz
1375 << " coeff = " << waves_[iw].coeff
1376 << "\n";
1377 std::cout << "-G = " << waves_[iwp].indicesBz
1378 << " coeff = " << waves_[iwp].coeff
1379 << "\n";
1380 std::cout << "Coefficients are not conjugates."
1381 << "\n";
1382 std::cout << "Waves in star " << is
1383 << " (starInvert ==1):" << "\n";
1384 for (j = begin1; j < end1; ++j) {
1385 std::cout << waves_[j].indicesBz << " "
1386 << waves_[j].coeff << "\n";
1387 }
1388 std::cout << "Waves in star " << is+1
1389 << " (starInvert == -1):" << "\n";
1390 for (j=begin2; j < end2; ++j) {
1391 std::cout << waves_[j].indicesBz << " "
1392 << waves_[j].coeff << "\n";
1393 }
1394 return false;
1395 }
1396 }
1397 }
1398
1399 // Finished processing a pair, increment star counter by 2
1400 is += 2;
1401
1402 } // end if (stars_[is].invertFlag == 0) ... else ...
1403
1404 } // end while (is < nStar_) loop over stars
1405
1406 // Loop over basis functions
1407 for (ib = 0; ib < nBasis_; ++ib) {
1408 is = starIds_[ib];
1409 if (stars_[is].cancel) {
1410 std::cout << "\n";
1411 std::cout << "Star referred to by starIds_ is cancelled\n";
1412 return false;
1413 }
1414 if (stars_[is].basisId != ib) {
1415 std::cout << "\n";
1416 std::cout << "Error: stars_[starIds_[ib]].basisId != ib\n";
1417 std::cout << "Basis function index ib = " << ib << "\n";
1418 std::cout << "is = starIds_[ib] = " << is << "\n";
1419 std::cout << "stars_[is].basisId = "
1420 << stars_[is].basisId << "\n";
1421 return false;
1422 }
1423 }
1424
1425 // The end of this function is reached iff all tests passed.
1426 return true;
1427 }
1428
1429 /*
1430 * Get the associated signal (notifies observers of basis initialization).
1431 */
1432 template <int D>
1434 {
1435 UTIL_CHECK(signalPtr_);
1436 return *signalPtr_;
1437 }
1438
1439}
1440}
1441#endif
An IntVec<D, T> is a D-component vector of elements of integer type T.
Definition IntVec.h:27
Iterator over points in a Mesh<D>.
void begin()
Set iterator to the first point in the mesh.
bool atEnd() const
Is this the end (i.e., one past the last point)?
IntVec< D > position() const
Get current position in the grid, as integer vector.
Description of a regular grid of points in a periodic domain.
Definition Mesh.h:61
A list of wavevectors that are related by space-group symmetries.
Definition Basis.h:482
Wave const & wave(int id) const
Get a specific Wave, access by integer index.
Definition Basis.h:876
int nBasis() const
Total number of nonzero symmetry-adapted basis functions.
Definition Basis.tpp:936
int waveId(IntVec< D > vector) const
Get the integer index of a wave, as required by wave(int id).
Definition Basis.h:890
Basis()
Constructor.
Definition Basis.tpp:35
~Basis()
Destructor.
Definition Basis.tpp:56
void makeBasis(Mesh< D > const &mesh, UnitCell< D > const &unitCell, SpaceGroup< D > const &group)
Construct basis for a specific mesh and space group.
Definition Basis.tpp:78
Signal< void > & signal()
Get a Signal that is triggered by basis initialization.
Definition Basis.tpp:1433
void outputWaves(std::ostream &out, bool outputAll=false) const
Print a list of all waves to an output stream.
Definition Basis.tpp:941
int nWave() const
Total number of wavevectors.
Definition Basis.h:863
bool isValid() const
Returns true if this basis is valid, false otherwise.
Definition Basis.tpp:1003
void outputStars(std::ostream &out, bool outputAll=false) const
Print a list of all stars to an output stream.
Definition Basis.tpp:970
Crystallographic space group.
Definition SpaceGroup.h:32
void checkMeshDimensions(IntVec< D > const &dimensions) const
Check if input mesh dimensions are compatible with space group.
Base template for UnitCell<D> classes, D=1, 2 or 3.
Definition UnitCell.h:56
Vec< D, T > & negate(const Vec< D, T > &v)
Return negative of vector v.
Definition Vec.h:520
static const double Pi
Trigonometric constant Pi.
Definition Constants.h:35
Wrapper for a double precision number, for formatted ostream output.
Definition Dbl.h:40
An automatically growable array, analogous to a std::vector.
Definition GArray.h:34
int size() const
Return logical size of this array (i.e., current number of elements).
Definition GArray.h:455
void clear()
Reset to empty state.
Definition GArray.h:299
void append(Data const &data)
Append an element to the end of the sequence.
Definition GArray.h:306
Wrapper for an int, for formatted ostream output.
Definition Int.h:37
Notifier (or subject) in the Observer design pattern.
Definition Signal.h:39
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
Definition global.h:49
RT imag(CT const &a)
Return the imaginary part of a complex number.
RT real(CT const &a)
Return the real part of a complex number.
Periodic fields and crystallography.
Definition CField.cpp:11
PSCF package top-level namespace.
Definition param_pc.dox:1
Comparator for TWave objects, based on TWave::indicesBz.
Definition TWave.h:76
Comparator for TWave objects, based on TWave::indicesDft.
Definition TWave.h:58
Simple wave struct for use within Basis construction.
Definition TWave.h:25