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/unitCellHeader.h>
14#include <prdc/crystal/replicateUnitCell.h>
15#include <prdc/crystal/shiftToMinimum.h>
17#include <pscf/mesh/Mesh.h>
18#include <pscf/mesh/MeshIterator.h>
19#include <pscf/mesh/MeshIteratorFortran.h>
20#include <pscf/math/IntVec.h>
21#include <pscf/math/complex.h>
23#include <util/containers/DArray.h>
24#include <util/misc/FileMaster.h>
25#include <util/misc/Log.h>
26#include <util/format/Int.h>
27#include <util/format/Dbl.h>
28#include <util/math/Constants.h>
40 template <
int D,
class FT>
44 if (field.isAllocated()) {
45 UTIL_CHECK(field.meshDimensions() == dimensions);
47 field.allocate(dimensions);
54 template <
int D,
class FT>
60 int nMonomerFields = fields.
capacity();
63 for (
int i = 0; i < nMonomer; ++i) {
65 UTIL_CHECK(fields[i].meshDimensions() == dimensions);
69 for (
int i = 0; i < nMonomer; ++i) {
78 template <
int D,
class FT>
87 dimensions = fields[0].meshDimensions();
88 for (
int i = 0; i < nMonomer; ++i) {
90 UTIL_CHECK(fields[i].meshDimensions() == dimensions);
103 int nMonomerArrays = arrays.
capacity();
106 for (
int i = 0; i < nMonomer; ++i) {
112 for (
int i = 0; i < nMonomer; ++i) {
132 for (
int i = 0; i < nMonomer; ++i) {
145 int& ver1,
int& ver2,
147 std::string& groupName,
172 if (label ==
"group_name") {
190 std::string
const & groupName,
194 out <<
"format " <<
Int(ver1,3) <<
" " <<
Int(ver2,3) << std::endl;
195 out <<
"dim" << std::endl
196 <<
" " << D << std::endl;
203 if (groupName !=
"") {
204 out <<
"group_name" << std::endl
205 <<
" " << groupName << std::endl;
211 out <<
"N_monomer" << std::endl
212 <<
" " << nMonomer << std::endl;
234 if (label !=
"mesh" && label !=
"ngrid") {
235 std::string msg =
"\n";
236 msg +=
"Error reading field file:\n";
237 msg +=
"Expected mesh or ngrid, but found [";
243 in >> meshDimensionsIn;
245 if (meshDimensionsIn != meshDimensions) {
247 <<
"Inconsistent mesh dimensions:\n"
248 <<
"meshDimensions (expected) = " << meshDimensions <<
"\n"
249 <<
"meshDimensions (from file) = " << meshDimensionsIn <<
"\n";
250 UTIL_THROW(
"Unexpected mesh dimensions in field file header");
258 out <<
"mesh " << std::endl
259 <<
" " << meshDimensions << std::endl;
262 template <
int D,
class ART>
273 for (iter.begin(); !iter.atEnd(); ++iter) {
275 for (
int k = 0; k < nMonomer; ++k) {
276 in >> fields[k][rank];
283 template <
int D,
class ART>
290 for (iter.begin(); !iter.atEnd(); ++iter) {
298 template <
int D,
class ART>
309 for (iter.begin(); !iter.atEnd(); ++iter) {
311 for (j = 0; j < nMonomer; ++j) {
312 out <<
" " <<
Dbl(fields[j][rank], 21, 13);
319 template <
int D,
class ART>
326 for (iter.begin(); !iter.atEnd(); ++iter) {
328 out <<
" " <<
Dbl(field[rank], 21, 13);
335 template <
int D,
class ACT>
341 typedef typename ACT::Complex CT;
342 typedef typename ACT::Real RT;
346 int rank, i, j, idum;
348 for (iter.begin(); !iter.atEnd(); ++iter) {
353 for (j = 0; j < nMonomer; ++j) {
364 template <
int D,
class ACT>
369 typedef typename ACT::Complex CT;
370 typedef typename ACT::Real RT;
376 for (iter.begin(); !iter.atEnd(); ++iter) {
390 template <
int D,
class ACT>
399 typedef typename ACT::Complex CT;
400 typedef typename ACT::Real RT;
406 for (iter.begin(); !iter.atEnd(); ++iter) {
410 for (
int j = 0; j < nMonomer; ++j) {
423 template <
int D,
class ACT>
428 typedef typename ACT::Complex CT;
429 typedef typename ACT::Real RT;
435 for (iter.begin(); !iter.atEnd(); ++iter) {
469 int fieldCapacity = fields[0].
capacity();
471 for (
int i = 0; i < nMonomer; ++i) {
472 UTIL_CHECK(fields[i].capacity() == fieldCapacity);
473 for (
int j = 0; j < fieldCapacity; ++j) {
479 int nBasis = nBasisIn;
480 if (fieldCapacity < nBasisIn) {
481 nBasis = fieldCapacity;
494 int basisId, basisId2;
497 std::complex<double> coeff, phasor;
499 int nReversedPair = 0;
500 bool waveExists, sizeMatches;
507 for (
int j = 0; j < nMonomer; ++j) {
519 waveBz = shiftToMinimum(waveIn, mesh.
dimensions(), unitCell);
520 waveExists = (waveIn == waveBz);
527 waveId = basis.
waveId(waveDft);
529 starPtr = &basis.
star(starId);
533 if (starPtr->
size == sizeIn) {
537 <<
"Warning: Inconsistent star size (line ignored)\n"
538 <<
"wave from file = " << waveIn <<
"\n"
539 <<
"size from file = " << sizeIn <<
"\n"
540 <<
"size of star = " << starPtr->
size
547 if (waveExists && sizeMatches) {
551 if (starPtr->
waveBz == waveIn) {
554 for (
int j = 0; j < nMonomer; ++j) {
555 fields[j][basisId] = temp[j];
561 <<
"Inconsistent wave of closed star on input\n"
562 <<
"wave from file = " << waveIn <<
"\n"
563 <<
"starId of wave = " << starId <<
"\n"
564 <<
"waveBz of star = " << starPtr->
waveBz
572 for (
int j = 0; j < nMonomer; ++j) {
582 shiftToMinimum(waveIn2, mesh.
dimensions(), unitCell);
588 waveId2 = basis.
waveId(waveDft);
590 starPtr2 = &basis.
star(starId2);
607 for (
int j = 0; j < nMonomer; ++j) {
608 fields[j][basisId] = temp[j];
609 fields[j][basisId2] = temp2[j];
626 shiftToMinimum(nVec, mesh.
dimensions(), unitCell);
675 phasor = phasor/std::abs(phasor);
676 for (
int j = 0; j < nMonomer; ++j) {
677 coeff = std::complex<double>(temp[j],-temp2[j]);
679 fields[j][basisId2] =
real(coeff);
680 fields[j][basisId ] =
imag(coeff);
696 if (nReversedPair > 0) {
698 Log::file() << nReversedPair <<
" reversed pairs of open stars"
699 <<
" detected in readFieldsBasis\n";
740 int fieldCapacity = fields[0].
capacity();
741 for (
int j = 0; j < nMonomer; ++j) {
742 UTIL_CHECK(fieldCapacity == fields[j].capacity());
748 int nStar = basis.
nStar();
749 for (
int i = 0; i < nStar; ++i) {
750 if (ib >= fieldCapacity)
break;
752 for (
int j = 0; j < nMonomer; ++j) {
753 out <<
Dbl(fields[j][ib], 20, 10);
756 for (
int j = 0; j < D; ++j) {
789 template <
int D,
class ACT>
792 Basis<D>
const& basis,
797 typedef typename ACT::Complex CT;
798 typedef typename ACT::Real RT;
801 Mesh<D> dftMesh(dftDimensions);
803 typename Basis<D>::Star
const* starPtr;
804 typename Basis<D>::Wave
const* wavePtr;
805 std::complex<double> component;
806 std::complex<double> coeff;
814 for (rank = 0; rank < dftMesh.size(); ++rank) {
820 while (is < basis.nStar()) {
821 starPtr = &(basis.star(is));
823 if (starPtr->cancel) {
829 ib = starPtr->basisId;
831 if (starPtr->invertFlag == 0) {
834 component = std::complex<double>(in[ib], 0.0);
837 for (iw = starPtr->beginId; iw < starPtr->endId; ++iw) {
838 wavePtr = &basis.wave(iw);
839 if (!wavePtr->implicit) {
840 coeff = component*(wavePtr->coeff);
841 indices = wavePtr->indicesDft;
842 rank = dftMesh.rank(indices);
851 if (starPtr->invertFlag == 1) {
854 component = std::complex<double>(in[ib], -in[ib+1]);
855 component /= sqrt(2.0);
856 starPtr = &(basis.star(is));
857 for (iw = starPtr->beginId; iw < starPtr->endId; ++iw) {
858 wavePtr = &basis.wave(iw);
859 if (!(wavePtr->implicit)) {
860 coeff = component*(wavePtr->coeff);
861 indices = wavePtr->indicesDft;
862 rank = dftMesh.rank(indices);
870 starPtr = &(basis.star(is+1));
872 component = std::complex<double>(in[ib], +in[ib+1]);
873 component /= sqrt(2.0);
874 for (iw = starPtr->beginId; iw < starPtr->endId; ++iw) {
875 wavePtr = &basis.wave(iw);
876 if (!(wavePtr->implicit)) {
877 coeff = component*(wavePtr->coeff);
878 indices = wavePtr->indicesDft;
879 rank = dftMesh.rank(indices);
899 template <
int D,
class ACT>
909 typedef typename ACT::Complex CT;
910 typedef typename ACT::Real RT;
915 symmetric =
hasSymmetry(in, basis, dftDimensions, epsilon,
true);
918 <<
"WARNING: Conversion of asymmetric field to"
919 <<
"symmetry-adapted basis in Prdc::convertKgridToBasis."
920 << std::endl << std::endl;
925 Mesh<D> dftMesh(dftDimensions);
929 std::complex<double> component;
936 for (is = 0; is < basis.nBasis(); ++is) {
942 while (is < basis.nStar()) {
943 starPtr = &(basis.star(is));
945 if (starPtr->cancel) {
951 ib = starPtr->basisId;
953 if (starPtr->invertFlag == 0) {
956 int beginId = starPtr->beginId;
957 int endId = starPtr->endId;
959 bool isImplicit =
true;
961 wavePtr = &basis.wave(beginId + iw);
962 if (!wavePtr->implicit) {
966 wavePtr = &basis.wave(endId - 1 - iw);
967 if (!wavePtr->implicit) {
976 rank = dftMesh.rank(wavePtr->indicesDft);
979 component /= wavePtr->coeff;
980 out[ib] = component.real();
984 if (starPtr->invertFlag == 1) {
989 wavePtr = &basis.wave(starPtr->beginId);
991 if (wavePtr->implicit) {
992 iw = wavePtr->inverseId;
993 starPtr = &(basis.star(is+1));
995 wavePtr = &basis.wave(iw);
999 rank = dftMesh.rank(wavePtr->indicesDft);
1002 UTIL_CHECK(std::abs(wavePtr->coeff) > 1.0E-8);
1003 component /= wavePtr->coeff;
1004 component *= sqrt(2.0);
1007 if (starPtr->invertFlag == 1) {
1008 out[ib] = component.real();
1009 out[ib+1] = -component.imag();
1011 out[ib] = component.real();
1012 out[ib+1] = component.imag();
1028 template <
int D,
class ACT>
1037 typedef typename ACT::Complex CT;
1038 typedef typename ACT::Real RT;
1042 std::complex<double> waveCoeff;
1043 std::complex<double> rootCoeff;
1044 std::complex<double> diff;
1050 double cancelledError(0.0);
1051 double uncancelledError(0.0);
1054 Mesh<D> dftMesh(dftDimensions);
1057 for (is = 0; is < basis.nStar(); ++is) {
1058 starPtr = &(basis.star(is));
1060 if (starPtr->cancel) {
1063 beginId = starPtr->beginId;
1064 endId = starPtr->endId;
1065 for (iw = beginId; iw < endId; ++iw) {
1066 wavePtr = &basis.wave(iw);
1067 if (!wavePtr->implicit) {
1068 rank = dftMesh.rank(wavePtr->indicesDft);
1072 if (std::abs(waveCoeff) > cancelledError) {
1073 cancelledError = std::abs(waveCoeff);
1074 if ((!verbose) && (cancelledError > epsilon)) {
1084 bool hasRoot =
false;
1085 beginId = starPtr->beginId;
1086 endId = starPtr->endId;
1087 for (iw = beginId; iw < endId; ++iw) {
1088 wavePtr = &basis.wave(iw);
1089 if (!(wavePtr->implicit)) {
1090 rank = dftMesh.rank(wavePtr->indicesDft);
1094 waveCoeff /= wavePtr->coeff;
1096 diff = waveCoeff - rootCoeff;
1097 if (std::abs(diff) > uncancelledError) {
1098 uncancelledError = std::abs(diff);
1099 if ((!verbose) && (uncancelledError > epsilon)) {
1104 rootCoeff = waveCoeff;
1114 if ((cancelledError < epsilon) && (uncancelledError < epsilon)) {
1116 }
else if (verbose) {
1119 <<
"Maximum coefficient of a cancelled star: "
1120 << cancelledError << std::endl
1121 <<
"Maximum error of coefficient for an uncancelled star: "
1122 << uncancelledError << std::endl;
1129 template <
int D,
class ART>
1143 for (
int i = 0; i < D; ++i) {
1146 repDimensions[i] = replicas[i] * meshDimensions[i];
1147 repSize *= repDimensions[i];
1153 for (
int i = 0; i < nMonomer; ++i) {
1160 Mesh<D> repMesh(repDimensions);
1168 for (
int i = 0; i < nMonomer; ++i) {
1169 for (iter.begin(); !iter.atEnd(); ++iter) {
1170 position = iter.position();
1171 value = fields[i][iter.rank()];
1172 for (cellIter.begin(); !cellIter.atEnd(); ++cellIter) {
1173 cellPosition = cellIter.position();
1174 for (
int j=0; j < D; ++j) {
1175 repPosition[j] = position[j]
1176 + meshDimensions[j]*cellPosition[j];
1178 repRank = repMesh.rank(repPosition);
1179 repFields[i][repRank] = value;
1192 std::string gName =
"";
1200 template <
int D,
class ART>
1212 template <
class ART>
1234 cellParameters.
append(unitCell.parameter(0));
1237 std::string gName =
"";
1245 dimensions[0] = meshDimensions[0];
1246 dimensions[1] = newGridDimensions[0];
1250 if (dimensions[0] == dimensions[1]) {
1253 cellParameters.
append((
double)dimensions[1]/dimensions[0]
1254 * cellParameters[0]);
1260 out <<
"mesh " << std::endl
1261 <<
" " << dimensions << std::endl;
1263 }
else if (d == 3) {
1269 dimensions[0] = meshDimensions[0];
1270 dimensions[1] = newGridDimensions[0];
1271 dimensions[2] = newGridDimensions[1];
1275 if (dimensions[2] == dimensions[1]) {
1276 if (dimensions[1] == dimensions[0]) {
1279 cellParameters.
append((
double)dimensions[1]/dimensions[0]
1280 * cellParameters[0]);
1281 cellParameters.append((
double)dimensions[2]/dimensions[0]
1282 * cellParameters[0]);
1286 if (dimensions[1] == dimensions[0]) {
1287 cellParameters.append((
double)dimensions[2]/dimensions[0]
1288 * cellParameters[0]);
1291 cellParameters.append((
double)dimensions[1]/dimensions[0]
1292 * cellParameters[0]);
1293 cellParameters.append((
double)dimensions[2]/dimensions[0]
1294 * cellParameters[0]);
1301 out <<
"mesh " << std::endl
1302 <<
" " << dimensions << std::endl;
1311 int nReplica = newGridDimensions[0];
1313 nReplica *= newGridDimensions[1];
1317 for (
int i = 0; i < nReplica; ++i) {
1318 for (iter.begin(); !iter.atEnd(); ++iter) {
1320 for (
int j = 0; j < nMonomer; ++j) {
1321 out <<
" " <<
Dbl(fields[j][rank], 21, 13);
1329 template <
class ART>
1349 dimensions[0] = meshDimensions[0];
1350 dimensions[1] = meshDimensions[1];
1351 dimensions[2] = newGridDimensions[0];
1356 cellParameters.
append(unitCell.parameter(0));
1358 if (newGridDimensions[0] == meshDimensions[0]){
1361 cellParameters.
append((
double)dimensions[2]/dimensions[0]
1362 * cellParameters[0]);
1366 cellParameters.
append(unitCell.parameter(1));
1367 cellParameters.
append((
double)dimensions[2]/dimensions[0]
1368 * cellParameters[0]);
1371 cellParameters.append((
double)dimensions[2]/dimensions[0]
1372 * cellParameters[0]);
1375 cellParameters.append(unitCell.parameter(0));
1376 cellParameters.append((
double)dimensions[2]/dimensions[0]
1377 * cellParameters[0]);
1379 cellParameters.append(0.0);
1380 cellParameters.append(unitCell.parameter(1));
1383 cellParameters.append(unitCell.parameter(1));
1384 cellParameters.append((
double)dimensions[2]/dimensions[0]
1385 * cellParameters[0]);
1387 cellParameters.append(0.0);
1388 cellParameters.append(unitCell.parameter(2));
1391 UTIL_THROW(
"Unrecognized 2D lattice system.");
1397 std::string gName =
"";
1399 out <<
"mesh " << std::endl
1400 <<
" " << dimensions << std::endl;
1403 int nReplica = newGridDimensions[0];
1406 for (
int i = 0; i < nReplica; ++i) {
1407 for (iter.begin(); !iter.atEnd(); ++iter) {
1409 for (
int j = 0; j < nMonomer; ++j) {
1410 out <<
" " <<
Dbl(fields[j][rank], 21, 13);
1418 template <
class ART>
1426 {
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.
Iterator over points in a Mesh<D>.
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?
Wavevector used to construct a basis 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.
#define UTIL_ASSERT(condition)
Assertion macro suitable for debugging serial or parallel code.
void readUnitCellHeader(std::istream &in, UnitCell< D > &cell)
Read UnitCell<D> from a field file header (fortran PSCF format).
void replicateUnitCell(IntVec< 1 > const &replicas, UnitCell< 1 > const &cellIn, UnitCell< 1 > &cellOut)
Create a replicated UnitCell<1>.
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 writeUnitCellHeader(std::ostream &out, UnitCell< D > const &cell)
Write UnitCell<D> to a field file header (fortran PSCF format).
void readFieldHeader(std::istream &in, int &ver1, int &ver2, UnitCell< D > &cell, std::string &groupName, int &nMonomer)
Read common part of field header (fortran PSCF format).
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.
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.
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 data section of an array of fields in basis format.
void readBasisData(std::istream &in, DArray< DArray< double > > &fields, UnitCell< D > const &unitCell, Mesh< D > const &mesh, Basis< D > const &basis, int nStarIn)
Read an array of fields in basis format, without a header.
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.
Periodic fields and crystallography.
PSCF package top-level namespace.