1#ifndef PRDC_FIELD_IO_UTIL_TPP
2#define PRDC_FIELD_IO_UTIL_TPP
11#include <prdc/crystal/Basis.h>
12#include <prdc/crystal/UnitCell.h>
13#include <prdc/crystal/shiftToMinimum.h>
14#include <prdc/crystal/fieldHeader.h>
16#include <pscf/mesh/Mesh.h>
17#include <pscf/mesh/MeshIterator.h>
18#include <pscf/mesh/MeshIteratorFortran.h>
19#include <pscf/math/IntVec.h>
20#include <pscf/math/complex.h>
22#include <util/containers/DArray.h>
23#include <util/misc/FileMaster.h>
24#include <util/misc/Log.h>
25#include <util/format/Int.h>
26#include <util/format/Dbl.h>
36 template <
int D,
class FT>
40 if (field.isAllocated()) {
41 UTIL_CHECK(field.meshDimensions() == dimensions);
43 field.allocate(dimensions);
50 template <
int D,
class FT>
56 int nMonomerFields = fields.
capacity();
59 for (
int i = 0; i < nMonomer; ++i) {
60 UTIL_CHECK(fields[i].meshDimensions() == dimensions);
64 for (
int i = 0; i < nMonomer; ++i) {
73 template <
int D,
class FT>
82 dimensions = fields[0].meshDimensions();
83 for (
int i = 0; i < nMonomer; ++i) {
84 UTIL_CHECK(fields[i].meshDimensions() == dimensions);
97 int nMonomerArrays = arrays.
capacity();
100 for (
int i = 0; i < nMonomer; ++i) {
105 for (
int i = 0; i < nMonomer; ++i) {
125 for (
int i = 0; i < nMonomer; ++i) {
139 if (label !=
"mesh" && label !=
"ngrid") {
140 std::string msg =
"\n";
141 msg +=
"Error reading field file:\n";
142 msg +=
"Expected mesh or ngrid, but found [";
148 in >> meshDimensionsIn;
149 if (meshDimensionsIn != meshDimensions) {
151 <<
"Inconsistent mesh dimensions:\n"
152 <<
"meshDimensions (expected) = " << meshDimensions <<
"\n"
153 <<
"meshDimensions (from file) = " << meshDimensionsIn <<
"\n";
154 UTIL_THROW(
"Unexpected mesh dimensions in field file header");
162 out <<
"mesh " << std::endl
163 <<
" " << meshDimensions << std::endl;
166 template <
int D,
class ART>
177 for (iter.begin(); !iter.atEnd(); ++iter) {
179 for (
int k = 0; k < nMonomer; ++k) {
180 in >> fields[k][rank];
185 template <
int D,
class ART>
188 IntVec<D>
const& dimensions)
190 MeshIteratorFortran<D> iter(dimensions);
192 for (iter.begin(); !iter.atEnd(); ++iter) {
198 template <
int D,
class ART>
202 IntVec<D>
const& dimensions)
207 MeshIteratorFortran<D> iter(dimensions);
209 for (iter.begin(); !iter.atEnd(); ++iter) {
211 for (j = 0; j < nMonomer; ++j) {
212 out <<
" " <<
Dbl(fields[j][rank], 21, 13);
219 template <
int D,
class ART>
222 IntVec<D>
const& dimensions)
224 MeshIteratorFortran<D> iter(dimensions);
226 for (iter.begin(); !iter.atEnd(); ++iter) {
228 out <<
" " <<
Dbl(field[rank], 21, 13);
235 template <
int D,
class ACT>
239 IntVec<D>
const& dftDimensions)
241 typedef typename ACT::Complex CT;
242 typedef typename ACT::Real RT;
245 MeshIterator<D> iter(dftDimensions);
246 int rank, i, j, idum;
248 for (iter.begin(); !iter.atEnd(); ++iter) {
253 for (j = 0; j < nMonomer; ++j) {
262 template <
int D,
class ACT>
265 IntVec<D>
const& dftDimensions)
267 typedef typename ACT::Complex CT;
268 typedef typename ACT::Real RT;
271 MeshIterator<D> iter(dftDimensions);
274 for (iter.begin(); !iter.atEnd(); ++iter) {
288 template <
int D,
class ACT>
292 IntVec<D>
const& dftDimensions)
297 typedef typename ACT::Complex CT;
298 typedef typename ACT::Real RT;
301 MeshIterator<D> iter(dftDimensions);
304 for (iter.begin(); !iter.atEnd(); ++iter) {
308 for (
int j = 0; j < nMonomer; ++j) {
321 template <
int D,
class ACT>
324 IntVec<D>
const& dftDimensions)
326 typedef typename ACT::Complex CT;
327 typedef typename ACT::Real RT;
330 MeshIterator<D> iter(dftDimensions);
333 for (iter.begin(); !iter.atEnd(); ++iter) {
358 if (label !=
"N_star" && label !=
"N_basis") {
359 std::string msg =
"\n";
360 msg +=
"Error reading field file:\n";
361 msg +=
"Expected N_basis or N_star, but found [";
380 out <<
"N_basis " << std::endl
381 <<
" " << nBasis << std::endl;
402 int fieldCapacity = fields[0].
capacity();
404 for (
int i = 0; i < nMonomer; ++i) {
405 UTIL_CHECK(fields[i].capacity() == fieldCapacity);
406 for (
int j = 0; j < fieldCapacity; ++j) {
412 int nBasis = nBasisIn;
413 if (fieldCapacity < nBasisIn) {
414 nBasis = fieldCapacity;
427 int basisId, basisId2;
430 std::complex<double> coeff, phasor;
432 int nReversedPair = 0;
433 bool waveExists, sizeMatches;
440 for (
int j = 0; j < nMonomer; ++j) {
451 waveBz = shiftToMinimum(waveIn, mesh.
dimensions(), unitCell);
452 waveExists = (waveIn == waveBz);
459 waveId = basis.
waveId(waveDft);
461 starPtr = &basis.
star(starId);
465 if (starPtr->
size == sizeIn) {
469 <<
"Warning: Inconsistent star size (line ignored)\n"
470 <<
"wave from file = " << waveIn <<
"\n"
471 <<
"size from file = " << sizeIn <<
"\n"
472 <<
"size of star = " << starPtr->
size
479 if (waveExists && sizeMatches) {
483 if (starPtr->
waveBz == waveIn) {
486 for (
int j = 0; j < nMonomer; ++j) {
487 fields[j][basisId] = temp[j];
493 <<
"Inconsistent wave of closed star on input\n"
494 <<
"wave from file = " << waveIn <<
"\n"
495 <<
"starId of wave = " << starId <<
"\n"
496 <<
"waveBz of star = " << starPtr->
waveBz
504 for (
int j = 0; j < nMonomer; ++j) {
513 shiftToMinimum(waveIn2, mesh.
dimensions(), unitCell);
519 waveId2 = basis.
waveId(waveDft);
521 starPtr2 = &basis.
star(starId2);
538 for (
int j = 0; j < nMonomer; ++j) {
539 fields[j][basisId] = temp[j];
540 fields[j][basisId2] = temp2[j];
557 shiftToMinimum(nVec, mesh.
dimensions(), unitCell);
606 phasor = phasor/std::abs(phasor);
607 for (
int j = 0; j < nMonomer; ++j) {
608 coeff = std::complex<double>(temp[j],-temp2[j]);
610 fields[j][basisId2] =
real(coeff);
611 fields[j][basisId ] =
imag(coeff);
627 if (nReversedPair > 0) {
629 Log::file() << nReversedPair <<
" reversed pairs of open stars"
630 <<
" detected in readFieldsBasis\n";
646 int fieldCapacity = fields[0].
capacity();
647 for (
int j = 0; j < nMonomer; ++j) {
648 UTIL_CHECK(fieldCapacity == fields[j].capacity());
654 int nStar = basis.
nStar();
655 for (
int i = 0; i < nStar; ++i) {
656 if (ib >= fieldCapacity)
break;
658 for (
int j = 0; j < nMonomer; ++j) {
659 out <<
Dbl(fields[j][ib], 20, 10);
662 for (
int j = 0; j < D; ++j) {
672 template <
int D,
class ACT>
675 Basis<D>
const& basis,
680 typedef typename ACT::Complex CT;
681 typedef typename ACT::Real RT;
684 Mesh<D> dftMesh(dftDimensions);
686 typename Basis<D>::Star
const* starPtr;
687 typename Basis<D>::Wave
const* wavePtr;
688 std::complex<double> component;
689 std::complex<double> coeff;
697 for (rank = 0; rank < dftMesh.size(); ++rank) {
705 while (is < basis.nStar()) {
706 starPtr = &(basis.star(is));
708 if (starPtr->cancel) {
714 ib = starPtr->basisId;
716 if (starPtr->invertFlag == 0) {
719 component = std::complex<double>(in[ib], 0.0);
722 for (iw = starPtr->beginId; iw < starPtr->endId; ++iw) {
723 wavePtr = &basis.wave(iw);
724 if (!wavePtr->implicit) {
725 coeff = component*(wavePtr->coeff);
726 indices = wavePtr->indicesDft;
727 rank = dftMesh.rank(indices);
736 if (starPtr->invertFlag == 1) {
739 component = std::complex<double>(in[ib], -in[ib+1]);
740 component /= sqrt(2.0);
741 starPtr = &(basis.star(is));
742 for (iw = starPtr->beginId; iw < starPtr->endId; ++iw) {
743 wavePtr = &basis.wave(iw);
744 if (!(wavePtr->implicit)) {
745 coeff = component*(wavePtr->coeff);
746 indices = wavePtr->indicesDft;
747 rank = dftMesh.rank(indices);
755 starPtr = &(basis.star(is+1));
757 component = std::complex<double>(in[ib], +in[ib+1]);
758 component /= sqrt(2.0);
759 for (iw = starPtr->beginId; iw < starPtr->endId; ++iw) {
760 wavePtr = &basis.wave(iw);
761 if (!(wavePtr->implicit)) {
762 coeff = component*(wavePtr->coeff);
763 indices = wavePtr->indicesDft;
764 rank = dftMesh.rank(indices);
784 template <
int D,
class ACT>
787 Basis<D>
const& basis,
788 IntVec<D>
const& dftDimensions,
794 typedef typename ACT::Complex CT;
795 typedef typename ACT::Real RT;
800 symmetric =
hasSymmetry(in, basis, dftDimensions, epsilon,
true);
803 <<
"WARNING: Conversion of asymmetric field to"
804 <<
"symmetry-adapted basis in Prdc::convertKgridToBasis."
805 << std::endl << std::endl;
810 Mesh<D> dftMesh(dftDimensions);
812 typename Basis<D>::Star
const* starPtr;
813 typename Basis<D>::Wave
const* wavePtr;
814 std::complex<double> component;
821 for (is = 0; is < basis.nBasis(); ++is) {
827 while (is < basis.nStar()) {
828 starPtr = &(basis.star(is));
830 if (starPtr->cancel) {
836 ib = starPtr->basisId;
838 if (starPtr->invertFlag == 0) {
841 int beginId = starPtr->beginId;
842 int endId = starPtr->endId;
844 bool isImplicit =
true;
846 wavePtr = &basis.wave(beginId + iw);
847 if (!wavePtr->implicit) {
851 wavePtr = &basis.wave(endId - 1 - iw);
852 if (!wavePtr->implicit) {
861 rank = dftMesh.rank(wavePtr->indicesDft);
864 component /= wavePtr->coeff;
865 out[ib] = component.real();
869 if (starPtr->invertFlag == 1) {
874 wavePtr = &basis.wave(starPtr->beginId);
876 if (wavePtr->implicit) {
877 iw = wavePtr->inverseId;
878 starPtr = &(basis.star(is+1));
880 wavePtr = &basis.wave(iw);
884 rank = dftMesh.rank(wavePtr->indicesDft);
887 UTIL_CHECK(std::abs(wavePtr->coeff) > 1.0E-8);
888 component /= wavePtr->coeff;
889 component *= sqrt(2.0);
892 if (starPtr->invertFlag == 1) {
893 out[ib] = component.real();
894 out[ib+1] = -component.imag();
896 out[ib] = component.real();
897 out[ib+1] = component.imag();
913 template <
int D,
class ACT>
915 Basis<D>
const& basis,
916 IntVec<D>
const& dftDimensions,
922 typedef typename ACT::Complex CT;
923 typedef typename ACT::Real RT;
925 typename Basis<D>::Star
const* starPtr;
926 typename Basis<D>::Wave
const* wavePtr;
927 std::complex<double> waveCoeff;
928 std::complex<double> rootCoeff;
929 std::complex<double> diff;
935 double cancelledError(0.0);
936 double uncancelledError(0.0);
939 Mesh<D> dftMesh(dftDimensions);
942 for (is = 0; is < basis.nStar(); ++is) {
943 starPtr = &(basis.star(is));
945 if (starPtr->cancel) {
948 beginId = starPtr->beginId;
949 endId = starPtr->endId;
950 for (iw = beginId; iw < endId; ++iw) {
951 wavePtr = &basis.wave(iw);
952 if (!wavePtr->implicit) {
953 rank = dftMesh.rank(wavePtr->indicesDft);
957 if (std::abs(waveCoeff) > cancelledError) {
958 cancelledError = std::abs(waveCoeff);
959 if ((!verbose) && (cancelledError > epsilon)) {
969 bool hasRoot =
false;
970 beginId = starPtr->beginId;
971 endId = starPtr->endId;
972 for (iw = beginId; iw < endId; ++iw) {
973 wavePtr = &basis.wave(iw);
974 if (!(wavePtr->implicit)) {
975 rank = dftMesh.rank(wavePtr->indicesDft);
979 waveCoeff /= wavePtr->coeff;
981 diff = waveCoeff - rootCoeff;
982 if (std::abs(diff) > uncancelledError) {
983 uncancelledError = std::abs(diff);
984 if ((!verbose) && (uncancelledError > epsilon)) {
989 rootCoeff = waveCoeff;
999 if ((cancelledError < epsilon) && (uncancelledError < epsilon)) {
1001 }
else if (verbose) {
1004 <<
"Maximum coefficient of a cancelled star: "
1005 << cancelledError << std::endl
1006 <<
"Maximum error of coefficient for an uncancelled star: "
1007 << uncancelledError << std::endl;
1014 template <
int D,
class ART>
1017 IntVec<D>
const & meshDimensions,
1018 UnitCell<D>
const & unitCell,
1019 IntVec<D>
const & replicas)
1026 IntVec<D> repDimensions;
1028 for (
int i = 0; i < D; ++i) {
1031 repDimensions[i] = replicas[i] * meshDimensions[i];
1032 repSize *= repDimensions[i];
1038 for (
int i = 0; i < nMonomer; ++i) {
1043 Mesh<D> mesh(meshDimensions);
1044 Mesh<D> cellMesh(replicas);
1045 Mesh<D> repMesh(repDimensions);
1046 MeshIterator<D> iter(meshDimensions);
1047 MeshIterator<D> cellIter(replicas);
1049 IntVec<D> cellPosition;
1050 IntVec<D> repPosition;
1053 for (
int i = 0; i < nMonomer; ++i) {
1054 for (iter.begin(); !iter.atEnd(); ++iter) {
1055 position = iter.position();
1056 value = fields[i][iter.rank()];
1057 for (cellIter.begin(); !cellIter.atEnd(); ++cellIter) {
1058 cellPosition = cellIter.position();
1059 for (
int j=0; j < D; ++j) {
1060 repPosition[j] = position[j]
1061 + meshDimensions[j]*cellPosition[j];
1063 repRank = repMesh.rank(repPosition);
1064 repFields[i][repRank] = value;
1072 int nParameter = unitCell.nParameter();
1073 for (
int i = 0; i < nParameter; i++) {
1074 parameters.
append(replicas[i] * unitCell.parameter(i));
1076 cell.set(unitCell.lattice(), parameters);
1081 std::string gName =
"";
1089 template <
int D,
class ART>
1093 IntVec<D>
const & meshDimensions,
1094 UnitCell<D>
const & unitCell,
1101 template <
class ART>
1105 IntVec<1>
const & meshDimensions,
1106 UnitCell<1>
const & unitCell,
1123 cellParameters.
append(unitCell.parameter(0));
1126 std::string gName =
"";
1133 IntVec<2> dimensions;
1134 dimensions[0] = meshDimensions[0];
1135 dimensions[1] = newGridDimensions[0];
1139 if (dimensions[0] == dimensions[1]) {
1140 cell.set(UnitCell<2>::Square, cellParameters);
1142 cellParameters.
append((
double)dimensions[1]/dimensions[0]
1143 * cellParameters[0]);
1144 cell.set(UnitCell<2>::Rectangular, cellParameters);
1149 out <<
"mesh " << std::endl
1150 <<
" " << dimensions << std::endl;
1152 }
else if (d == 3) {
1157 IntVec<3> dimensions;
1158 dimensions[0] = meshDimensions[0];
1159 dimensions[1] = newGridDimensions[0];
1160 dimensions[2] = newGridDimensions[1];
1164 if (dimensions[2] == dimensions[1]) {
1165 if (dimensions[1] == dimensions[0]) {
1166 cell.set(UnitCell<3>::Cubic, cellParameters);
1168 cellParameters.
append((
double)dimensions[1]/dimensions[0]
1169 * cellParameters[0]);
1170 cellParameters.append((
double)dimensions[2]/dimensions[0]
1171 * cellParameters[0]);
1172 cell.set(UnitCell<3>::Orthorhombic, cellParameters);
1175 if (dimensions[1] == dimensions[0]) {
1176 cellParameters.append((
double)dimensions[2]/dimensions[0]
1177 * cellParameters[0]);
1178 cell.set(UnitCell<3>::Tetragonal, cellParameters);
1180 cellParameters.append((
double)dimensions[1]/dimensions[0]
1181 * cellParameters[0]);
1182 cellParameters.append((
double)dimensions[2]/dimensions[0]
1183 * cellParameters[0]);
1184 cell.set(UnitCell<3>::Orthorhombic, cellParameters);
1190 out <<
"mesh " << std::endl
1191 <<
" " << dimensions << std::endl;
1200 int nReplica = newGridDimensions[0];
1202 nReplica *= newGridDimensions[1];
1204 MeshIteratorFortran<1> iter(meshDimensions);
1206 for (
int i = 0; i < nReplica; ++i) {
1207 for (iter.begin(); !iter.atEnd(); ++iter) {
1209 for (
int j = 0; j < nMonomer; ++j) {
1210 out <<
" " <<
Dbl(fields[j][rank], 21, 13);
1218 template <
class ART>
1221 IntVec<2>
const & meshDimensions,
1222 UnitCell<2>
const & unitCell,
1237 IntVec<3> dimensions;
1238 dimensions[0] = meshDimensions[0];
1239 dimensions[1] = meshDimensions[1];
1240 dimensions[2] = newGridDimensions[0];
1245 cellParameters.
append(unitCell.parameter(0));
1246 if (unitCell.lattice() == UnitCell<2>::Square) {
1247 if (newGridDimensions[0] == meshDimensions[0]){
1248 cell.set(UnitCell<3>::Cubic, cellParameters);
1250 cellParameters.
append((
double)dimensions[2]/dimensions[0]
1251 * cellParameters[0]);
1252 cell.set(UnitCell<3>::Tetragonal, cellParameters);
1254 }
else if (unitCell.lattice() == UnitCell<2>::Rectangular) {
1255 cellParameters.
append(unitCell.parameter(1));
1256 cellParameters.
append((
double)dimensions[2]/dimensions[0]
1257 * cellParameters[0]);
1258 cell.set(UnitCell<3>::Orthorhombic, cellParameters);
1259 }
else if (unitCell.lattice() == UnitCell<2>::Hexagonal){
1260 cellParameters.append((
double)dimensions[2]/dimensions[0]
1261 * cellParameters[0]);
1262 cell.set(UnitCell<3>::Hexagonal, cellParameters);
1263 }
else if (unitCell.lattice() == UnitCell<2>::Rhombic) {
1264 cellParameters.append(unitCell.parameter(0));
1265 cellParameters.append((
double)dimensions[2]/dimensions[0]
1266 * cellParameters[0]);
1268 cellParameters.append(0.0);
1269 cellParameters.append(unitCell.parameter(1));
1270 cell.set(UnitCell<3>::Triclinic, cellParameters);
1271 }
else if (unitCell.lattice() == UnitCell<2>::Oblique) {
1272 cellParameters.append(unitCell.parameter(1));
1273 cellParameters.append((
double)dimensions[2]/dimensions[0]
1274 * cellParameters[0]);
1276 cellParameters.append(0.0);
1277 cellParameters.append(unitCell.parameter(2));
1278 cell.set(UnitCell<3>::Triclinic, cellParameters);
1280 UTIL_THROW(
"Unrecognized 2D lattice system.");
1286 std::string gName =
"";
1288 out <<
"mesh " << std::endl
1289 <<
" " << dimensions << std::endl;
1292 int nReplica = newGridDimensions[0];
1293 MeshIteratorFortran<2> iter(meshDimensions);
1295 for (
int i = 0; i < nReplica; ++i) {
1296 for (iter.begin(); !iter.atEnd(); ++iter) {
1298 for (
int j = 0; j < nMonomer; ++j) {
1299 out <<
" " <<
Dbl(fields[j][rank], 21, 13);
1307 template <
class ART>
1310 IntVec<3>
const & meshDimensions,
1311 UnitCell<3>
const & unitCell,
1315 {
UTIL_THROW(
"expandRGridDimension is invalid when D = 3."); }
An IntVec<D, T> is a D-component vector of elements of integer type T.
Iterator over points in a mesh in "Fortran" order.
Description of a regular grid of points in a periodic domain.
IntVec< D > dimensions() const
Get an IntVec<D> of the grid dimensions.
int shift(int &coordinate, int i) const
Shift a periodic coordinate into range.
A list of wavevectors that are related by space-group symmetries.
int basisId
Index of basis function associated with this star.
IntVec< D > waveBz
Integer indices indicesBz of a characteristic wave of this star.
int invertFlag
Index for inversion symmetry of star.
int size
Number of wavevectors in this star.
int beginId
Wave index of first wavevector in star.
bool cancel
Is this star cancelled, i.e., associated with a zero function?
std::complex< double > coeff
Coefficient of wave within the associated star basis function.
int starId
Index of the star that contains this wavevector.
Symmetry-adapted Fourier basis for pseudo-spectral scft.
Wave const & wave(int id) const
Get a specific Wave, access by integer index.
int waveId(IntVec< D > vector) const
Get the integer index of a wave, as required by wave(int id).
int nStar() const
Total number of stars.
bool isInitialized() const
Returns true iff this basis is fully initialized.
Star const & star(int id) const
Get a Star, accessed by integer star index.
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.
int capacity() const
Return allocated size.
static const double Pi
Trigonometric constant Pi.
Dynamically allocatable contiguous array template.
void allocate(int capacity)
Allocate the underlying C array.
bool isAllocated() const
Return true if this DArray has been allocated, false otherwise.
Wrapper for a double precision number, for formatted ostream output.
A fixed capacity (static) contiguous array with a variable logical size.
void append(Data const &data)
Append data to the end of the array.
Wrapper for an int, for formatted ostream output.
static std::ostream & file()
Get log ostream by reference.
#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.
void writeFieldHeader(std::ostream &out, int ver1, int ver2, UnitCell< D > const &cell, std::string const &groupName, int nMonomer)
Write common part of field header (fortran PSCF format).
void replicateUnitCell(std::ostream &out, DArray< AT > const &fields, IntVec< D > const &meshDimensions, UnitCell< D > const &unitCell, IntVec< D > const &replicas)
Write r-grid fields in a replicated unit cell to std::ostream.
void writeKGridData(std::ostream &out, DArray< AT > const &fields, int nMonomer, IntVec< D > const &dftDimensions)
Write data for array of k-grid fields, with no header section.
void writeNBasis(std::ostream &out, int nBasis)
Write the number of basis functions to a basis field file header.
bool hasSymmetry(AT const &in, Basis< D > const &basis, IntVec< D > const &dftDimensions, double epsilon=1.0e-8, bool verbose=true)
Check if a k-grid field has the declared space group symmetry.
void inspectArrays(DArray< AT > const &arrays, int &nMonomer, int &capacity)
Inspect dimensions of a DArray of 1D arrays, each of type AT.
void inspectFields(DArray< FT > const &fields, int &nMonomer, IntVec< D > &dimensions)
Inspect dimensions of a DArray of fields, each of type FT.
void checkAllocateArrays(DArray< AT > &arrays, int nMonomer, int capacity)
Check allocation of a DArray of 1D arrays, allocate if necessary.
void readKGridData(std::istream &in, DArray< AT > &fields, int nMonomer, IntVec< D > const &dftDimensions)
Read data for array of k-grid fields, with no header section.
void convertBasisToKGrid(DArray< double > const &components, AT &dft, Basis< D > const &basis, IntVec< D > const &dftDimensions)
Convert a real field from symmetrized basis to Fourier grid.
void writeMeshDimensions(std::ostream &out, IntVec< D > const &meshDimensions)
Write mesh dimensions to a field file header.
void expandRGridDimension(std::ostream &out, DArray< AT > const &fields, IntVec< D > const &meshDimensions, UnitCell< D > const &unitCell, int d, DArray< int > newGridDimensions)
Expand the dimensionality of space from D to d.
void writeRGridData(std::ostream &out, DArray< AT > const &fields, int nMonomer, IntVec< D > const &dimensions)
Write data for array of r-grid fields, with no header section.
void convertKGridToBasis(AT const &in, DArray< double > &out, Basis< D > const &basis, IntVec< D > const &dftDimensions, bool checkSymmetry=true, double epsilon=1.0e-8)
Convert a real field from Fourier grid to symmetrized basis.
int readNBasis(std::istream &in)
Read the number of basis functions from a basis field file header.
void readRGridData(std::istream &in, DArray< AT > &fields, int nMonomer, IntVec< D > const &dimensions)
Read data for array of r-grid fields, with no header section.
void readMeshDimensions(std::istream &in, IntVec< D > const &meshDimensions)
Read mesh dimensions from a field file header.
void writeBasisData(std::ostream &out, DArray< DArray< double > > const &fields, Basis< D > const &basis)
Write array of fields in basis format, without a header.
void readBasisData(std::istream &in, DArray< DArray< double > > &fields, UnitCell< D > const &unitCell, Mesh< D > const &mesh, Basis< D > const &basis, int nStarIn)
Read a set of fields in basis format.
void checkAllocateFields(DArray< FT > &fields, int nMonomer, IntVec< D > const &dimensions)
Check allocation of an array of fields, allocate if necessary.
void checkAllocateField(FT &field, IntVec< D > const &dimensions)
Check allocation of a single field, allocate if necessary.
RT imag(CT const &a)
Return the imaginary part of a complex number.
void assign(CT &z, RT const &a, RT const &b)
Create a complex number from real and imaginary parts, z = a + ib.
RT real(CT const &a)
Return the real part of a complex number.
PSCF package top-level namespace.