PSCF v1.3
FieldIoReal.tpp
1#ifndef RPC_FIELD_IO_REAL_TPP
2#define RPC_FIELD_IO_REAL_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 "FieldIoReal.h"
12
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>
20
21#include <pscf/mesh/Mesh.h>
22#include <pscf/mesh/MeshIterator.h>
23#include <pscf/mesh/MeshIteratorFortran.h>
24#include <pscf/math/IntVec.h>
25
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>
33
34#include <iomanip>
35#include <string>
36
37namespace Pscf {
38namespace Prdc {
39
40 using namespace Util;
41
42 /*
43 * Constructor.
44 */
45 template <int D, class RFT, class KFT, class FFT>
47 : meshPtr_(nullptr),
48 fftPtr_(nullptr),
49 hasGroupPtr_(nullptr),
50 groupNamePtr_(nullptr),
51 groupPtr_(nullptr),
52 basisPtr_(nullptr),
53 fileMasterPtr_(nullptr),
54 nMonomer_(0),
55 isAllocatedBasis_(false),
56 isAllocatedRGrid_(false),
57 isAllocatedKGrid_(false)
58 {}
59
60 /*
61 * Destructor.
62 */
63 template <int D, class RFT, class KFT, class FFT>
66
67 // Initialization functions
68
69 /*
70 * Create associations with other members of a parent Domain<D> object.
71 */
72 template <int D, class RFT, class KFT, class FFT>
73 void
75 Mesh<D> const & mesh,
76 FFT const & fft,
77 typename UnitCell<D>::LatticeSystem const & lattice,
78 bool const & hasGroup,
79 std::string const & groupName,
80 SpaceGroup<D> const & group,
82 {
83 meshPtr_ = &mesh;
84 fftPtr_ = &fft;
85 latticePtr_ = &lattice;
86 hasGroupPtr_ = &hasGroup;
87 groupNamePtr_ = &groupName;
88 groupPtr_ = &group;
89 basisPtr_ = &basis;
90 }
91
92 /*
93 * Create an association with a FileMaster.
94 */
95 template <int D, class RFT, class KFT, class FFT>
99
100 /*
101 * Set nMonomer, the number of monomer types.
102 */
103 template <int D, class RFT, class KFT, class FFT>
105 {
106 // Preconditions - require that function is only called once
107 UTIL_CHECK(nMonomer_ == 0);
108 UTIL_CHECK(nMonomer > 0);
109
110 nMonomer_ = nMonomer;
111 }
112
113 // Field File IO - Symmetry-Adapted Basis Format
114
115 /*
116 * Read an array of fields in basis format from an input stream.
117 */
118 template <int D, class RFT, class KFT, class FFT>
120 std::istream& in,
121 DArray< DArray<double> >& fields,
122 UnitCell<D>& unitCell) const
123 {
124 // Precondition
126
127 // Read header (checks compatibility with space group)
128 int nMonomer;
129 bool isSymmetric;
130 readFieldHeader(in, nMonomer, unitCell, isSymmetric);
131 UTIL_CHECK(isSymmetric);
132 UTIL_CHECK(basis().isInitialized());
133 // Note: readFieldHeader can initialize basis if not done previously
134 int nBasisIn = readNBasis(in);
135
136 // Check allocation of fields container
137 if (fields.isAllocated()) {
138 int nMonomerFields, fieldCapacity;
139 inspectArrays(fields, nMonomerFields, fieldCapacity);
140 UTIL_CHECK(nMonomerFields == nMonomer);
141 UTIL_CHECK(fieldCapacity == basis().nBasis());
142 } else {
143 fields.allocate(nMonomer);
144 for (int i = 0; i < nMonomer; ++i) {
145 fields[i].allocate(basis().nBasis());
146 }
147 }
148
149 // Read field data (components in basis)
150 Prdc::readBasisData(in, fields, unitCell, mesh(), basis(), nBasisIn);
151 }
152
153 /*
154 * Read a single field in basis format from an input stream
155 */
156 template <int D, class RFT, class KFT, class FFT>
158 std::istream& in,
159 DArray<double>& field,
160 UnitCell<D>& unitCell) const
161 {
162 // Local array of fields, of type required by readFieldsBasis
163 DArray< DArray<double> > fields;
164
165 // If single field is allocated, allocate local array fields
166 if (field.isAllocated()) {
167 int capacity = field.capacity();
168 UTIL_CHECK(basis().nBasis() == capacity);
169 fields.allocate(1);
170 fields[0].allocate(capacity);
171 }
172 // Otherwise, pass unallocated fields array to readFieldsBasis,
173 // which will allocate it as necessary.
174
175 // Read file containing a single field
176 readFieldsBasis(in, fields, unitCell);
177
178 // Check that only one field was read from file
179 UTIL_CHECK(fields.capacity() == 1);
180
181 // Copy data from local array fields to function parameter field.
182 // The assignment operator will allocate the field if necessary.
183 field = fields[0];
184 }
185
186 /*
187 * Write an array of fields in basis format to an output stream.
188 */
189 template <int D, class RFT, class KFT, class FFT>
191 std::ostream &out,
192 DArray< DArray<double> > const & fields,
193 UnitCell<D> const & unitCell) const
194 {
195 // Inspect fields to obtain nMonomer and fieldCapacity
196 int nMonomer;
197 int fieldCapacity;
198 inspectArrays(fields, nMonomer, fieldCapacity);
199
200 // Preconditions
201 UTIL_CHECK(basis().isInitialized());
202 UTIL_CHECK(fieldCapacity <= basis().nBasis());
203 int nBasis = fieldCapacity;
204
205 // Write header
206 bool isSymmetric = true;
207 writeFieldHeader(out, nMonomer, unitCell, isSymmetric);
208 writeNBasis(out, nBasis);
209
210 // Write data (field components)
211 Prdc::writeBasisData(out, fields, basis());
212 }
213
214 /*
215 * Write a single field in basis format to an output stream.
216 */
217 template <int D, class RFT, class KFT, class FFT>
219 std::ostream& out,
220 DArray<double> const & field,
221 UnitCell<D> const & unitCell) const
222 {
223 // Create local array of type required by writeFieldsBasis
224 DArray<DArray<double> > fields;
225 fields.allocate(1);
226 fields[0].allocate(field.capacity());
227
228 // Copy data from input parameter to local array
229 fields[0] = field;
230
231 writeFieldsBasis(out, fields, unitCell);
232 }
233
234 /*
235 * File IO wrapper functions:
236 *
237 * These functions take a file name as an argument, and simply wrap
238 * file open and close operations around a function of the same name
239 * that takes an io stream argument. These functions can use the same
240 * implementation in Rpc and Rpg.
241 */
242
243 /*
244 * Open-close a file and read a set of fields in basis format.
245 */
246 template <int D, class RFT, class KFT, class FFT>
248 std::string filename,
249 DArray<DArray<double> >& fields,
250 UnitCell<D>& unitCell) const
251 {
252
253 std::ifstream file;
254 fileMaster().openInputFile(filename, file);
255 readFieldsBasis(file, fields, unitCell);
256 file.close();
257 }
258
259 /*
260 * Open-close a file and read a single field in basis format.
261 */
262 template <int D, class RFT, class KFT, class FFT>
264 std::string filename,
265 DArray<double>& field,
266 UnitCell<D>& unitCell) const
267 {
268 std::ifstream file;
269 fileMaster().openInputFile(filename, file);
270 readFieldBasis(file, field, unitCell);
271 file.close();
272 }
273
274 /*
275 * Open-close a file, and write an array of fields in basis format.
276 */
277 template <int D, class RFT, class KFT, class FFT>
279 std::string filename,
280 DArray<DArray<double> > const & fields,
281 UnitCell<D> const & unitCell) const
282 {
283 std::ofstream file;
284 fileMaster().openOutputFile(filename, file);
285 writeFieldsBasis(file, fields, unitCell);
286 file.close();
287 }
288
289 /*
290 * Open-close a file, and write a single field in basis format.
291 */
292 template <int D, class RFT, class KFT, class FFT>
294 std::string filename,
295 DArray<double> const & field,
296 UnitCell<D> const & unitCell) const
297 {
298 std::ofstream file;
299 fileMaster().openOutputFile(filename, file);
300 writeFieldBasis(file, field, unitCell);
301 file.close();
302 }
303
304 /*
305 * Open-close a file, and write an array of fields in r-grid format
306 */
307 template <int D, class RFT, class KFT, class FFT>
309 std::string filename,
310 DArray< RFT >& fields,
311 UnitCell<D>& unitCell) const
312 {
313 std::ifstream file;
314 fileMaster().openInputFile(filename, file);
315 bool isSymmetric;
316 isSymmetric = readFieldsRGrid(file, fields, unitCell);
317 file.close();
318 return isSymmetric;
319 }
320
321 /*
322 * Open-close a file, and read a single field in r-grid format
323 */
324 template <int D, class RFT, class KFT, class FFT>
326 std::string filename,
327 RFT & field,
328 UnitCell<D>& unitCell) const
329 {
330 std::ifstream file;
331 fileMaster().openInputFile(filename, file);
332 bool isSymmetric;
333 isSymmetric = readFieldRGrid(file, field, unitCell);
334 file.close();
335 return isSymmetric;
336 }
337
338 /*
339 * Open-close a file, and write an array of fields in r-grid format
340 */
341 template <int D, class RFT, class KFT, class FFT>
343 std::string filename,
344 DArray< RFT > const & fields,
345 UnitCell<D> const & unitCell,
346 bool isSymmetric) const
347 {
348 std::ofstream file;
349 fileMaster().openOutputFile(filename, file);
350 bool writeHeader = true;
351 bool writeMeshSize = true;
352 writeFieldsRGrid(file, fields, unitCell,
353 writeHeader, isSymmetric, writeMeshSize);
354 file.close();
355 }
356
357 /*
358 * Open-close a file, and write a single field in r-grid format
359 */
360 template <int D, class RFT, class KFT, class FFT>
362 std::string filename,
363 RFT const & field,
364 UnitCell<D> const & unitCell,
365 bool isSymmetric) const
366 {
367 std::ofstream file;
368 fileMaster().openOutputFile(filename, file);
369 writeFieldRGrid(file, field, unitCell, isSymmetric);
370 file.close();
371 }
372
373 /*
374 * Open-close a file, and read an array of fields in k-grid format
375 */
376 template <int D, class RFT, class KFT, class FFT>
378 std::string filename,
379 DArray< KFT >& fields,
380 UnitCell<D>& unitCell) const
381 {
382 std::ifstream file;
383 fileMaster().openInputFile(filename, file);
384 readFieldsKGrid(file, fields, unitCell);
385 file.close();
386 }
387
388 /*
389 * Open-close a file, and read an array of fields in k-grid format
390 */
391 template <int D, class RFT, class KFT, class FFT>
393 std::string filename,
394 DArray< KFT > const & fields,
395 UnitCell<D> const & unitCell,
396 bool isSymmetric) const
397 {
398 std::ofstream file;
399 fileMaster().openOutputFile(filename, file);
400 writeFieldsKGrid(file, fields, unitCell, isSymmetric);
401 file.close();
402 }
403
404 // Field Format Conversion Functions - Basis <-> KGrid
405
406 /*
407 * Convert array of fields from basis to k-grid format.
408 */
409 template <int D, class RFT, class KFT, class FFT>
411 DArray< DArray <double> > const & in,
412 DArray< KFT >& out) const
413 {
414 // Inspect input and output field containers
415 int nMonomer, nMonomerOut, capacity;
416 inspectArrays(in, nMonomer, capacity);
417 IntVec<D> dimensions;
418 inspectFields(out, nMonomerOut, dimensions);
419 UTIL_CHECK(nMonomer == nMonomerOut);
420 UTIL_CHECK(capacity == basis().nBasis());
421 UTIL_CHECK(dimensions == mesh().dimensions());
422
423 // Convert fields for all monomer types
424 for (int i = 0; i < nMonomer; ++i) {
425 convertBasisToKGrid(in[i], out[i]);
426 }
427 }
428
429 /*
430 * Convert array of fields from k-grid format to basis format.
431 */
432 template <int D, class RFT, class KFT, class FFT>
434 DArray< KFT > const & in,
435 DArray< DArray <double> > & out,
436 bool checkSymmetry,
437 double epsilon) const
438 {
439 // Inspect input and output containers
440 int nMonomer, nMonomerOut, capacity;
441 IntVec<D> dimensions;
442 inspectFields(in, nMonomer, dimensions);
443 inspectArrays(out, nMonomerOut, capacity);
444 UTIL_CHECK(nMonomer == nMonomerOut);
445
446 // Convert fields for all monomer types
447 for (int i = 0; i < nMonomer; ++i) {
448 convertKGridToBasis(in[i], out[i], checkSymmetry, epsilon);
449 }
450 }
451
452 /*
453 * Convert a field file from k-grid to basis format.
454 */
455 template <int D, class RFT, class KFT, class FFT>
457 std::string const & inFileName,
458 std::string const & outFileName) const
459 {
461 checkAllocateBasis(inFileName);
462 UnitCell<D> tmpUnitCell;
463 readFieldsKGrid(inFileName, tmpFieldsKGrid_, tmpUnitCell);
464 convertKGridToBasis(tmpFieldsKGrid_, tmpFieldsBasis_);
465 writeFieldsBasis(outFileName, tmpFieldsBasis_, tmpUnitCell);
466 }
467
468 /*
469 * Convert a field file from basis to k-grid format.
470 */
471 template <int D, class RFT, class KFT, class FFT>
473 std::string const & inFileName,
474 std::string const & outFileName) const
475 {
477 checkAllocateBasis(inFileName);
478 UnitCell<D> tmpUnitCell;
479 readFieldsBasis(inFileName, tmpFieldsBasis_, tmpUnitCell);
480 convertBasisToKGrid(tmpFieldsBasis_, tmpFieldsKGrid_);
481 writeFieldsKGrid(outFileName, tmpFieldsKGrid_, tmpUnitCell);
482 }
483
484 // Field Format Conversion Functions - Basis <-> RGrid
485
486 /*
487 * Convert a single field from basis to r-grid format.
488 */
489 template <int D, class RFT, class KFT, class FFT>
491 DArray<double> const & in,
492 RFT& out) const
493 {
495 UTIL_CHECK(out.isAllocated());
496 checkAllocateField(workDft_, mesh().dimensions());
497 convertBasisToKGrid(in, workDft_);
498 fft().inverseTransformUnsafe(workDft_, out);
499 }
500
501 /*
502 * Convert an array of fields from basis to r-grid format.
503 */
504 template <int D, class RFT, class KFT, class FFT>
506 DArray< DArray <double> > const & in,
507 DArray< RFT >& out) const
508 {
509 UTIL_CHECK(in.isAllocated());
510 UTIL_CHECK(out.isAllocated());
511 UTIL_CHECK(in.capacity() == out.capacity());
512 checkAllocateField(workDft_, mesh().dimensions());
513
514 int n = in.capacity();
515 for (int i = 0; i < n; ++i) {
516 convertBasisToKGrid(in[i], workDft_);
517 fft().inverseTransformUnsafe(workDft_, out[i]);
518 }
519 }
520
521 /*
522 * Convert a single field from r-grid to basis format.
523 */
524 template <int D, class RFT, class KFT, class FFT>
526 RFT const & in,
527 DArray<double> & out,
528 bool checkSymmetry,
529 double epsilon) const
530 {
531 checkAllocateField(workDft_, mesh().dimensions());
532 fft().forwardTransform(in, workDft_);
533 convertKGridToBasis(workDft_, out, checkSymmetry, epsilon);
534 }
535
536 /*
537 * Convert an array of fields from r-grid to basis format.
538 */
539 template <int D, class RFT, class KFT, class FFT>
541 DArray< RFT > const & in,
542 DArray< DArray <double> > & out,
543 bool checkSymmetry,
544 double epsilon) const
545 {
546 // Inspect input and output containers
547 int nMonomer, nMonomerOut, capacity;
548 IntVec<D> dimensions;
549 inspectFields(in, nMonomer, dimensions);
550 UTIL_CHECK(mesh().dimensions() == dimensions);
551 inspectArrays(out, nMonomerOut, capacity);
552 UTIL_CHECK(nMonomer == nMonomerOut);
553 checkAllocateField(workDft_, dimensions);
554
555 // Convert RGrid -> KGrid -> Basis for each field
556 for (int i = 0; i < nMonomer; ++i) {
557 fft().forwardTransform(in[i], workDft_);
558 convertKGridToBasis(workDft_, out[i], checkSymmetry, epsilon);
559 }
560 }
561
562 /*
563 * Convert a field file from basis to r-grid format.
564 */
565 template <int D, class RFT, class KFT, class FFT>
567 std::string const & inFileName,
568 std::string const & outFileName) const
569 {
571 checkAllocateBasis(inFileName);
572 UnitCell<D> tmpUnitCell;
573
574 readFieldsBasis(inFileName, tmpFieldsBasis_, tmpUnitCell);
575 convertBasisToRGrid(tmpFieldsBasis_, tmpFieldsRGrid_);
576 writeFieldsRGrid(outFileName, tmpFieldsRGrid_, tmpUnitCell);
577 }
578
579 /*
580 * Convert a field file from r-grid to basis format.
581 */
582 template <int D, class RFT, class KFT, class FFT>
584 std::string const & inFileName,
585 std::string const & outFileName) const
586 {
588 checkAllocateBasis(inFileName);
589 UnitCell<D> tmpUnitCell;
590 readFieldsRGrid(inFileName, tmpFieldsRGrid_, tmpUnitCell);
591 convertRGridToBasis(tmpFieldsRGrid_, tmpFieldsBasis_);
592 writeFieldsBasis(outFileName, tmpFieldsBasis_, tmpUnitCell);
593 }
594
595 // FFT conversion functions (KGrid <-> RGrid) [Same in Rpc and Rpg]
596
597 /*
598 * Apply inverse FFT to an array of k-grid fields, converting to r-grid.
599 */
600 template <int D, class RFT, class KFT, class FFT>
602 DArray< KFT > const & in,
603 DArray< RFT >& out) const
604 {
605 UTIL_ASSERT(in.capacity() == out.capacity());
606 int n = in.capacity();
607 for (int i = 0; i < n; ++i) {
608 fft().inverseTransformSafe(in[i], out[i]);
609 }
610 }
611
612 /*
613 * Apply inverse FFT to a single k-grid field, converting to r-grid.
614 */
615 template <int D, class RFT, class KFT, class FFT>
617 KFT const & in, RFT& out) const
618 {
619 fft().inverseTransformSafe(in, out);
620 }
621
622 /*
623 * Apply forward FFT to an array of r-grid fields, converting to k-grid.
624 */
625 template <int D, class RFT, class KFT, class FFT>
627 DArray< RFT > const & in,
628 DArray< KFT >& out) const
629 {
630 UTIL_ASSERT(in.capacity() == out.capacity());
631 int n = in.capacity();
632 for (int i = 0; i < n; ++i) {
633 fft().forwardTransform(in[i], out[i]);
634 }
635 }
636
637 /*
638 * Apply forward FFT to a single r-grid field, converting to k-grid.
639 */
640 template <int D, class RFT, class KFT, class FFT>
642 RFT const & in,
643 KFT& out) const
644 { fft().forwardTransform(in, out); }
645
646 /*
647 * Convert a field file from k-grid to r-grid format.
648 */
649 template <int D, class RFT, class KFT, class FFT>
651 std::string const & inFileName,
652 std::string const & outFileName) const
653 {
656 UnitCell<D> tmpUnitCell;
657 readFieldsKGrid(inFileName, tmpFieldsKGrid_, tmpUnitCell);
658 for (int i = 0; i < nMonomer_; ++i) {
659 fft().inverseTransformUnsafe(tmpFieldsKGrid_[i],
660 tmpFieldsRGrid_[i]);
661 }
662 writeFieldsRGrid(outFileName, tmpFieldsRGrid_, tmpUnitCell);
663 }
664
665 /*
666 * Convert a field file from r-grid to k-grid format.
667 */
668 template <int D, class RFT, class KFT, class FFT>
670 std::string const & inFileName,
671 std::string const & outFileName) const
672 {
675 UnitCell<D> tmpUnitCell;
676 readFieldsRGrid(inFileName, tmpFieldsRGrid_, tmpUnitCell);
677 convertRGridToKGrid(tmpFieldsRGrid_, tmpFieldsKGrid_);
678 writeFieldsKGrid(outFileName, tmpFieldsKGrid_, tmpUnitCell);
679 }
680
681 // Field Inspection
682
683 /*
684 * Test if a single r-grid field has declared space group symmetry.
685 * Return true if symmetric, false otherwise. Print error values
686 * if verbose == true and hasSymmetry == false.
687 */
688 template <int D, class RFT, class KFT, class FFT>
690 RFT const & in,
691 double epsilon,
692 bool verbose) const
693 {
694 checkAllocateField(workDft_, mesh().dimensions());
695 fft().forwardTransform(in, workDft_);
696 return hasSymmetry(workDft_, epsilon, verbose);
697 }
698
699 /*
700 * Check if r-grid fields have declared space group symmetry.
701 */
702 template <int D, class RFT, class KFT, class FFT>
704 std::string const & inFileName,
705 double epsilon) const
706 {
708
709 // Read fields
710 UnitCell<D> tmpUnitCell;
711 readFieldsRGrid(inFileName, tmpFieldsRGrid_, tmpUnitCell);
712
713 // Check symmetry for all fields
714 for (int i = 0; i < nMonomer_; ++i) {
715 bool symmetric;
716 symmetric = hasSymmetry(tmpFieldsRGrid_[i], epsilon);
717 if (!symmetric) {
718 return false;
719 }
720 }
721 return true;
722 }
723
724 /*
725 * Compare two fields in basis format, write report to Log file.
726 */
727 template <int D, class RFT, class KFT, class FFT>
729 DArray< DArray<double> > const & field1,
730 DArray< DArray<double> > const & field2) const
731 {
732 BFieldComparison comparison(1);
733 comparison.compare(field1, field2);
734
735 Log::file() << "\n Basis expansion field comparison results"
736 << std::endl;
737 Log::file() << " Maximum Absolute Difference: "
738 << comparison.maxDiff() << std::endl;
739 Log::file() << " Root-Mean-Square Difference: "
740 << comparison.rmsDiff() << "\n" << std::endl;
741 }
742
743 template <int D, class RFT, class KFT, class FFT>
745 std::string const & filename1,
746 std::string const & filename2) const
747 {
748 DArray< DArray<double> > fields1, fields2;
749 UnitCell<D> tmpUnitCell;
750 // Unallocated arrays will be allocated in readFieldsBasis
751 readFieldsBasis(filename1, fields1, tmpUnitCell);
752 readFieldsBasis(filename2, fields2, tmpUnitCell);
753 compareFieldsBasis(fields1, fields2);
754 }
755
756 template <int D, class RFT, class KFT, class FFT>
758 std::string const & filename1,
759 std::string const & filename2) const
760 {
761 DArray< RFT > fields1, fields2;
762 UnitCell<D> tmpUnitCell;
763 // Unallocated arrays will be allocated in readFieldsRGrid
764 readFieldsRGrid(filename1, fields1, tmpUnitCell);
765 readFieldsRGrid(filename2, fields2, tmpUnitCell);
766 compareFieldsRGrid(fields1, fields2);
767 }
768
769 // Field Scaling
770
771 /*
772 * Multiply a single field in basis format by a constant factor.
773 */
774 template <int D, class RFT, class KFT, class FFT>
776 DArray<double> & field,
777 double factor) const
778 {
779 UTIL_CHECK(field.isAllocated());
780 int n = field.capacity();
781 for (int i = 0; i < n; ++i) {
782 field[i] *= factor;
783 }
784 }
785
786 /*
787 * Rescale an array of fields in basis format by a constant factor.
788 */
789 template <int D, class RFT, class KFT, class FFT>
791 DArray< DArray<double> >& fields,
792 double factor) const
793 {
794 int n = fields.capacity();
795 UTIL_CHECK(n > 0);
796 int m = fields[0].capacity();
797 for (int i = 0; i < n; ++i) {
798 UTIL_CHECK(fields[i].capacity() == m);
799 scaleFieldBasis(fields[i], factor);
800 }
801 }
802
803 /*
804 * Rescale fields in files by a constant factor.
805 */
806 template <int D, class RFT, class KFT, class FFT>
808 std::string const & inFileName,
809 std::string const & outFileName,
810 double factor) const
811 {
812 checkAllocateBasis(inFileName);
813 UnitCell<D> tmpUnitCell;
814 readFieldsBasis(inFileName, tmpFieldsBasis_, tmpUnitCell);
815 scaleFieldsBasis(tmpFieldsBasis_, factor);
816 writeFieldsBasis(outFileName, tmpFieldsBasis_, tmpUnitCell);
817 }
818
819 /*
820 * Rescale fields in r-grid format by a constant factor.
821 */
822 template <int D, class RFT, class KFT, class FFT>
824 DArray< RFT > & fields,
825 double factor) const
826 {
827 int n = fields.capacity();
828 for (int i = 0; i < n; ++i) {
829 scaleFieldRGrid(fields[i], factor);
830 }
831 }
832
833 /*
834 * Rescale fields by a constant factor, read and write to file.
835 */
836 template <int D, class RFT, class KFT, class FFT>
838 std::string const & inFileName,
839 std::string const & outFileName,
840 double factor) const
841 {
843 UnitCell<D> tmpUnitCell;
844 bool isSymmetric;
845 isSymmetric = readFieldsRGrid(inFileName, tmpFieldsRGrid_,
846 tmpUnitCell);
847 scaleFieldsRGrid(tmpFieldsRGrid_, factor);
848 writeFieldsRGrid(outFileName, tmpFieldsRGrid_, tmpUnitCell,
849 isSymmetric);
850 }
851
852 // Computation of Estimated W Fields
853
854 /*
855 * Convert an array of c fields to estimated w fields, in basis form.
856 */
857 template <int D, class RFT, class KFT, class FFT>
859 DMatrix<double> const & chi,
860 DArray< DArray<double> > & fields) const
861 {
862 // Check dimensions
863 UTIL_CHECK(fields.capacity() == nMonomer_);
864 UTIL_CHECK(chi.capacity1() == nMonomer_);
865 UTIL_CHECK(chi.capacity2() == nMonomer_);
866 const int nb = fields[0].capacity();
867 for (int i = 0; i < nMonomer_; ++i) {
868 UTIL_CHECK(fields[i].capacity() == nb);
869 }
870
871 // Allocate small work array (one element per monomer type)
872 DArray<double> wtmp;
873 wtmp.allocate(nMonomer_);
874
875 // Compute estimated w fields from c fields - convert in place
876 int i, j, k;
877 for (i = 0; i < nb; ++i) {
878 for (j = 0; j < nMonomer_; ++j) {
879 wtmp[j] = 0.0;
880 for (k = 0; k < nMonomer_; ++k) {
881 wtmp[j] += chi(j,k)*fields[k][i];
882 }
883 }
884 for (j = 0; j < nMonomer_; ++j) {
885 fields[j][i] = wtmp[j];
886 }
887 }
888
889 }
890
891 /*
892 * Convert a file of c fields to estimated w fields, in basis format.
893 */
894 template <int D, class RFT, class KFT, class FFT>
896 std::string const & inFileName,
897 std::string const & outFileName,
898 DMatrix<double> const & chi) const
899 {
900 // Preconditions
901 UTIL_CHECK(chi.capacity1() == nMonomer_);;
902 UTIL_CHECK(chi.capacity2() == nMonomer_);;
903
904 // Set and check nb = number of basis functions
905 checkAllocateBasis(inFileName);
906 UTIL_CHECK(basis().isInitialized());
907 const int nb = basis().nBasis();
908 UTIL_CHECK(nb > 0);
909
910 // Open input file
911 std::ifstream in;
912 fileMaster().openInputFile(inFileName, in);
913
914 // Read input file header
915 int nMonomerIn;
916 bool isSymmetric;
917 UnitCell<D> tmpUnitCell;
918 readFieldHeader(in, nMonomerIn, tmpUnitCell, isSymmetric);
919 int nBasisIn = readNBasis(in);
920 UTIL_CHECK(nMonomer_ == nMonomerIn);
921 UTIL_CHECK(isSymmetric);
922
923 // Read c field data into tmpFieldsBasis_
924 Prdc::readBasisData(in, tmpFieldsBasis_,
925 tmpUnitCell, mesh(), basis(), nBasisIn);
926 in.close();
927
928 // Compute estimated w, stored in tmpFilesBasis_
929 estimateWBasis(chi, tmpFieldsBasis_);
930
931 // Open and write output file
932 std::ofstream out;
933 fileMaster().openOutputFile(outFileName, out);
934 writeFieldsBasis(out, tmpFieldsBasis_, tmpUnitCell);
935 }
936
937
938 // Grid manipulation utilities
939
940 /*
941 * Replicate unit cell a specified number of times in each direction.
942 */
943 template <int D, class RFT, class KFT, class FFT>
945 std::string filename,
946 DArray<RFT> const & fields,
947 UnitCell<D> const & unitCell,
948 IntVec<D> const & replicas) const
949 {
950 std::ofstream file;
951 fileMaster().openOutputFile(filename, file);
952 replicateUnitCell(file, fields, unitCell, replicas);
953 file.close();
954 }
955
956 /*
957 * Replicate unit cell a specified number of times in each direction.
958 */
959 template <int D, class RFT, class KFT, class FFT>
961 std::string const & inFileName,
962 std::string const & outFileName,
963 IntVec<D> const & replicas) const
964 {
966 UnitCell<D> tmpUnitCell;
967 readFieldsRGrid(inFileName, tmpFieldsRGrid_, tmpUnitCell);
968 replicateUnitCell(outFileName, tmpFieldsRGrid_, tmpUnitCell,
969 replicas);
970 }
971
972 /*
973 * Expand the number of spatial dimensions of an RField.
974 */
975 template <int D, class RFT, class KFT, class FFT>
977 std::string filename,
978 DArray<RFT> const & fields,
979 UnitCell<D> const & unitCell, int d,
980 DArray<int> newGridDimensions) const
981 {
982 std::ofstream file;
983 fileMaster().openOutputFile(filename, file);
984 expandRGridDimension(file, fields, unitCell, d, newGridDimensions);
985 file.close();
986 }
987
988 /*
989 * Expand the number of spatial dimensions of an RField.
990 */
991 template <int D, class RFT, class KFT, class FFT>
993 std::string const & inFileName,
994 std::string const & outFileName,
995 int d,
996 DArray<int> newGridDimensions) const
997 {
999 UnitCell<D> tmpUnitCell;
1000 readFieldsRGrid(inFileName, tmpFieldsRGrid_, tmpUnitCell);
1001 expandRGridDimension(outFileName, tmpFieldsRGrid_, tmpUnitCell,
1002 d, newGridDimensions);
1003 }
1004
1005 // File Header IO Utilities
1006
1007 /*
1008 * Read common part of field header.
1009 *
1010 * Extracts number of monomers (i.e., number of fields) and unitCell
1011 * data (lattice type and parameters) from the file, and returns these
1012 * via non-const reference parameters nMonomer and unitCell.
1013 *
1014 * Also validates lattice type and groupName (if any), and constructs
1015 * a symmetry-adapted basis if there is a group but no basis.
1016 */
1017 template <int D, class RFT, class KFT, class FFT>
1019 std::istream& in,
1020 int& nMonomer,
1021 UnitCell<D>& unitCell,
1022 bool & isSymmetric) const
1023 {
1024 // Preconditions
1025 UTIL_CHECK(latticePtr_);
1026 if (unitCell.lattice() == UnitCell<D>::Null) {
1027 UTIL_CHECK(unitCell.nParameter() == 0);
1028 } else {
1029 UTIL_CHECK(unitCell.nParameter() > 0);
1030 UTIL_CHECK(unitCell.lattice() == lattice());
1031 }
1032
1033 // Read field header to set unitCell, groupNameIn, nMonomer
1034 int ver1, ver2;
1035 std::string groupNameIn;
1036
1037 Pscf::Prdc::readFieldHeader(in, ver1, ver2, unitCell,
1038 groupNameIn, nMonomer);
1039 // Note: Function definition in prdc/crystal/fieldHeader.tpp
1040
1041 // Checks of data from header
1042 UTIL_CHECK(ver1 == 1);
1043 //UTIL_CHECK(ver2 == 0);
1044 UTIL_CHECK(unitCell.isInitialized());
1045 UTIL_CHECK(unitCell.lattice() != UnitCell<D>::Null);
1046 UTIL_CHECK(unitCell.nParameter() > 0);
1047
1048 // Validate or initialize lattice type
1049 if (lattice() != unitCell.lattice()) {
1050 Log::file() << std::endl
1051 << "Error - Mismatched lattice types "
1052 << "in FieldIo function readFieldHeader:\n"
1053 << " FieldIo::lattice :" << lattice() << "\n"
1054 << " Unit cell lattice :" << unitCell.lattice()
1055 << "\n";
1056 UTIL_THROW("Mismatched lattice types");
1057 }
1058
1059 // Check for presence of group name in the field file header
1060 isSymmetric = false;
1061 if (groupNameIn != "") {
1062 isSymmetric = true;
1063 }
1064
1065 // Process group and basis (if any)
1066 if (hasGroup()) {
1067
1068 // Check consistency of groupName values
1069 if (isSymmetric) {
1070 UTIL_CHECK(groupNamePtr_);
1071 if (groupNameIn != groupName()) {
1072 Log::file() << std::endl
1073 << "Error - Mismatched group names in "
1074 << "FieldIo member function readFieldHeader:\n"
1075 << " FieldIo::groupName :" << groupName() << "\n"
1076 << " Field file header :" << groupNameIn << "\n";
1077 UTIL_THROW("Mismatched group names");
1078 }
1079 }
1080
1081 // If there is a group but no basis, construct a basis
1082 UTIL_CHECK(basisPtr_);
1083 if (!basis().isInitialized()) {
1084 basisPtr_->makeBasis(mesh(), unitCell, group());
1085 }
1086 UTIL_CHECK(basis().isInitialized());
1087
1088 } else {
1089
1090 if (isSymmetric) {
1091 Log::file() << std::endl
1092 << "Warning: Group name found in a field file header"
1093 << "but no group declared in the parameter file.\n";
1094 }
1095
1096 }
1097
1098 }
1099
1100 /*
1101 * Write a field file header.
1102 */
1103 template <int D, class RFT, class KFT, class FFT>
1105 std::ostream &out,
1106 int nMonomer,
1107 UnitCell<D> const & unitCell,
1108 bool isSymmetric) const
1109 {
1110 int v1 = 1;
1111 int v2 = 0;
1112 std::string gName = "";
1113 if (isSymmetric) {
1115 gName = groupName();
1116 }
1117 Pscf::Prdc::writeFieldHeader(out, v1, v2, unitCell,
1118 gName, nMonomer);
1119 // Note: This function is defined in prdc/crystal/fieldHeader.tpp
1120 }
1121
1122 // Protected functions to check and allocate private workspace arrays
1123
1124 /*
1125 * If necessary, allocate r-grid workspace.
1126 */
1127 template <int D, class RFT, class KFT, class FFT>
1129 {
1130 if (isAllocatedRGrid_) return;
1131
1132 UTIL_CHECK(nMonomer_ > 0);
1133 UTIL_CHECK(mesh().size() > 0);
1134 IntVec<D> const & meshDimensions = mesh().dimensions();
1135 tmpFieldsRGrid_.allocate(nMonomer_);
1136 for (int i = 0; i < nMonomer_; ++i) {
1137 tmpFieldsRGrid_[i].allocate(meshDimensions);
1138 }
1139 isAllocatedRGrid_ = true;
1140 }
1141
1142 /*
1143 * If necessary, allocate k-grid field workspace.
1144 */
1145 template <int D, class RFT, class KFT, class FFT>
1147 {
1148 if (isAllocatedKGrid_) return;
1149
1150 UTIL_CHECK(nMonomer_ > 0);
1151 UTIL_CHECK(mesh().size() > 0);
1152 IntVec<D> const & meshDimensions = mesh().dimensions();
1153 tmpFieldsKGrid_.allocate(nMonomer_);
1154 for (int i = 0; i < nMonomer_; ++i) {
1155 tmpFieldsKGrid_[i].allocate(meshDimensions);
1156 }
1157 isAllocatedKGrid_ = true;
1158 }
1159
1160 /*
1161 * If necessary, allocate basis field workspace.
1162 */
1163 template <int D, class RFT, class KFT, class FFT>
1164 void
1166 std::string const & inFileName) const
1167 {
1168 if (isAllocatedBasis_) return;
1169
1170 UTIL_CHECK(nMonomer_ > 0);
1172 if (!basis().isInitialized()) {
1173 // Peek at field header to initialize basis
1174 std::ifstream file;
1175 fileMaster().openInputFile(inFileName, file);
1176 int nMonomerIn;
1177 bool isSymmetricIn;
1178 UnitCell<D> tmpUnitCell;
1179 readFieldHeader(file, nMonomerIn, tmpUnitCell, isSymmetricIn);
1180 // Note: readFieldHeader can initialize basis if needed
1181 file.close();
1182 }
1183 UTIL_CHECK(basis().isInitialized());
1184 int nBasis = basis().nBasis();
1185 UTIL_CHECK(nBasis > 0);
1186 tmpFieldsBasis_.allocate(nMonomer_);
1187 for (int i = 0; i < nMonomer_; ++i) {
1188 DArray<double>& field = tmpFieldsBasis_[i];
1189 field.allocate(nBasis);
1190 for (int j = 0; j < nBasis; ++j) {
1191 field[j] = 0.0;
1192 }
1193 }
1194 isAllocatedBasis_ = true;
1195 }
1196
1197} // namespace Prdc
1198} // namespace Pscf
1199#endif
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.
Definition IntVec.h:27
Description of a regular grid of points in a periodic domain.
Definition Mesh.h:61
Comparator for fields in symmetry-adapted basis format.
Symmetry-adapted Fourier basis for pseudo-spectral SCFT.
Definition Basis.h:383
Fourier transform wrapper.
Definition cpu/FFT.h:38
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.
Definition SpaceGroup.h:32
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.
Definition UnitCell.h:56
int capacity() const
Return allocated size.
Definition Array.h:159
Dynamically allocatable contiguous array template.
Definition DArray.h:32
void allocate(int capacity)
Allocate the underlying C array.
Definition DArray.h:199
bool isAllocated() const
Return true if this DArray has been allocated, false otherwise.
Definition DArray.h:247
Dynamically allocated Matrix.
Definition DMatrix.h:25
A FileMaster manages input and output files for a simulation.
Definition FileMaster.h:143
static std::ostream & file()
Get log ostream by reference.
Definition Log.cpp:59
int capacity2() const
Get number of columns (range of the second array index).
Definition Matrix.h:143
int capacity1() const
Get number of rows (range of the first array index).
Definition Matrix.h:136
#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
#define UTIL_ASSERT(condition)
Assertion macro suitable for debugging serial or parallel code.
Definition global.h:75
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.
Definition CField.cpp:11
PSCF package top-level namespace.
Definition param_pc.dox:1