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>
67 std::string groupName)
70 readGroup(groupName, group);
87 unitCellPtr_ = &unitCell;
91 waves_.allocate(nWave_);
92 waveIds_.allocate(nWave_);
103 UTIL_THROW(
"Basis failed validity test suite");
107 isInitialized_ =
true;
122 void Basis<D>::makeWaves()
124 IntVec<D> meshDimensions = mesh().dimensions();
125 std::vector< TWave<D> > twaves;
126 twaves.reserve(nWave_);
132 for (itr.begin(); !itr.atEnd(); ++itr) {
133 w.indicesDft = itr.position();
134 v = shiftToMinimum(w.indicesDft, meshDimensions, *unitCellPtr_);
136 w.sqNorm = unitCell().ksq(v);
141 TWaveNormComp<D> comp;
142 std::sort(twaves.begin(), twaves.end(), comp);
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;
187 std::vector< TWave<D> > tempStar;
195 std::complex<double> coeff;
200 const double epsilon = 1.0E-8;
201 IntVec<D> meshDimensions = mesh().dimensions();
360 Gsq_max = waves_[0].sqNorm;
361 for (i = 1; i <= nWave_; ++i) {
364 bool newList =
false;
367 listSize = listEnd - listBegin;
370 Gsq = waves_[i].sqNorm;
371 if (Gsq > Gsq_max + epsilon) {
374 listSize = listEnd - listBegin;
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;
390 UTIL_CHECK( std::abs(wave.sqNorm-waves_[j].sqNorm)
405 rootItr = list.begin();
413 while (list.size() > 0) {
415 rootVecBz = rootItr->indicesBz;
416 rootVecDft = rootItr->indicesDft;
417 Gsq = rootItr->sqNorm;
423 for (j = 0; j < group.size(); ++j) {
427 vec = rootVecBz*group[j];
430 UTIL_CHECK(std::abs(Gsq - unitCell().ksq(vec)) < epsilon);
434 wave.indicesBz = shiftToMinimum(vec, meshDimensions,
436 wave.indicesDft = vec;
437 mesh().shift(wave.indicesDft);
442 for (k = 0; k < D; ++k) {
443 wave.phase += rootVecBz[k]*(group[j].t(k));
445 while (wave.phase > 0.5) {
448 while (wave.phase <= -0.5) {
459 if (wave.indicesDft == rootVecDft) {
460 if (std::abs(wave.phase) > 1.0E-6) {
467 setItr = star.find(wave);
469 if (setItr == star.end()) {
481 phase_diff = setItr->phase - wave.phase;
482 while (phase_diff > 0.5) {
485 while (phase_diff <= -0.5) {
488 if (std::abs(phase_diff) > 1.0E-6) {
498 setItr = star.begin();
499 for ( ; setItr != star.end(); ++setItr) {
500 tempStar.push_back(*setItr);
505 std::sort(tempStar.begin(), tempStar.end(), waveBzComp);
508 int tempStarSize = tempStar.size();
509 for (j = 0; j < tempStarSize; ++j) {
510 list.erase(tempStar[j]);
511 tempList.
append(tempStar[j]);
519 nBasisWave_ += star.size();
524 newStar.beginId = starBegin;
525 newStar.endId = newStar.beginId + star.size();
526 newStar.size = star.size();
527 newStar.cancel = cancel;
531 if (nextInvert == -1) {
536 newStar.invertFlag = -1;
537 rootItr = list.begin();
546 nVec.negate(rootVecBz);
549 (*meshPtr_).shift(nVec);
552 bool inverseFound =
false;
553 setItr = star.begin();
554 for ( ; setItr != star.end(); ++setItr) {
555 if (nVec == setItr->indicesDft) {
566 newStar.invertFlag = 0;
567 rootItr = list.begin();
575 newStar.invertFlag = 1;
582 setItr = list.begin();
583 for ( ; setItr != list.end(); ++setItr) {
584 if (nVec == setItr->indicesDft) {
596 std::cout <<
"Inverse not found for: " <<
"\n";
597 std::cout <<
" vec (ft):"
598 << rootVecDft <<
"\n";
599 std::cout <<
" vec (bz):"
601 std::cout <<
"-vec (dft):" << nVec <<
"\n";
609 stars_.append(newStar);
610 starBegin = newStar.endId;
621 for (j = 0; j < tempList.
size(); ++j) {
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);
628 if (std::abs(
imag(coeff)) < 1.0E-6) {
629 coeff = std::complex<double>(
real(coeff), 0.0);
631 if (std::abs(
real(coeff)) < 1.0E-6) {
632 coeff = std::complex<double>(0.0,
imag(coeff));
634 waves_[k].coeff = coeff;
648 nStar_ = stars_.size();
665 std::complex<double> rootCoeff;
666 std::complex<double> partCoeff;
667 std::complex<double> d;
669 for (i = 0; i < nStar_; ++i) {
673 if (stars_[i].invertFlag == 0) {
676 for (j = stars_[i].beginId; j < stars_[i].endId; ++j) {
677 if (waves_[j].inverseId < 0) {
680 nVec.negate(waves_[j].indicesBz);
681 (*meshPtr_).shift(nVec);
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;
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;
706 if (waves_[j].inverseId < 0) {
708 std::cout <<
"Inverse not found in closed star"
710 std::cout <<
"G = " << waves_[j].indicesBz
711 <<
", coeff = " << waves_[j].coeff
713 std::cout <<
"All waves in star " << i
715 for (k=stars_[i].beginId; k < stars_[i].endId; ++k) {
716 std::cout << waves_[k].indicesBz <<
" "
717 << waves_[k].coeff << std::endl;
728 rootId = stars_[i].beginId;
729 stars_[i].waveBz = waves_[rootId].indicesBz;
731 if (stars_[i].cancel) {
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;
742 partId = waves_[rootId].inverseId;
745 rootCoeff = waves_[rootId].coeff;
746 for (j = stars_[i].beginId; j < stars_[i].endId; ++j) {
747 waves_[j].coeff /= rootCoeff;
749 rootCoeff = waves_[rootId].coeff;
754 if (partId != rootId) {
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) {
760 if (
real(d) < -1.0E-4) {
763 if (std::abs(
real(d)) <= 1.0E-4) {
768 for (j=stars_[i].beginId; j < stars_[i].endId; ++j){
769 waves_[j].coeff /= d;
779 if (stars_[i].invertFlag == 1) {
784 UTIL_CHECK(stars_[i].size == stars_[i+1].size);
785 UTIL_CHECK(stars_[i].cancel == stars_[i+1].cancel);
789 for (j = stars_[i].beginId; j < stars_[i].endId; ++j) {
790 if (waves_[j].inverseId < 0) {
793 nVec.negate(waves_[j].indicesBz);
794 (*meshPtr_).shift(nVec);
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;
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;
827 for (j = stars_[i+1].beginId; j < stars_[i+1].endId; ++j) {
833 rootId = stars_[i].beginId;
834 stars_[i].waveBz = waves_[rootId].indicesBz;
838 partId = waves_[rootId].inverseId;
839 stars_[i+1].waveBz = waves_[partId].indicesBz;
841 if (stars_[i].cancel) {
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;
847 for (j = stars_[i+1].beginId; j < stars_[i+1].endId; ++j) {
848 waves_[j].coeff = czero;
854 rootCoeff = waves_[rootId].coeff;
855 for (j = stars_[i].beginId; j < stars_[i].endId; ++j) {
856 waves_[j].coeff /= rootCoeff;
860 partCoeff = waves_[partId].coeff;
861 for (j = stars_[i+1].beginId; j < stars_[i+1].endId; ++j) {
862 waves_[j].coeff /= partCoeff;
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;
884 for (i = 0; i < nWave_; ++i) {
885 if (std::abs(
real(waves_[i].coeff)) < 1.0E-8) {
887 = std::complex<double>(0.0,
imag(waves_[i].coeff));
889 if (std::abs(
imag(waves_[i].coeff)) < 1.0E-8) {
891 = std::complex<double>(
real(waves_[i].coeff), 0.0);
896 for (i = 0; i < nWave_; ++i) {
897 vec = waves_[i].indicesDft;
900 for (j = 0; j < D; ++j) {
906 if ((vec[D-1] + 1) > (meshDimensions[D-1]/2 + 1)) {
907 waves_[i].implicit =
true;
909 waves_[i].implicit =
false;
913 waveIds_[mesh().rank(vec)] = i;
917 starIds_.allocate(nBasis_);
919 for (i = 0; i < nStar_; ++i) {
920 stars_[i].starId = i;
921 if (stars_[i].cancel) {
922 stars_[i].basisId = -1;
924 stars_[i].basisId = j;
943 out <<
"N_wave" << std::endl;
945 out <<
" " << nWave_ << std::endl;
947 out <<
" " << nBasisWave_ << std::endl;
951 for (i = 0; i < nWave_; ++i) {
952 starId = waves_[i].starId;
953 if (outputAll || (!stars_[starId].cancel)) {
956 for (j = 0; j < D; ++j) {
957 out <<
Int(waves_[i].indicesBz[j], 5);
959 out <<
Int(waves_[i].starId, 6);
960 out <<
" " <<
Dbl(waves_[i].coeff.real(), 15);
961 out <<
" " <<
Dbl(waves_[i].coeff.imag(), 15);
974 out <<
"N_star" << std::endl
975 <<
" " << nStar_ << std::endl;
977 out <<
"N_basis" << std::endl
978 <<
" " << nBasis_ << std::endl;
983 for (i = 0; i < nStar_; ++i) {
984 if (outputAll || (!stars_[i].cancel)) {
985 out <<
Int(stars_[i].basisId, 6);
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);
992 out <<
Int(stars_[i].cancel, 4);
994 for (j = 0; j < D; ++j) {
995 out <<
Int(stars_[i].waveBz[j], 6);
1007 int is, ib, iw, iwp, j;
1010 if (nWave_ != mesh().size()) {
1011 std::cout <<
"nWave != size of mesh" << std::endl;
1020 if (
wave(iw).indicesDft != v) {
1021 std::cout <<
"Inconsistent waveId and Wave::indicesDft"
1028 for (iw = 0; iw < nWave_; ++iw) {
1031 v = waves_[iw].indicesBz;
1032 Gsq = unitCell().ksq(v);
1033 if (std::abs(Gsq - waves_[iw].sqNorm) > 1.0E-8) {
1035 std::cout <<
"Incorrect sqNorm:" <<
"\n"
1036 <<
"wave.indicesBz = " <<
"\n"
1037 <<
"wave.sqNorm = " << waves_[iw].sqNorm <<
"\n"
1038 <<
"|v|^{2} = " << Gsq <<
"\n";
1044 if (v != waves_[iw].indicesDft) {
1046 std::cout <<
"shift(indicesBz) != indicesDft" << std::endl;
1051 is = waves_[iw].starId;
1052 if (iw < stars_[is].beginId) {
1054 std::cout <<
"Wave::starId < Star::beginId" << std::endl;
1057 if (iw >= stars_[is].endId) {
1059 std::cout <<
"Wave::starId >= Star::endId" << std::endl;
1064 if (waves_[iw].inverseId < 0) {
1066 std::cout <<
"Wave::inverseId not assigned\n";
1067 std::cout <<
"G = " << waves_[iw].indicesBz << std::endl;
1072 v.
negate(waves_[iw].indicesBz);
1074 iwp = waves_[iw].inverseId;
1075 if (waves_[iwp].indicesDft != v) {
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;
1085 if (waves_[iwp].inverseId != iw) {
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;
1094 if (waves_[iw].implicit ==
true && waves_[iwp].implicit ==
true)
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;
1106 for (is = 0; is < nStar_; ++is) {
1109 nWave += stars_[is].size;
1110 if (stars_[is].size != stars_[is].endId - stars_[is].beginId) {
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;
1120 if (stars_[is].beginId != stars_[is-1].endId) {
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;
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) {
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;
1147 v = waves_[stars_[is].beginId].indicesBz;
1148 if (stars_[is].waveBz != v) {
1150 std::cout <<
"waveBz of star != first wave of star"
1152 std::cout <<
"star id " << is << std::endl;
1153 std::cout <<
"waveBz " << stars_[is].waveBz
1155 std::cout <<
"first wave " << v << std::endl;
1161 for (iw = stars_[is].beginId; iw < stars_[is].endId; ++iw) {
1162 if (waves_[iw].starId != is) {
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";
1175 if (stars_[is].starId != is) {
1177 std::cout <<
"stars_[is].starId != is for "
1178 <<
"is = " << is <<
"\n";
1183 ib = stars_[is].basisId;
1184 if (stars_[is].cancel) {
1187 std::cout <<
"basisId != -1 for cancelled star\n";
1188 std::cout <<
"star id = " << is <<
"\n";
1192 if (starIds_[ib] != is) {
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]
1204 for (iw = stars_[is].beginId + 1; iw < stars_[is].endId; ++iw) {
1205 if (waves_[iw].indicesBz > waves_[iw-1].indicesBz) {
1207 std::cout <<
"Failure of ordering by indicesB within star"
1211 if (waves_[iw].indicesBz == waves_[iw-1].indicesBz) {
1213 std::cout <<
"Equal values of indicesBz within star"
1220 if (stars_[is].cancel) {
1221 for (iw = stars_[is].beginId + 1; iw < stars_[is].endId; ++iw)
1223 if (std::abs(waves_[iw].coeff) > 1.0E-8) {
1225 std::cout <<
"Nonzero coefficient in a cancelled star"
1227 std::cout <<
"G = " << waves_[iw].indicesBz
1228 <<
" coeff = " << waves_[iw].coeff
1238 if (stars_[nStar_-1].endId != mesh().size()) {
1240 std::cout <<
"Star endId of last star != mesh size" << std::endl;
1243 if (
nWave != mesh().size()) {
1245 std::cout <<
"Sum of star sizes != mesh size" << std::endl;
1251 std::complex<double> cdel;
1254 while (is < nStar_) {
1255 cancel = stars_[is].cancel;
1257 if (stars_[is].invertFlag == 0) {
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) {
1266 std::cout <<
"Inverse not found in closed star"
1268 std::cout <<
"G = " << waves_[iw].indicesBz
1269 <<
"coeff = " << waves_[iw].coeff
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";
1279 cdel = std::conj(waves_[iwp].coeff);
1280 cdel -= waves_[iw].coeff;
1281 if (std::abs(cdel) > 1.0E-8) {
1283 std::cout <<
"Function for closed star is not real:"
1285 std::cout <<
"+G = " << waves_[iw].indicesBz
1286 <<
" coeff = " << waves_[iw].coeff
1288 std::cout <<
"-G = " << waves_[iwp].indicesBz
1289 <<
" coeff = " << waves_[iwp].coeff
1291 std::cout <<
"Coefficients are not conjugates."
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";
1310 if (stars_[is].invertFlag != 1) {
1312 std::cout <<
"Expected invertFlag == 1" << std::endl;
1315 if (stars_[is+1].invertFlag != -1) {
1317 std::cout <<
"Expected invertFlag == -1" << std::endl;
1320 if (stars_[is+1].size != stars_[is].size) {
1322 std::cout <<
"Partner stars of different size" << std::endl;
1325 if (stars_[is+1].cancel != stars_[is].cancel) {
1327 std::cout <<
"Partners stars with different cancel flags"
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;
1342 for (iw = begin1; iw < end1; ++iw) {
1343 iwp = waves_[iw].inverseId;
1344 if (waves_[iwp].starId != is + 1) {
1346 std::cout <<
"Inverse not found for G in open star"
1348 std::cout <<
"First star id = " << is << std::endl;
1349 std::cout <<
"+G = " << waves_[iw].indicesBz
1350 <<
"coeff = " << waves_[iw].coeff
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";
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";
1367 cdel = std::conj(waves_[iwp].coeff);
1368 cdel -= waves_[iw].coeff;
1369 if (std::abs(cdel) > 1.0E-8) {
1371 std::cout <<
"Error of coefficients in open stars:"
1373 std::cout <<
"First star id = " << is << std::endl;
1374 std::cout <<
"+G = " << waves_[iw].indicesBz
1375 <<
" coeff = " << waves_[iw].coeff
1377 std::cout <<
"-G = " << waves_[iwp].indicesBz
1378 <<
" coeff = " << waves_[iwp].coeff
1380 std::cout <<
"Coefficients are not conjugates."
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";
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";
1407 for (ib = 0; ib < nBasis_; ++ib) {
1409 if (stars_[is].cancel) {
1411 std::cout <<
"Star referred to by starIds_ is cancelled\n";
1414 if (stars_[is].basisId != ib) {
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";
An IntVec<D, T> is a D-component vector of elements of integer type T.
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.
A list of wavevectors that are related by space-group symmetries.
Wave const & wave(int id) const
Get a specific Wave, access by integer index.
int nBasis() const
Total number of nonzero symmetry-adapted basis functions.
int waveId(IntVec< D > vector) const
Get the integer index of a wave, as required by wave(int id).
void makeBasis(Mesh< D > const &mesh, UnitCell< D > const &unitCell, SpaceGroup< D > const &group)
Construct basis for a specific mesh and space group.
Signal< void > & signal()
Get a Signal that is triggered by basis initialization.
void outputWaves(std::ostream &out, bool outputAll=false) const
Print a list of all waves to an output stream.
int nWave() const
Total number of wavevectors.
bool isValid() const
Returns true if this basis is valid, false otherwise.
void outputStars(std::ostream &out, bool outputAll=false) const
Print a list of all stars to an output stream.
Crystallographic space group.
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.
Vec< D, T > & negate(const Vec< D, T > &v)
Return negative of vector v.
static const double Pi
Trigonometric constant Pi.
Wrapper for a double precision number, for formatted ostream output.
An automatically growable array, analogous to a std::vector.
int size() const
Return logical size of this array (i.e., current number of elements).
void clear()
Reset to empty state.
void append(Data const &data)
Append an element to the end of the sequence.
Wrapper for an int, for formatted ostream output.
Notifier (or subject) in the Observer design pattern.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
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.
PSCF package top-level namespace.
Comparator for TWave objects, based on TWave::indicesBz.
Comparator for TWave objects, based on TWave::indicesDft.
Simple wave struct for use within Basis construction.