1#ifndef RPC_FIELD_IO_REAL_TPP
2#define RPC_FIELD_IO_REAL_TPP
11#include "FieldIoReal.h"
13#include <prdc/field/fieldIoUtil.h>
14#include <prdc/crystal/shiftToMinimum.h>
15#include <prdc/crystal/UnitCell.h>
16#include <prdc/crystal/Basis.h>
17#include <prdc/crystal/SpaceGroup.h>
18#include <prdc/crystal/UnitCell.h>
19#include <prdc/crystal/BFieldComparison.h>
21#include <pscf/mesh/Mesh.h>
22#include <pscf/mesh/MeshIterator.h>
23#include <pscf/mesh/MeshIteratorFortran.h>
24#include <pscf/math/IntVec.h>
26#include <util/containers/DArray.h>
27#include <util/containers/DMatrix.h>
28#include <util/misc/FileMaster.h>
29#include <util/misc/Log.h>
30#include <util/format/Str.h>
31#include <util/format/Int.h>
32#include <util/format/Dbl.h>
45 template <
int D,
class RFT,
class KFT,
class FFT>
49 hasGroupPtr_(nullptr),
50 groupNamePtr_(nullptr),
53 fileMasterPtr_(nullptr),
55 isAllocatedBasis_(false),
56 isAllocatedRGrid_(false),
57 isAllocatedKGrid_(false)
63 template <
int D,
class RFT,
class KFT,
class FFT>
72 template <
int D,
class RFT,
class KFT,
class FFT>
95 template <
int D,
class RFT,
class KFT,
class FFT>
103 template <
int D,
class RFT,
class KFT,
class FFT>
110 nMonomer_ = nMonomer;
118 template <
int D,
class RFT,
class KFT,
class FFT>
137 if (fields.isAllocated()) {
138 int nMonomerFields, fieldCapacity;
143 fields.allocate(nMonomer);
144 for (
int i = 0; i < nMonomer; ++i) {
145 fields[i].allocate(
basis().nBasis());
156 template <
int D,
class RFT,
class KFT,
class FFT>
189 template <
int D,
class RFT,
class KFT,
class FFT>
203 int nBasis = fieldCapacity;
206 bool isSymmetric =
true;
217 template <
int D,
class RFT,
class KFT,
class FFT>
246 template <
int D,
class RFT,
class KFT,
class FFT>
248 std::string filename,
262 template <
int D,
class RFT,
class KFT,
class FFT>
264 std::string filename,
277 template <
int D,
class RFT,
class KFT,
class FFT>
279 std::string filename,
292 template <
int D,
class RFT,
class KFT,
class FFT>
294 std::string filename,
307 template <
int D,
class RFT,
class KFT,
class FFT>
309 std::string filename,
324 template <
int D,
class RFT,
class KFT,
class FFT>
326 std::string filename,
341 template <
int D,
class RFT,
class KFT,
class FFT>
343 std::string filename,
346 bool isSymmetric)
const
350 bool writeHeader =
true;
351 bool writeMeshSize =
true;
353 writeHeader, isSymmetric, writeMeshSize);
360 template <
int D,
class RFT,
class KFT,
class FFT>
362 std::string filename,
365 bool isSymmetric)
const
376 template <
int D,
class RFT,
class KFT,
class FFT>
378 std::string filename,
391 template <
int D,
class RFT,
class KFT,
class FFT>
393 std::string filename,
396 bool isSymmetric)
const
409 template <
int D,
class RFT,
class KFT,
class FFT>
411 DArray< DArray <double> >
const & in,
415 int nMonomer, nMonomerOut, capacity;
424 for (
int i = 0; i < nMonomer; ++i) {
432 template <
int D,
class RFT,
class KFT,
class FFT>
435 DArray< DArray <double> > & out,
437 double epsilon)
const
440 int nMonomer, nMonomerOut, capacity;
447 for (
int i = 0; i < nMonomer; ++i) {
455 template <
int D,
class RFT,
class KFT,
class FFT>
457 std::string
const & inFileName,
458 std::string
const & outFileName)
const
471 template <
int D,
class RFT,
class KFT,
class FFT>
473 std::string
const & inFileName,
474 std::string
const & outFileName)
const
489 template <
int D,
class RFT,
class KFT,
class FFT>
498 fft().inverseTransformUnsafe(workDft_, out);
504 template <
int D,
class RFT,
class KFT,
class FFT>
506 DArray< DArray <double> >
const & in,
514 int n = in.capacity();
515 for (
int i = 0; i < n; ++i) {
517 fft().inverseTransformUnsafe(workDft_, out[i]);
524 template <
int D,
class RFT,
class KFT,
class FFT>
529 double epsilon)
const
532 fft().forwardTransform(in, workDft_);
539 template <
int D,
class RFT,
class KFT,
class FFT>
542 DArray< DArray <double> > & out,
544 double epsilon)
const
547 int nMonomer, nMonomerOut, capacity;
556 for (
int i = 0; i < nMonomer; ++i) {
557 fft().forwardTransform(in[i], workDft_);
565 template <
int D,
class RFT,
class KFT,
class FFT>
567 std::string
const & inFileName,
568 std::string
const & outFileName)
const
582 template <
int D,
class RFT,
class KFT,
class FFT>
584 std::string
const & inFileName,
585 std::string
const & outFileName)
const
600 template <
int D,
class RFT,
class KFT,
class FFT>
607 for (
int i = 0; i < n; ++i) {
608 fft().inverseTransformSafe(in[i], out[i]);
615 template <
int D,
class RFT,
class KFT,
class FFT>
617 KFT
const & in, RFT& out)
const
619 fft().inverseTransformSafe(in, out);
625 template <
int D,
class RFT,
class KFT,
class FFT>
632 for (
int i = 0; i < n; ++i) {
633 fft().forwardTransform(in[i], out[i]);
640 template <
int D,
class RFT,
class KFT,
class FFT>
644 {
fft().forwardTransform(in, out); }
649 template <
int D,
class RFT,
class KFT,
class FFT>
651 std::string
const & inFileName,
652 std::string
const & outFileName)
const
658 for (
int i = 0; i < nMonomer_; ++i) {
659 fft().inverseTransformUnsafe(tmpFieldsKGrid_[i],
668 template <
int D,
class RFT,
class KFT,
class FFT>
670 std::string
const & inFileName,
671 std::string
const & outFileName)
const
688 template <
int D,
class RFT,
class KFT,
class FFT>
695 fft().forwardTransform(in, workDft_);
702 template <
int D,
class RFT,
class KFT,
class FFT>
704 std::string
const & inFileName,
705 double epsilon)
const
714 for (
int i = 0; i < nMonomer_; ++i) {
716 symmetric =
hasSymmetry(tmpFieldsRGrid_[i], epsilon);
727 template <
int D,
class RFT,
class KFT,
class FFT>
733 comparison.
compare(field1, field2);
735 Log::file() <<
"\n Basis expansion field comparison results"
737 Log::file() <<
" Maximum Absolute Difference: "
738 << comparison.
maxDiff() << std::endl;
739 Log::file() <<
" Root-Mean-Square Difference: "
740 << comparison.
rmsDiff() <<
"\n" << std::endl;
743 template <
int D,
class RFT,
class KFT,
class FFT>
745 std::string
const & filename1,
746 std::string
const & filename2)
const
756 template <
int D,
class RFT,
class KFT,
class FFT>
758 std::string
const & filename1,
759 std::string
const & filename2)
const
774 template <
int D,
class RFT,
class KFT,
class FFT>
781 for (
int i = 0; i < n; ++i) {
789 template <
int D,
class RFT,
class KFT,
class FFT>
794 int n = fields.capacity();
796 int m = fields[0].capacity();
797 for (
int i = 0; i < n; ++i) {
806 template <
int D,
class RFT,
class KFT,
class FFT>
808 std::string
const & inFileName,
809 std::string
const & outFileName,
822 template <
int D,
class RFT,
class KFT,
class FFT>
828 for (
int i = 0; i < n; ++i) {
836 template <
int D,
class RFT,
class KFT,
class FFT>
838 std::string
const & inFileName,
839 std::string
const & outFileName,
857 template <
int D,
class RFT,
class KFT,
class FFT>
866 const int nb = fields[0].capacity();
867 for (
int i = 0; i < nMonomer_; ++i) {
877 for (i = 0; i < nb; ++i) {
878 for (j = 0; j < nMonomer_; ++j) {
880 for (k = 0; k < nMonomer_; ++k) {
881 wtmp[j] += chi(j,k)*fields[k][i];
884 for (j = 0; j < nMonomer_; ++j) {
885 fields[j][i] = wtmp[j];
894 template <
int D,
class RFT,
class KFT,
class FFT>
896 std::string
const & inFileName,
897 std::string
const & outFileName,
907 const int nb =
basis().nBasis();
933 fileMaster().openOutputFile(outFileName, out);
943 template <
int D,
class RFT,
class KFT,
class FFT>
945 std::string filename,
959 template <
int D,
class RFT,
class KFT,
class FFT>
961 std::string
const & inFileName,
962 std::string
const & outFileName,
975 template <
int D,
class RFT,
class KFT,
class FFT>
977 std::string filename,
991 template <
int D,
class RFT,
class KFT,
class FFT>
993 std::string
const & inFileName,
994 std::string
const & outFileName,
1002 d, newGridDimensions);
1017 template <
int D,
class RFT,
class KFT,
class FFT>
1022 bool & isSymmetric)
const
1035 std::string groupNameIn;
1038 groupNameIn, nMonomer);
1049 if (
lattice() != unitCell.lattice()) {
1051 <<
"Error - Mismatched lattice types "
1052 <<
"in FieldIo function readFieldHeader:\n"
1053 <<
" FieldIo::lattice :" <<
lattice() <<
"\n"
1054 <<
" Unit cell lattice :" << unitCell.lattice()
1060 isSymmetric =
false;
1061 if (groupNameIn !=
"") {
1073 <<
"Error - Mismatched group names in "
1074 <<
"FieldIo member function readFieldHeader:\n"
1075 <<
" FieldIo::groupName :" <<
groupName() <<
"\n"
1076 <<
" Field file header :" << groupNameIn <<
"\n";
1083 if (!
basis().isInitialized()) {
1084 basisPtr_->makeBasis(
mesh(), unitCell,
group());
1092 <<
"Warning: Group name found in a field file header"
1093 <<
"but no group declared in the parameter file.\n";
1103 template <
int D,
class RFT,
class KFT,
class FFT>
1108 bool isSymmetric)
const
1112 std::string gName =
"";
1127 template <
int D,
class RFT,
class KFT,
class FFT>
1130 if (isAllocatedRGrid_)
return;
1135 tmpFieldsRGrid_.allocate(nMonomer_);
1136 for (
int i = 0; i < nMonomer_; ++i) {
1137 tmpFieldsRGrid_[i].allocate(meshDimensions);
1139 isAllocatedRGrid_ =
true;
1145 template <
int D,
class RFT,
class KFT,
class FFT>
1148 if (isAllocatedKGrid_)
return;
1153 tmpFieldsKGrid_.allocate(nMonomer_);
1154 for (
int i = 0; i < nMonomer_; ++i) {
1155 tmpFieldsKGrid_[i].allocate(meshDimensions);
1157 isAllocatedKGrid_ =
true;
1163 template <
int D,
class RFT,
class KFT,
class FFT>
1166 std::string
const & inFileName)
const
1168 if (isAllocatedBasis_)
return;
1172 if (!
basis().isInitialized()) {
1175 fileMaster().openInputFile(inFileName, file);
1184 int nBasis =
basis().nBasis();
1186 tmpFieldsBasis_.allocate(nMonomer_);
1187 for (
int i = 0; i < nMonomer_; ++i) {
1190 for (
int j = 0; j < nBasis; ++j) {
1194 isAllocatedBasis_ =
true;
double rmsDiff() const
Return the precomputed root-mean-squared difference.
double compare(FT const &a, FT const &b)
Compare individual fields.
double maxDiff() const
Return the precomputed maximum element-by-element difference.
An IntVec<D, T> is a D-component vector of elements of integer type T.
Description of a regular grid of points in a periodic domain.
Comparator for fields in symmetry-adapted basis format.
Symmetry-adapted Fourier basis for pseudo-spectral SCFT.
Fourier transform wrapper.
void checkAllocateBasis(std::string const &inFileName) const
Check if basis workspace is allocated, allocate if necessary.
virtual void replicateUnitCell(std::ostream &out, DArray< RFT > const &fields, UnitCell< D > const &unitCell, IntVec< D > const &replicas) const =0
Write r-grid fields in a replicated unit cell to std::ostream.
UnitCell< D >::LatticeSystem const & lattice() const
Get the lattice type enum value by const reference.
virtual void writeFieldRGrid(std::ostream &out, RFT const &field, UnitCell< D > const &unitCell, bool writeHeader=true, bool isSymmetric=true) const =0
Write a single r-grid field to an an output stream.
void convertRGridToBasis(RFT const &in, DArray< double > &out, bool checkSymmetry=true, double epsilon=1.0e-8) const
Convert a single field from r-grid to basis form.
virtual void writeFieldsKGrid(std::ostream &out, DArray< KFT > const &fields, UnitCell< D > const &unitCell, bool isSymmetric=true) const =0
Write an array of k-grid fields to a output stream.
virtual ~FieldIoReal()
Destructor.
void scaleFieldsBasis(DArray< DArray< double > > &fields, double factor) const
Multiply an array of fields in basis format by a real scalar.
Mesh< D > const & mesh() const
Get spatial discretization mesh by const reference.
FieldIoReal()
Constructor.
virtual bool hasSymmetry(KFT const &in, double epsilon=1.0e-8, bool verbose=true) const =0
Check if a k-grid field has the declared space group symmetry.
virtual void scaleFieldBasis(DArray< double > &field, double factor) const
Multiply a single field in basis form by a real scalar.
FileMaster const & fileMaster() const
Get associated FileMaster by const reference.
void estimateWBasis(DMatrix< double > const &chi, DArray< DArray< double > > &fields) const
Convert c fields to estimated w fields, in basis format.
virtual void convertBasisToKGrid(DArray< double > const &components, KFT &dft) const =0
Convert a single field from basis to Fourier (k-grid) form.
void scaleFieldsRGrid(DArray< RFT > &fields, double factor) const
Scale an array of r-grid fields by a scalar.
void writeFieldHeader(std::ostream &out, int nMonomer, UnitCell< D > const &unitCell, bool isSymmetric=true) const
Write header for field file (fortran pscf format).
void readFieldHeader(std::istream &in, int &nMonomer, UnitCell< D > &unitCell, bool &isSymmetric) const
Reader header of field file (fortran PSCF format)
virtual void scaleFieldRGrid(RFT &field, double factor) const =0
Multiply a single field in r-grid format by a scalar.
bool hasGroup() const
Has a space group been declared externally ?
SpaceGroup< D > const & group() const
Get an associated SpaceGroup<D> by const reference.
void checkAllocateKGrid() const
Check if k-grid workspace is allocated, allocate if necessary.
void convertKGridToRGrid(DArray< KFT > const &in, DArray< RFT > &out) const
Convert an array of field from k-grid to r-grid format.
Basis< D > const & basis() const
Get the associated Basis by const reference.
void associate(Mesh< D > const &mesh, FFT const &fft, typename UnitCell< D >::LatticeSystem const &lattice, bool const &hasGroup, std::string const &groupName, SpaceGroup< D > const &group, Basis< D > &basis)
Create associations with other members of the parent Domain.
FFT const & fft() const
Get FFT object by const reference.
std::string const & groupName() const
Get an associated group name string by const reference.
virtual void convertKGridToBasis(KFT const &in, DArray< double > &out, bool checkSymmetry=true, double epsilon=1.0e-8) const =0
Convert a single field from Fourier (k-grid) to basis form.
void writeFieldsBasis(std::ostream &out, DArray< DArray< double > > const &fields, UnitCell< D > const &unitCell) const
Write an array of fields in basis format to an output stream.
void convertBasisToRGrid(DArray< double > const &in, RFT &out) const
Convert a single field from basis to r-grid format.
void setFileMaster(FileMaster const &fileMaster)
Create an association with a FileMaster.
virtual void readFieldsKGrid(std::istream &in, DArray< KFT > &fields, UnitCell< D > &unitCell) const =0
Read an array of k-grid fields from an input stream.
void setNMonomer(int nMonomer)
Set the number of monomer types.
virtual bool readFieldsRGrid(std::istream &in, DArray< RFT > &fields, UnitCell< D > &unitCell) const =0
Read array of r-grid fields from an input stream.
virtual void compareFieldsRGrid(DArray< RFT > const &field1, DArray< RFT > const &field2) const =0
Compare two fields in r-grid form, write a report to Log file.
void compareFieldsBasis(DArray< DArray< double > > const &field1, DArray< DArray< double > > const &field2) const
Compare array of fields in basis form, write a report to Log file.
void checkAllocateRGrid() const
Check if r-grid workspace is allocated, allocate if necessary.
void readFieldsBasis(std::istream &in, DArray< DArray< double > > &fields, UnitCell< D > &unitCell) const
Read an array of fields in basis format from an input stream.
void convertRGridToKGrid(DArray< RFT > const &in, DArray< KFT > &out) const
Convert an array of fields from r-grid to k-grid (Fourier) format.
void writeFieldBasis(std::ostream &out, DArray< double > const &field, UnitCell< D > const &unitCell) const
Write a single field in basis format to an output stream.
void readFieldBasis(std::istream &in, DArray< double > &field, UnitCell< D > &unitCell) const
Read a single field in basis format from an input stream.
virtual void writeFieldsRGrid(std::ostream &out, DArray< RFT > const &fields, UnitCell< D > const &unitCell, bool writeHeader=true, bool isSymmetric=true, bool writeMeshSize=true) const =0
Write an array of r-grid fields to an output stream.
virtual void expandRGridDimension(std::ostream &out, DArray< RFT > const &fields, UnitCell< D > const &unitCell, int d, DArray< int > const &newGridDimensions) const =0
Increase D for an array of r-grid fields, write to ostream.
virtual bool readFieldRGrid(std::istream &in, RFT &field, UnitCell< D > &unitCell) const =0
Read a single r-grid field from an input stream.
Crystallographic space group.
bool isInitialized() const
Has this unit cell been initialized?
int nParameter() const
Get the number of parameters in the unit cell.
Base template for UnitCell<D> classes, D=1, 2 or 3.
int capacity() const
Return allocated size.
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.
Dynamically allocated Matrix.
A FileMaster manages input and output files for a simulation.
static std::ostream & file()
Get log ostream by reference.
int capacity2() const
Get number of columns (range of the second array index).
int capacity1() const
Get number of rows (range of the first array index).
#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 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 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 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 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 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 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 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 checkAllocateField(FT &field, IntVec< D > const &dimensions)
Check allocation of a single field, allocate if necessary.
Periodic fields and crystallography.
PSCF package top-level namespace.