PSCF v1.1
pspc/field/FieldIo.tpp
1#ifndef PSPC_FIELD_IO_TPP
2#define PSPC_FIELD_IO_TPP
3
4/*
5* PSCF - Polymer Self-Consistent Field Theory
6*
7* Copyright 2016 - 2022, The Regents of the University of Minnesota
8* Distributed under the terms of the GNU General Public License.
9*/
10
11#include "FieldIo.h"
12
13#include <pscf/crystal/shiftToMinimum.h>
14#include <pscf/mesh/MeshIterator.h>
15#include <pscf/math/IntVec.h>
16
17#include <util/misc/Log.h>
18#include <util/format/Str.h>
19#include <util/format/Int.h>
20#include <util/format/Dbl.h>
21
22#include <iomanip>
23#include <string>
24
25namespace Pscf {
26namespace Pspc
27{
28
29 using namespace Util;
30
31 /*
32 * Constructor.
33 */
34 template <int D>
35 FieldIo<D>::FieldIo()
36 : meshPtr_(0),
37 fftPtr_(0),
38 groupNamePtr_(0),
39 groupPtr_(0),
40 basisPtr_(0),
41 fileMasterPtr_()
42 {}
43
44 /*
45 * Destructor.
46 */
47 template <int D>
48 FieldIo<D>::~FieldIo()
49 {}
50
51 /*
52 * Get and store addresses of associated objects.
53 */
54 template <int D>
55 void
56 FieldIo<D>::associate(Mesh<D> const & mesh,
57 FFT<D> const & fft,
58 typename UnitCell<D>::LatticeSystem & lattice,
59 std::string & groupName,
60 SpaceGroup<D> & group,
61 Basis<D> & basis,
62 FileMaster const & fileMaster)
63 {
64 meshPtr_ = &mesh;
65 fftPtr_ = &fft;
66 latticePtr_ = &lattice;
67 groupNamePtr_ = &groupName;
68 groupPtr_ = &group;
69 basisPtr_ = &basis;
70 fileMasterPtr_ = &fileMaster;
71 }
72
73 template <int D>
74 void FieldIo<D>::readFieldBasis(std::istream& in, DArray<double>& field,
75 UnitCell<D>& unitCell)
76 const
77 {
78 DArray<DArray<double> > fields;
79 if (field.isAllocated()) {
80 fields.allocate(1);
81 fields[0].allocate(field.capacity());
82 } // otherwise pass unallocated field into readFieldsBasis
83
84 readFieldsBasis(in, fields, unitCell);
85 UTIL_CHECK(fields.capacity() == 1); // Check that it only read 1 field
86 field = fields[0];
87 }
88
89 template <int D>
90 void FieldIo<D>::readFieldBasis(std::string filename,
91 DArray<double>& field,
92 UnitCell<D>& unitCell)
93 const
94 {
95 DArray<DArray<double> > fields;
96 if (field.isAllocated()) {
97 fields.allocate(1);
98 fields[0].allocate(field.capacity());
99 } // otherwise pass unallocated field into readFieldsBasis
100 readFieldsBasis(filename, fields, unitCell);
101 UTIL_CHECK(fields.capacity() == 1); // Check that it only read 1 field
102 field = fields[0];
103 }
104
105 template <int D>
106 void FieldIo<D>::writeFieldBasis(std::ostream& out,
107 DArray<double> const & field,
108 UnitCell<D> const & unitCell)
109 const
110 {
111 DArray<DArray<double> > fields;
112 fields.allocate(1);
113 fields[0].allocate(field.capacity());
114 fields[0] = field;
115 writeFieldsBasis(out, fields, unitCell);
116 }
117
118 template <int D>
119 void FieldIo<D>::writeFieldBasis(std::string filename,
120 DArray<double> const & field,
121 UnitCell<D> const & unitCell)
122 const
123 {
124 DArray<DArray<double> > fields;
125 fields.allocate(1);
126 fields[0].allocate(field.capacity());
127 fields[0] = field;
128 writeFieldsBasis(filename, fields, unitCell);
129 }
130
131 template <int D>
132 void FieldIo<D>::readFieldsBasis(std::istream& in,
133 DArray< DArray<double> >& fields,
134 UnitCell<D>& unitCell)
135 const
136 {
137 int nMonomer;
138 FieldIo<D>::readFieldHeader(in, nMonomer, unitCell);
139 UTIL_CHECK(basis().isInitialized());
140
141 // Read the number of stars into nStarIn
142 std::string label;
143 in >> label;
144 if (label != "N_star" && label != "N_basis") {
145 std::string msg = "\n";
146 msg += "Error reading field file:\n";
147 msg += "Expected N_basis or N_star, but found [";
148 msg += label;
149 msg += "]";
150 UTIL_THROW(msg.c_str());
151 }
152 int nStarIn;
153 in >> nStarIn;
154 UTIL_CHECK(nStarIn > 0);
155
156 // If "fields" parameter is allocated, check if dimensions
157 // match those of the system's mesh. Otherwise, allocate.
158
159 int fieldCapacity;
160 if (fields.isAllocated()) {
161 // If outer DArray is allocated, require that it matches the
162 // number of inputted fields and that internal DArrays are also
163 // allocated all with the same dimensions.
164 int nMonomerFields = fields.capacity();
165
166 UTIL_CHECK(nMonomerFields > 0);
167 UTIL_CHECK(nMonomerFields == nMonomer);
168
169 // Check that all internal DArrays have same dimension.
170 fieldCapacity = fields[0].capacity();
171 for (int i = 0; i < nMonomer; ++i) {
172 UTIL_CHECK( fields[i].capacity() == fieldCapacity );
173 }
174 } else {
175 // Else, allocate fields to the number of inputted fields
176 // and the internal dimensions to the number of stars.
177 fields.allocate(nMonomer);
178 fieldCapacity = nStarIn;
179
180 for (int i = 0; i < nMonomer; ++i) {
181 fields[i].allocate(fieldCapacity);
182 }
183 }
184
185 // Initialize all field array elements to zero
186 for (int i = 0; i < nMonomer; ++i) {
187 for (int j = 0; j < fieldCapacity; ++j) {
188 fields[i][j] = 0.0;
189 }
190 }
191
192 // Reset nStarIn = min(nStarIn, fieldCapacity)
193 if (fieldCapacity < nStarIn) {
194 nStarIn = fieldCapacity;
195 }
196
197 // Allocate temp arrays used to read in components
198 DArray<double> temp, temp2;
199 temp.allocate(nMonomer);
200 temp2.allocate(nMonomer);
201
202 typename Basis<D>::Star const * starPtr;
203 typename Basis<D>::Star const * starPtr2;
204 IntVec<D> waveIn, waveIn2;
205 int sizeIn, sizeIn2;
206 int starId, starId2;
207 int basisId, basisId2;
208 int waveId, waveId2;
209
210 std::complex<double> coeff, phasor;
211 IntVec<D> waveBz, waveDft;
212 int nReversedPair = 0;
213 bool waveExists, sizeMatches;
214
215 // Loop over stars in input file to read field components
216 int i = 0;
217 while (i < nStarIn) {
218
219 // Read next line of data
220 for (int j = 0; j < nMonomer; ++j) {
221 in >> temp[j]; // field components
222 }
223 in >> waveIn; // wave of star
224 in >> sizeIn; // # of waves in star
225 ++i;
226
227 sizeMatches = false;
228 waveExists = false;
229
230 // Check if waveIn is in first Brillouin zone (FBZ) for the mesh.
231 waveBz = shiftToMinimum(waveIn, mesh().dimensions(), unitCell);
232 waveExists = (waveIn == waveBz);
233
234 if (waveExists) {
235
236 // Find the star containing waveIn
237 waveDft = waveIn;
238 mesh().shift(waveDft);
239 waveId = basis().waveId(waveDft);
240 starId = basis().wave(waveId).starId;
241 starPtr = &basis().star(starId);
242 UTIL_CHECK(!(starPtr->cancel));
243 basisId = starPtr->basisId;
244
245 if (starPtr->size == sizeIn) {
246 sizeMatches = true;
247 } else {
248 Log::file()
249 << "Warning: Inconsistent star size (line ignored)\n"
250 << "wave from file = " << waveIn << "\n"
251 << "size from file = " << sizeIn << "\n"
252 << "size of star = " << starPtr->size
253 << "\n";
254 sizeMatches = false;
255 }
256
257 }
258
259 if (waveExists && sizeMatches) { // Attempt to process wave
260
261 if (starPtr->invertFlag == 0) {
262
263 if (starPtr->waveBz == waveIn) {
264
265 // Copy components of closed star to fields array
266 for (int j = 0; j < nMonomer; ++j) {
267 fields[j][basisId] = temp[j];
268 }
269
270 } else {
271
272 Log::file()
273 << "Inconsistent wave of closed star on input\n"
274 << "wave from file = " << waveIn << "\n"
275 << "starId of wave = " << starId << "\n"
276 << "waveBz of star = " << starPtr->waveBz
277 << "\n";
278
279 }
280
281 } else {
282
283 // Read the next line
284 for (int j = 0; j < nMonomer; ++j) {
285 in >> temp2[j]; // components of field
286 }
287 in >> waveIn2; // wave of star
288 in >> sizeIn2; // # of wavevectors in star
289 ++i;
290
291 // Check that waveIn2 is also in the 1st BZ
292 waveBz =
293 shiftToMinimum(waveIn2, mesh().dimensions(), unitCell);
294 UTIL_CHECK(waveIn2 == waveBz);
295
296 // Identify the star containing waveIn2
297 waveDft = waveIn2;
298 mesh().shift(waveDft);
299 waveId2 = basis().waveId(waveDft);
300 starId2 = basis().wave(waveId2).starId;
301 starPtr2 = &basis().star(starId2);
302 UTIL_CHECK(!(starPtr2->cancel));
303 basisId2 = starPtr2->basisId;
304 UTIL_CHECK(starPtr2->size == sizeIn2);
305 UTIL_CHECK(sizeIn == sizeIn2);
306
307 if (starPtr->invertFlag == 1) {
308
309 // This is a pair of open stars written in the same
310 // order as in this basis. Check preconditions:
311 UTIL_CHECK(starPtr2->invertFlag == -1);
312 UTIL_CHECK(starId2 = starId + 1);
313 UTIL_CHECK(basisId2 = basisId + 1);
314 UTIL_CHECK(starPtr->waveBz == waveIn);
315 UTIL_CHECK(starPtr2->waveBz == waveIn2);
316
317 // Copy components for both stars into fields array
318 for (int j = 0; j < nMonomer; ++j) {
319 fields[j][basisId] = temp[j];
320 fields[j][basisId2] = temp2[j];
321 }
322
323 } else
324 if (starPtr->invertFlag == -1) {
325
326 // This is a pair of open stars written in opposite
327 // order from in this basis. Check preconditions:
328 UTIL_CHECK(starPtr2->invertFlag == 1);
329 UTIL_CHECK(starId == starId2 + 1);
330 UTIL_CHECK(basisId == basisId2 + 1);
331 UTIL_CHECK(waveId == starPtr->beginId);
332
333 // Check that waveIn2 is negation of waveIn
334 IntVec<D> nVec;
335 nVec.negate(waveIn);
336 nVec =
337 shiftToMinimum(nVec, mesh().dimensions(), unitCell);
338 UTIL_CHECK(waveIn2 == nVec);
339
340 /*
341 * Consider two related stars, C and D, that are listed
342 * in the order (C,D) in the basis used in this code (the
343 * reading program), but that were listed in the opposite
344 * order (D,C) in the program that wrote the file (the
345 * writing program). In the basis of the reading program,
346 * star C has star index starId2, while star D has index
347 * starId = starid2 + 1.
348 *
349 * Let f(r) and f^{*}(r) denote the basis functions used
350 * by the reading program for stars C and D, respectively.
351 * Let u(r) and u^{*}(r) denote the corresponding basis
352 * functions used by the writing program for stars C
353 * and D. Let exp(i phi) denote the unit magnitude
354 * coefficient (i.e., phasor) within f(r) of the wave
355 * with wave index waveId2, which was the characteristic
356 * wave for star C in the writing program. The
357 * coefficient of this wave within the basis function
358 * u(r) used by the writing program must instead be real
359 * and positive. This implies that
360 * u(r) = exp(-i phi) f(r).
361 *
362 * Focus on the contribution to the field for a specific
363 * monomer type j. Let a and b denote the desired
364 * coefficients of stars C and D in the reading program,
365 * for which the total contribution of both stars to the
366 * field is:
367 *
368 * (a - ib) f(r) + (a + ib) f^{*}(r)
369 *
370 * Let A = temp[j] and B = temp2[j] denote the
371 * coefficients read from file in order (A,B). Noting
372 * that the stars were listed in the order (D,C) in the
373 * basis used by the writing program, the contribution
374 * of both stars must be (A-iB)u^{*}(r)+(A+iB)u(r), or:
375 *
376 * (A+iB) exp(-i phi)f(r) + (A-iB) exp(i phi) f^{*}(r)
377 *
378 * Comparing coefficients of f^{*}(r), we find that
379 *
380 * (a + ib) = (A - iB) exp(i phi)
381 *
382 * This equality is implemented below, where the
383 * variable "phasor" is set equal to exp(i phi).
384 */
385 phasor = basis().wave(waveId2).coeff;
386 phasor = phasor/std::abs(phasor);
387 for (int j = 0; j < nMonomer; ++j) {
388 coeff = std::complex<double>(temp[j],-temp2[j]);
389 coeff *= phasor;
390 fields[j][basisId2] = real(coeff);
391 fields[j][basisId ] = imag(coeff);
392 }
393
394 // Increment count of number of reversed open pairs
395 ++nReversedPair;
396
397 } else {
398 UTIL_THROW("Invalid starInvert value");
399 }
400
401 } // if (wavePtr->invertFlag == 0) ... else ...
402
403 } // if (waveExists && sizeMatches)
404
405 } // end while (i < nStarIn)
406
407 if (nReversedPair > 0) {
408 Log::file() << "\n";
409 Log::file() << nReversedPair << " reversed pairs of open stars"
410 << " detected in FieldIo::readFieldsBasis\n";
411 }
412
413 }
414
415 template <int D>
416 void FieldIo<D>::readFieldsBasis(std::string filename,
417 DArray<DArray<double> >& fields,
418 UnitCell<D>& unitCell)
419 const
420 {
421
422 std::ifstream file;
423 fileMaster().openInputFile(filename, file);
424 readFieldsBasis(file, fields, unitCell);
425 file.close();
426 }
427
428 template <int D>
429 void
430 FieldIo<D>::writeFieldsBasis(std::ostream &out,
431 DArray<DArray<double> > const & fields,
432 UnitCell<D> const & unitCell)
433 const
434 {
435 int nMonomer = fields.capacity();
436 UTIL_CHECK(nMonomer > 0);
437 UTIL_CHECK(basis().isInitialized());
438
439 // Write header
440 writeFieldHeader(out, nMonomer, unitCell);
441 int nStar = basis().nStar();
442 int nBasis = basis().nBasis();
443 out << "N_basis " << std::endl
444 << " " << nBasis << std::endl;
445
446 // Write fields
447 int ib = 0;
448 for (int i = 0; i < nStar; ++i) {
449 if (!basis().star(i).cancel) {
450 for (int j = 0; j < nMonomer; ++j) {
451 out << Dbl(fields[j][ib], 20, 10);
452 }
453 out << " ";
454 for (int j = 0; j < D; ++j) {
455 out << Int(basis().star(i).waveBz[j], 5);
456 }
457 out << Int(basis().star(i).size, 5) << std::endl;
458 ++ib;
459 }
460 }
461
462 }
463
464 template <int D>
465 void
466 FieldIo<D>::writeFieldsBasis(std::string filename,
467 DArray<DArray<double> > const & fields,
468 UnitCell<D> const & unitCell)
469 const
470 {
471 std::ofstream file;
472 fileMaster().openOutputFile(filename, file);
473 writeFieldsBasis(file, fields, unitCell);
474 file.close();
475 }
476
477 template <int D>
478 void FieldIo<D>::readFieldsRGrid(std::istream &in,
479 DArray<RField<D> >& fields,
480 UnitCell<D>& unitCell)
481 const
482 {
483 int nMonomer;
484 FieldIo<D>::readFieldHeader(in, nMonomer, unitCell);
485
486 // If "fields" parameter is allocated, check if dimensions match
487 // those of the system's mesh. Otherwise, allocate.
488
489 if (fields.isAllocated()) {
490 int nMonomerFields = fields.capacity();
491
492 UTIL_CHECK(nMonomerFields > 0);
493 UTIL_CHECK(nMonomerFields == nMonomer);
494
495 for (int i = 0; i < nMonomer; ++i) {
496 UTIL_CHECK(fields[i].meshDimensions() == mesh().dimensions());
497 }
498 } else {
499 fields.allocate(nMonomer);
500 for (int i = 0; i < nMonomer; ++i) {
501 fields[i].allocate(mesh().dimensions());
502 }
503 }
504
505 // Read and check input stream mesh dimensions
506 std::string label;
507 in >> label;
508 if (label != "mesh" && label != "ngrid") {
509 std::string msg = "\n";
510 msg += "Error reading field file:\n";
511 msg += "Expected mesh or ngrid, but found [";
512 msg += label;
513 msg += "]";
514 UTIL_THROW(msg.c_str());
515 }
516 IntVec<D> nGrid;
517 in >> nGrid;
518 UTIL_CHECK(nGrid == mesh().dimensions());
519
520 // Setup temporary workspace array.
521 DArray<RField<D> > temp;
522 temp.allocate(nMonomer);
523 for (int i = 0; i < nMonomer; ++i) {
524 temp[i].allocate(mesh().dimensions());
525 }
526
527 // Read Fields;
528 MeshIterator<D> itr(mesh().dimensions());
529 for (itr.begin(); !itr.atEnd(); ++itr) {
530 for (int i = 0; i < nMonomer; ++i) {
531 in >> std::setprecision(15) >> temp[i][itr.rank()];
532 }
533 }
534
535 int p = 0;
536 int q = 0;
537 int r = 0;
538 int s = 0;
539 int n1 = 0;
540 int n2 = 0;
541 int n3 = 0;
542
543 if (D==3) {
544 while (n1 < mesh().dimension(0)) {
545 q = p;
546 n2 = 0;
547 while (n2 < mesh().dimension(1)) {
548 r = q;
549 n3 = 0;
550 while (n3 < mesh().dimension(2)) {
551 for (int i = 0; i < nMonomer; ++i) {
552 fields[i][s] = temp[i][r];
553 }
554 r = r + (mesh().dimension(0) * mesh().dimension(1));
555 ++s;
556 ++n3;
557 }
558 q = q + mesh().dimension(0);
559 ++n2;
560 }
561 ++n1;
562 ++p;
563 }
564 }
565
566 else if (D==2) {
567 while (n1 < mesh().dimension(0)) {
568 r =q;
569 n2 = 0;
570 while (n2 < mesh().dimension(1)) {
571 for (int i = 0; i < nMonomer; ++i) {
572 fields[i][s] = temp[i][r];
573 }
574 r = r + (mesh().dimension(0));
575 ++s;
576 ++n2;
577 }
578 ++q;
579 ++n1;
580 }
581 }
582
583 else if (D==1) {
584
585 while (n1 < mesh().dimension(0)) {
586 for (int i = 0; i < nMonomer; ++i) {
587 fields[i][s] = temp[i][r];
588 }
589 ++r;
590 ++s;
591 ++n1;
592 }
593 }
594
595 else{
596 Log::file() << "Invalid Dimensions";
597 }
598
599 }
600
601 template <int D>
602 void FieldIo<D>::readFieldsRGrid(std::string filename,
603 DArray< RField<D> >& fields,
604 UnitCell<D>& unitCell)
605 const
606 {
607 std::ifstream file;
608 fileMaster().openInputFile(filename, file);
609 readFieldsRGrid(file, fields, unitCell);
610 file.close();
611 }
612
613 template <int D>
614 void FieldIo<D>::writeFieldsRGrid(std::ostream &out,
615 DArray<RField<D> > const & fields,
616 UnitCell<D> const & unitCell)
617 const
618 {
619 int nMonomer = fields.capacity();
620 UTIL_CHECK(nMonomer > 0);
621
622 writeFieldHeader(out, nMonomer, unitCell);
623 out << "mesh " << std::endl
624 << " " << mesh().dimensions() << std::endl;
625
626 DArray<RField<D> > temp;
627 temp.allocate(nMonomer);
628 for (int i = 0; i < nMonomer; ++i) {
629 temp[i].allocate(mesh().dimensions());
630 }
631
632 int p = 0;
633 int q = 0;
634 int r = 0;
635 int s = 0;
636 int n1 =0;
637 int n2 =0;
638 int n3 =0;
639
640 if (D==3) {
641 while (n3 < mesh().dimension(2)) {
642 q = p;
643 n2 = 0;
644 while (n2 < mesh().dimension(1)) {
645 r =q;
646 n1 = 0;
647 while (n1 < mesh().dimension(0)) {
648 for (int i = 0; i < nMonomer; ++i) {
649 temp[i][s] = fields[i][r];
650 }
651 r = r + (mesh().dimension(1) * mesh().dimension(2));
652 ++s;
653 ++n1;
654 }
655 q = q + mesh().dimension(2);
656 ++n2;
657 }
658 ++n3;
659 ++p;
660 }
661 }
662 else if (D==2) {
663 while (n2 < mesh().dimension(1)) {
664 r =q;
665 n1 = 0;
666 while (n1 < mesh().dimension(0)) {
667 for (int i = 0; i < nMonomer; ++i) {
668 temp[i][s] = fields[i][r];
669 }
670 r = r + (mesh().dimension(1));
671 ++s;
672 ++n1;
673 }
674 ++q;
675 ++n2;
676 }
677 }
678 else if (D==1) {
679 while (n1 < mesh().dimension(0)) {
680 for (int i = 0; i < nMonomer; ++i) {
681 temp[i][s] = fields[i][r];
682 }
683 ++r;
684 ++s;
685 ++n1;
686 }
687 } else {
688 Log::file() << "Invalid Dimensions";
689 }
690
691 // Write fields
692 MeshIterator<D> itr(mesh().dimensions());
693 for (itr.begin(); !itr.atEnd(); ++itr) {
694 for (int j = 0; j < nMonomer; ++j) {
695 out << " " << Dbl(temp[j][itr.rank()], 18, 15);
696 }
697 out << std::endl;
698 }
699
700 }
701
702 template <int D>
703 void FieldIo<D>::writeFieldsRGrid(std::string filename,
704 DArray< RField<D> > const & fields,
705 UnitCell<D> const & unitCell)
706 const
707 {
708 std::ofstream file;
709 fileMaster().openOutputFile(filename, file);
710 writeFieldsRGrid(file, fields, unitCell);
711 file.close();
712 }
713
714 template <int D>
715 void FieldIo<D>::readFieldRGrid(std::istream &in,
716 RField<D> & field,
717 UnitCell<D>& unitCell)
718 const
719 {
720 int nMonomer;
721 FieldIo<D>::readFieldHeader(in, nMonomer, unitCell);
722
723 // Only reading in a file with a single field.
724 UTIL_CHECK(nMonomer == 1);
725
726 // If "field" parameter is allocated, check if dimensions match
727 // those of the system's mesh. Otherwise, allocate.
728
729 if (field.isAllocated()) {
730 UTIL_CHECK(field.meshDimensions() == mesh().dimensions());
731 } else {
732 field.allocate(mesh().dimensions());
733 }
734
735 // Read and check input stream mesh dimensions
736 std::string label;
737 in >> label;
738 if (label != "mesh" && label != "ngrid") {
739 std::string msg = "\n";
740 msg += "Error reading field file:\n";
741 msg += "Expected mesh or ngrid, but found [";
742 msg += label;
743 msg += "]";
744 UTIL_THROW(msg.c_str());
745 }
746 IntVec<D> nGrid;
747 in >> nGrid;
748 UTIL_CHECK(nGrid == mesh().dimensions());
749
750 // Setup temporary workspace.
751 RField<D> temp;
752 temp.allocate(mesh().dimensions());
753
754 // Read Field;
755 MeshIterator<D> itr(mesh().dimensions());
756 for (itr.begin(); !itr.atEnd(); ++itr) {
757 in >> std::setprecision(15) >> temp[itr.rank()];
758 }
759
760 int p = 0;
761 int q = 0;
762 int r = 0;
763 int s = 0;
764 int n1 = 0;
765 int n2 = 0;
766 int n3 = 0;
767
768 if (D==3) {
769 while (n1 < mesh().dimension(0)) {
770 q = p;
771 n2 = 0;
772 while (n2 < mesh().dimension(1)) {
773 r = q;
774 n3 = 0;
775 while (n3 < mesh().dimension(2)) {
776 field[s] = temp[r];
777 r = r + (mesh().dimension(0) * mesh().dimension(1));
778 ++s;
779 ++n3;
780 }
781 q = q + mesh().dimension(0);
782 ++n2;
783 }
784 ++n1;
785 ++p;
786 }
787 }
788
789 else if (D==2) {
790 while (n1 < mesh().dimension(0)) {
791 r =q;
792 n2 = 0;
793 while (n2 < mesh().dimension(1)) {
794 field[s] = temp[r];
795 r = r + (mesh().dimension(0));
796 ++s;
797 ++n2;
798 }
799 ++q;
800 ++n1;
801 }
802 }
803
804 else if (D==1) {
805
806 while (n1 < mesh().dimension(0)) {
807 field[s] = temp[r];
808 ++r;
809 ++s;
810 ++n1;
811 }
812 }
813
814 else{
815 Log::file() << "Invalid Dimensions";
816 }
817 }
818
819 template <int D>
820 void FieldIo<D>::readFieldRGrid(std::string filename,
821 RField<D> & field,
822 UnitCell<D>& unitCell)
823 const
824 {
825 std::ifstream file;
826 fileMaster().openInputFile(filename, file);
827 readFieldRGrid(file, field, unitCell);
828 file.close();
829 }
830
831 template <int D>
832 void FieldIo<D>::writeFieldRGrid(std::ostream &out,
833 RField<D> const & field,
834 UnitCell<D> const & unitCell,
835 bool writeHeader)
836 const
837 {
838 if (writeHeader) {
839 writeFieldHeader(out, 1, unitCell);
840 out << "mesh " << std::endl
841 << " " << mesh().dimensions() << std::endl;
842 }
843
844 RField<D> temp;
845 temp.allocate(mesh().dimensions());
846
847 int p = 0;
848 int q = 0;
849 int r = 0;
850 int s = 0;
851 int n1 =0;
852 int n2 =0;
853 int n3 =0;
854
855 if (D==3) {
856 while (n3 < mesh().dimension(2)) {
857 q = p;
858 n2 = 0;
859 while (n2 < mesh().dimension(1)) {
860 r =q;
861 n1 = 0;
862 while (n1 < mesh().dimension(0)) {
863 temp[s] = field[r];
864 r = r + (mesh().dimension(1) * mesh().dimension(2));
865 ++s;
866 ++n1;
867 }
868 q = q + mesh().dimension(2);
869 ++n2;
870 }
871 ++n3;
872 ++p;
873 }
874 }
875 else if (D==2) {
876 while (n2 < mesh().dimension(1)) {
877 r =q;
878 n1 = 0;
879 while (n1 < mesh().dimension(0)) {
880 temp[s] = field[r];
881 r = r + (mesh().dimension(1));
882 ++s;
883 ++n1;
884 }
885 ++q;
886 ++n2;
887 }
888 }
889 else if (D==1) {
890 while (n1 < mesh().dimension(0)) {
891 temp[s] = field[r];
892 ++r;
893 ++s;
894 ++n1;
895 }
896 } else {
897 Log::file() << "Invalid Dimensions";
898 }
899
900 // Write field
901 MeshIterator<D> itr(mesh().dimensions());
902 for (itr.begin(); !itr.atEnd(); ++itr) {
903 out << " " << Dbl(temp[itr.rank()], 18, 15);
904 out << std::endl;
905 }
906 }
907
908 template <int D>
909 void FieldIo<D>::writeFieldRGrid(std::string filename,
910 RField<D> const & field,
911 UnitCell<D> const & unitCell)
912 const
913 {
914 std::ofstream file;
915 fileMaster().openOutputFile(filename, file);
916 writeFieldRGrid(file, field, unitCell);
917 file.close();
918 }
919
920 template <int D>
921 void FieldIo<D>::readFieldsKGrid(std::istream &in,
922 DArray<RFieldDft<D> >& fields,
923 UnitCell<D>& unitCell)
924 const
925 {
926 int nMonomer;
927 FieldIo<D>::readFieldHeader(in, nMonomer, unitCell);
928
929 // If "fields" parameter is allocated, check if dimensions match
930 // those of the system's mesh. Otherwise, allocate.
931
932 if (fields.isAllocated()) {
933
934 int nMonomerFields = fields.capacity();
935
936 UTIL_CHECK(nMonomerFields > 0)
937 UTIL_CHECK(nMonomerFields == nMonomer)
938
939 for (int i = 0; i < nMonomer; ++i) {
940 UTIL_CHECK(fields[i].meshDimensions() == mesh().dimensions());
941 }
942
943 } else {
944
945 fields.allocate(nMonomer);
946 for (int i = 0; i < nMonomer; ++i) {
947 fields[i].allocate(mesh().dimensions());
948 }
949
950 }
951
952 // Read and check input stream mesh dimensions
953 std::string label;
954 in >> label;
955 if (label != "mesh" && label != "ngrid") {
956 std::string msg = "\n";
957 msg += "Error reading field file:\n";
958 msg += "Expected mesh or ngrid, but found [";
959 msg += label;
960 msg += "]";
961 UTIL_THROW(msg.c_str());
962 }
963 IntVec<D> nGrid;
964 in >> nGrid;
965 UTIL_CHECK(nGrid == mesh().dimensions());
966
967 // Read fields;
968 int i, idum;
969 MeshIterator<D> itr(fields[0].dftDimensions());
970 i = 0;
971 for (itr.begin(); !itr.atEnd(); ++itr) {
972 in >> idum;
973 UTIL_CHECK(i == idum);
974 UTIL_CHECK(i == itr.rank());
975 for (int i = 0; i < nMonomer; ++i) {
976 for (int j = 0; j < 2; ++j) {
977 in >> fields[i][itr.rank()][j];
978 }
979 }
980 ++i;
981 }
982 }
983
984 template <int D>
985 void FieldIo<D>::readFieldsKGrid(std::string filename,
986 DArray< RFieldDft<D> >& fields,
987 UnitCell<D>& unitCell)
988 const
989 {
990 std::ifstream file;
991 fileMaster().openInputFile(filename, file);
992 readFieldsKGrid(file, fields, unitCell);
993 file.close();
994 }
995
996 template <int D>
997 void FieldIo<D>::writeFieldsKGrid(std::ostream &out,
998 DArray<RFieldDft<D> > const & fields,
999 UnitCell<D> const & unitCell)
1000 const
1001 {
1002 // Inspect fields array
1003 int nMonomer = fields.capacity();
1004 UTIL_CHECK(nMonomer > 0);
1005 for (int i = 0; i < nMonomer; ++i) {
1006 UTIL_CHECK(fields[i].meshDimensions() == mesh().dimensions());
1007 }
1008
1009 // Write header
1010 writeFieldHeader(out, nMonomer, unitCell);
1011 out << "mesh " << std::endl
1012 << " " << mesh().dimensions() << std::endl;
1013
1014 // Write fields
1015 MeshIterator<D> itr(fields[0].dftDimensions());
1016 int i = 0;
1017 for (itr.begin(); !itr.atEnd(); ++itr) {
1018 UTIL_CHECK(i == itr.rank());
1019 out << Int(itr.rank(), 5);
1020 for (int j = 0; j < nMonomer; ++j) {
1021 out << " "
1022 << Dbl(fields[j][itr.rank()][0], 20, 12)
1023 << Dbl(fields[j][itr.rank()][1], 20, 12);
1024 }
1025 out << std::endl;
1026 ++i;
1027 }
1028 }
1029
1030 template <int D>
1031 void FieldIo<D>::writeFieldsKGrid(std::string filename,
1032 DArray< RFieldDft<D> > const & fields,
1033 UnitCell<D> const & unitCell)
1034 const
1035 {
1036 std::ofstream file;
1037 fileMaster().openOutputFile(filename, file);
1038 writeFieldsKGrid(file, fields, unitCell);
1039 file.close();
1040 }
1041
1042 /*
1043 * Read common part of field header and extract
1044 * the number of monomers (number of fields) in the file.
1045 */
1046 template <int D>
1047 void FieldIo<D>::readFieldHeader(std::istream& in,
1048 int& nMonomer,
1049 UnitCell<D>& unitCell)
1050 const
1051 {
1052 // Preconditions
1053 UTIL_CHECK(latticePtr_);
1054 UTIL_CHECK(groupNamePtr_);
1055 if (unitCell.lattice() == UnitCell<D>::Null) {
1056 UTIL_CHECK(unitCell.nParameter() == 0);
1057 } else {
1058 UTIL_CHECK(unitCell.nParameter() > 0);
1059 if (lattice() != UnitCell<D>::Null) {
1060 UTIL_CHECK(unitCell.lattice() == lattice());
1061 }
1062 }
1063
1064 // Read field header to set unitCell, groupNameIn, nMonomer
1065 int ver1, ver2;
1066 std::string groupNameIn;
1067
1068 Pscf::readFieldHeader(in, ver1, ver2, unitCell,
1069 groupNameIn, nMonomer);
1070 // Note: Function definition in pscf/crystal/UnitCell.tpp
1071
1072 // Checks of data from header
1073 UTIL_CHECK(ver1 == 1);
1074 UTIL_CHECK(ver2 == 0);
1075 UTIL_CHECK(unitCell.isInitialized());
1076 UTIL_CHECK(unitCell.lattice() != UnitCell<D>::Null);
1077 UTIL_CHECK(unitCell.nParameter() > 0);
1078
1079 // Validate or initialize lattice type
1080 if (lattice() == UnitCell<D>::Null) {
1081 lattice() = unitCell.lattice();
1082 } else {
1083 if (lattice() != unitCell.lattice()) {
1084 Log::file() << std::endl
1085 << "Error - "
1086 << "Mismatched lattice types, FieldIo::readFieldHeader:\n"
1087 << " FieldIo::lattice :" << lattice() << "\n"
1088 << " Unit cell lattice :" << unitCell.lattice()
1089 << "\n";
1090 UTIL_THROW("Mismatched lattice types");
1091 }
1092 }
1093
1094 // Validate or initialize group name
1095 if (groupName() == "") {
1096 groupName() = groupNameIn;
1097 } else {
1098 if (groupNameIn != groupName()) {
1099 Log::file() << std::endl
1100 << "Error - "
1101 << "Mismatched group names in FieldIo::readFieldHeader:\n"
1102 << " FieldIo::groupName :" << groupName() << "\n"
1103 << " Field file header :" << groupNameIn << "\n";
1104 UTIL_THROW("Mismatched group names");
1105 }
1106 }
1107
1108 // Check group, read from file if necessary
1109 UTIL_CHECK(groupPtr_);
1110 if (group().size() == 1) {
1111 if (groupName() != "I") {
1112 readGroup(groupName(), group());
1113 }
1114 }
1115
1116 // Check basis, construct if not initialized
1117 UTIL_CHECK(basisPtr_);
1118 if (!basis().isInitialized()) {
1119 basisPtr_->makeBasis(mesh(), unitCell, group());
1120 }
1121 UTIL_CHECK(basis().isInitialized());
1122
1123 }
1124
1125 template <int D>
1126 void FieldIo<D>::writeFieldHeader(std::ostream &out, int nMonomer,
1127 UnitCell<D> const & unitCell) const
1128 {
1129 int ver1 = 1;
1130 int ver2 = 0;
1131 Pscf::writeFieldHeader(out, ver1, ver2, unitCell,
1132 groupName(), nMonomer);
1133 // Note: This function is defined in pscf/crystal/UnitCell.tpp
1134 }
1135
1136 template <int D>
1137 void FieldIo<D>::convertBasisToKGrid(DArray<double> const & in,
1138 RFieldDft<D>& out) const
1139 {
1140 UTIL_CHECK(basis().isInitialized());
1141
1142 // Create Mesh<D> with dimensions of DFT Fourier grid.
1143 Mesh<D> dftMesh(out.dftDimensions());
1144
1145 typename Basis<D>::Star const* starPtr; // pointer to current star
1146 typename Basis<D>::Wave const* wavePtr; // pointer to current wave
1147 std::complex<double> component; // coefficient for star
1148 std::complex<double> coeff; // coefficient for wave
1149 IntVec<D> indices; // dft grid indices of wave
1150 int rank; // dft grid rank of wave
1151 int is; // star index
1152 int ib; // basis index
1153 int iw; // wave index
1154
1155 // Initialize all dft coponents to zero
1156 for (rank = 0; rank < dftMesh.size(); ++rank) {
1157 out[rank][0] = 0.0;
1158 out[rank][1] = 0.0;
1159 }
1160
1161 // Loop over stars, skipping cancelled stars
1162 is = 0;
1163 while (is < basis().nStar()) {
1164 starPtr = &(basis().star(is));
1165
1166 if (starPtr->cancel) {
1167 ++is;
1168 continue;
1169 }
1170
1171 // Set basisId for uncancelled star
1172 ib = starPtr->basisId;
1173
1174 if (starPtr->invertFlag == 0) {
1175
1176 // Make complex coefficient for star basis function
1177 component = std::complex<double>(in[ib], 0.0);
1178
1179 // Loop over waves in closed star
1180 for (iw = starPtr->beginId; iw < starPtr->endId; ++iw) {
1181 wavePtr = &basis().wave(iw);
1182 if (!wavePtr->implicit) {
1183 coeff = component*(wavePtr->coeff);
1184 indices = wavePtr->indicesDft;
1185 rank = dftMesh.rank(indices);
1186 out[rank][0] = coeff.real();
1187 out[rank][1] = coeff.imag();
1188 }
1189 }
1190 ++is;
1191
1192 } else
1193 if (starPtr->invertFlag == 1) {
1194
1195 // Loop over waves in first star
1196 component = std::complex<double>(in[ib], -in[ib+1]);
1197 component /= sqrt(2.0);
1198 starPtr = &(basis().star(is));
1199 for (iw = starPtr->beginId; iw < starPtr->endId; ++iw) {
1200 wavePtr = &basis().wave(iw);
1201 if (!(wavePtr->implicit)) {
1202 coeff = component*(wavePtr->coeff);
1203 indices = wavePtr->indicesDft;
1204 rank = dftMesh.rank(indices);
1205 out[rank][0] = coeff.real();
1206 out[rank][1] = coeff.imag();
1207 }
1208 }
1209
1210 // Loop over waves in second star
1211 starPtr = &(basis().star(is+1));
1212 UTIL_CHECK(starPtr->invertFlag == -1);
1213 component = std::complex<double>(in[ib], +in[ib+1]);
1214 component /= sqrt(2.0);
1215 for (iw = starPtr->beginId; iw < starPtr->endId; ++iw) {
1216 wavePtr = &basis().wave(iw);
1217 if (!(wavePtr->implicit)) {
1218 coeff = component*(wavePtr->coeff);
1219 indices = wavePtr->indicesDft;
1220 rank = dftMesh.rank(indices);
1221 out[rank][0] = coeff.real();
1222 out[rank][1] = coeff.imag();
1223 }
1224 }
1225
1226 // Increment is by 2 (two stars were processed)
1227 is += 2;
1228
1229 } else {
1230
1231 UTIL_THROW("Invalid invertFlag value");
1232
1233 }
1234
1235 }
1236
1237 }
1238
1239 template <int D>
1240 void FieldIo<D>::convertKGridToBasis(RFieldDft<D> const & in,
1241 DArray<double>& out,
1242 bool checkSymmetry,
1243 double epsilon) const
1244 {
1245 UTIL_CHECK(basis().isInitialized());
1246
1247 if (checkSymmetry) {
1248 // Check if kgrid has symmetry
1249 bool symmetric = hasSymmetry(in, epsilon, true);
1250 if (!symmetric) {
1251 Log::file() << std::endl
1252 << "WARNING: non-negligible error in conversion to "
1253 << "symmetry-adapted basis format." << std::endl
1254 << " See error values printed above for each "
1255 << "asymmetric field." << std::endl
1256 << " The field that is output by the above operation "
1257 << "will be a" << std::endl
1258 << " symmetrized version of the input field."
1259 << std::endl << std::endl;
1260 }
1261 }
1262
1263 // Create Mesh<D> with dimensions of DFT Fourier grid.
1264 Mesh<D> dftMesh(in.dftDimensions());
1265
1266 typename Basis<D>::Star const* starPtr; // pointer to current star
1267 typename Basis<D>::Wave const* wavePtr; // pointer to current wave
1268 std::complex<double> component; // coefficient for star
1269 int rank; // dft grid rank of wave
1270 int is; // star index
1271 int ib; // basis index
1272
1273 // Initialize all components to zero
1274 for (is = 0; is < basis().nBasis(); ++is) {
1275 out[is] = 0.0;
1276 }
1277
1278 // Loop over stars
1279 is = 0;
1280 while (is < basis().nStar()) {
1281 starPtr = &(basis().star(is));
1282
1283 if (starPtr->cancel) {
1284 ++is;
1285 continue;
1286 }
1287
1288 // Set basis id for uncancelled star
1289 ib = starPtr->basisId;
1290
1291 if (starPtr->invertFlag == 0) {
1292
1293 // Choose a wave in the star that is not implicit
1294 int beginId = starPtr->beginId;
1295 int endId = starPtr->endId;
1296 int iw = 0;
1297 bool isImplicit = true;
1298 while (isImplicit) {
1299 wavePtr = &basis().wave(beginId + iw);
1300 if (!wavePtr->implicit) {
1301 isImplicit = false;
1302 } else {
1303 UTIL_CHECK(beginId + iw < endId - 1 - iw);
1304 wavePtr = &basis().wave(endId - 1 - iw);
1305 if (!wavePtr->implicit) {
1306 isImplicit = false;
1307 }
1308 }
1309 ++iw;
1310 }
1311 UTIL_CHECK(wavePtr->starId == is);
1312
1313 // Compute component value
1314 rank = dftMesh.rank(wavePtr->indicesDft);
1315 component = std::complex<double>(in[rank][0], in[rank][1]);
1316 component /= wavePtr->coeff;
1317 out[ib] = component.real();
1318 ++is;
1319
1320 } else
1321 if (starPtr->invertFlag == 1) {
1322
1323 // Identify a characteristic wave that is not implicit:
1324 // Either the first wave of the 1st star or last wave of 2nd
1325 wavePtr = &basis().wave(starPtr->beginId);
1326 UTIL_CHECK(wavePtr->starId == is);
1327 if (wavePtr->implicit) {
1328 starPtr = &(basis().star(is+1));
1329 UTIL_CHECK(starPtr->invertFlag == -1);
1330 wavePtr = &basis().wave(starPtr->endId - 1);
1331 UTIL_CHECK(!(wavePtr->implicit));
1332 UTIL_CHECK(wavePtr->starId == is+1);
1333 }
1334 rank = dftMesh.rank(wavePtr->indicesDft);
1335 component = std::complex<double>(in[rank][0], in[rank][1]);
1336 UTIL_CHECK(std::abs(wavePtr->coeff) > 1.0E-8);
1337 component /= wavePtr->coeff;
1338 component *= sqrt(2.0);
1339
1340 // Compute basis function coefficient values
1341 if (starPtr->invertFlag == 1) {
1342 out[ib] = component.real();
1343 out[ib+1] = -component.imag();
1344 } else {
1345 out[ib] = component.real();
1346 out[ib+1] = component.imag();
1347 }
1348
1349 is += 2;
1350 } else {
1351 UTIL_THROW("Invalid invertFlag value");
1352 }
1353
1354 } // loop over star index is
1355 }
1356
1357 template <int D>
1358 void
1359 FieldIo<D>::convertBasisToKGrid(DArray< DArray <double> > const & in,
1360 DArray< RFieldDft<D> >& out) const
1361 {
1362 UTIL_ASSERT(in.capacity() == out.capacity());
1363 int n = in.capacity();
1364 for (int i = 0; i < n; ++i) {
1365 convertBasisToKGrid(in[i], out[i]);
1366 }
1367 }
1368
1369 template <int D>
1370 void FieldIo<D>::convertKGridToBasis(DArray< RFieldDft<D> > const & in,
1371 DArray< DArray <double> > & out,
1372 bool checkSymmetry,
1373 double epsilon) const
1374 {
1375 UTIL_ASSERT(in.capacity() == out.capacity());
1376 int n = in.capacity();
1377
1378 bool symmetric(true);
1379 for (int i = 0; i < n; ++i) {
1380 if (checkSymmetry) {
1381 // Check if kgrid has symmetry
1382 bool tmp_sym = hasSymmetry(in[i], epsilon, true);
1383 if (!tmp_sym) symmetric = false;
1384 }
1385 convertKGridToBasis(in[i], out[i], false);
1386 }
1387
1388 // Print warning if any input field is assymmetric
1389 if (!symmetric) {
1390 Log::file() << std::endl
1391 << "WARNING: non-negligible error in conversion to "
1392 << "symmetry-adapted basis format." << std::endl
1393 << "See error values printed above for each asymmetric field."
1394 << std::endl
1395 << "The field that is output by this operation will be "
1396 << "a symmetrized version of" << std::endl
1397 << "the input field." << std::endl << std::endl;
1398 }
1399 }
1400
1401 template <int D>
1402 void
1403 FieldIo<D>::convertBasisToRGrid(DArray<double> const & in,
1404 RField<D>& out) const
1405 {
1406 checkWorkDft();
1407 convertBasisToKGrid(in, workDft_);
1408 fft().inverseTransformSafe(workDft_, out);
1409 }
1410
1411 template <int D>
1412 void
1413 FieldIo<D>::convertBasisToRGrid(DArray< DArray <double> > const & in,
1414 DArray< RField<D> >& out) const
1415 {
1416 UTIL_ASSERT(in.capacity() == out.capacity());
1417 checkWorkDft();
1418
1419 int n = in.capacity();
1420 for (int i = 0; i < n; ++i) {
1421 convertBasisToKGrid(in[i], workDft_);
1422 fft().inverseTransformSafe(workDft_, out[i]);
1423 }
1424 }
1425
1426 template <int D>
1427 void
1428 FieldIo<D>::convertRGridToBasis(RField<D> const & in,
1429 DArray<double> & out,
1430 bool checkSymmetry,
1431 double epsilon) const
1432 {
1433 checkWorkDft();
1434 fft().forwardTransform(in, workDft_);
1435 convertKGridToBasis(workDft_, out, checkSymmetry, epsilon);
1436 }
1437
1438 template <int D>
1439 void
1440 FieldIo<D>::convertRGridToBasis(DArray< RField<D> > const & in,
1441 DArray< DArray <double> > & out,
1442 bool checkSymmetry,
1443 double epsilon) const
1444 {
1445 UTIL_ASSERT(in.capacity() == out.capacity());
1446 checkWorkDft();
1447
1448 int n = in.capacity();
1449
1450 bool symmetric(true);
1451 for (int i = 0; i < n; ++i) {
1452 fft().forwardTransform(in[i], workDft_);
1453 if (checkSymmetry) {
1454 // Check if kgrid has symmetry
1455 bool tmp_sym = hasSymmetry(workDft_, epsilon, true);
1456 if (!tmp_sym) symmetric = false;
1457 }
1458 convertKGridToBasis(workDft_, out[i], false);
1459 }
1460
1461 // Print warning if any input fields is asymmetric
1462 if (!symmetric) {
1463 Log::file() << std::endl
1464 << "WARNING: non-negligible error in conversion to "
1465 << "symmetry-adapted basis format." << std::endl
1466 << " See error values printed above for each "
1467 << "asymmetric field." << std::endl
1468 << " The field that is output by the above operation "
1469 << "will be a" << std::endl
1470 << " symmetrized version of the input field."
1471 << std::endl << std::endl;
1472 }
1473 }
1474
1475 /*
1476 * Apply inverse FFT to an array of k-grid fields.
1477 */
1478 template <int D>
1479 void
1480 FieldIo<D>::convertKGridToRGrid(DArray< RFieldDft<D> > & in,
1481 DArray< RField<D> >& out) const
1482 {
1483 UTIL_ASSERT(in.capacity() == out.capacity());
1484 int n = in.capacity();
1485 for (int i = 0; i < n; ++i) {
1486 fft().inverseTransformSafe(in[i], out[i]);
1487 }
1488 }
1489
1490 /*
1491 * Apply inverse FFT to a single k-grid field.
1492 */
1493 template <int D>
1494 void
1495 FieldIo<D>::convertKGridToRGrid(RFieldDft<D>& in, RField<D>& out) const
1496 {
1497 fft().inverseTransformSafe(in, out);
1498 }
1499
1500 /*
1501 * Apply forward FFT to an array of r-grid fields.
1502 */
1503 template <int D>
1504 void
1505 FieldIo<D>::convertRGridToKGrid(DArray< RField<D> > const & in,
1506 DArray< RFieldDft<D> >& out) const
1507 {
1508 UTIL_ASSERT(in.capacity() == out.capacity());
1509 int n = in.capacity();
1510 for (int i = 0; i < n; ++i) {
1511 fft().forwardTransform(in[i], out[i]);
1512 }
1513 }
1514
1515 /*
1516 * Apply forward FFT to a single r-grid field.
1517 */
1518 template <int D>
1519 void
1520 FieldIo<D>::convertRGridToKGrid(RField<D> const & in,
1521 RFieldDft<D>& out) const
1522 {
1523 fft().forwardTransform(in, out);
1524 }
1525
1526 /*
1527 * Test if an RField<D> has declared space group symmetry.
1528 * Return true if symmetric, false otherwise. Print error values
1529 * if verbose == true and hasSymmetry == false.
1530 */
1531 template <int D>
1532 bool FieldIo<D>::hasSymmetry(RField<D> const & in, double epsilon,
1533 bool verbose) const
1534 {
1535 checkWorkDft();
1536 fft().forwardTransform(in, workDft_);
1537 return hasSymmetry(workDft_, epsilon, verbose);
1538 }
1539
1540 /*
1541 * Test if an RFieldDft has the declared space group symmetry.
1542 * Return true if symmetric, false otherwise. Print error values
1543 * if verbose == true and hasSymmetry == false.
1544 */
1545 template <int D>
1546 bool FieldIo<D>::hasSymmetry(RFieldDft<D> const & in, double epsilon,
1547 bool verbose) const
1548 {
1549 UTIL_CHECK(basis().isInitialized());
1550
1551 typename Basis<D>::Star const* starPtr; // pointer to current star
1552 typename Basis<D>::Wave const* wavePtr; // pointer to current wave
1553 std::complex<double> waveCoeff; // coefficient from wave
1554 std::complex<double> rootCoeff; // coefficient from root
1555 std::complex<double> diff; // coefficient difference
1556 int is; // star index
1557 int iw; // wave index
1558 int beginId, endId; // star begin, end ids
1559 int rank; // dft grid rank of wave
1560
1561 double cancelledError(0.0); // max error from cancelled stars
1562 double uncancelledError(0.0); // max error from uncancelled stars
1563
1564 // Create Mesh<D> with dimensions of DFT Fourier grid.
1565 Mesh<D> dftMesh(in.dftDimensions());
1566
1567 // Loop over all stars
1568 for (is = 0; is < basis().nStar(); ++is) {
1569 starPtr = &(basis().star(is));
1570
1571 if (starPtr->cancel) {
1572
1573 // Check that coefficients are zero for all waves in star
1574 beginId = starPtr->beginId;
1575 endId = starPtr->endId;
1576 for (iw = beginId; iw < endId; ++iw) {
1577 wavePtr = &basis().wave(iw);
1578 if (!wavePtr->implicit) {
1579 rank = dftMesh.rank(wavePtr->indicesDft);
1580 waveCoeff = std::complex<double>(in[rank][0], in[rank][1]);
1581 if (std::abs(waveCoeff) > cancelledError) {
1582 cancelledError = std::abs(waveCoeff);
1583 if ((!verbose) && (cancelledError > epsilon)) {
1584 return false;
1585 }
1586 }
1587 }
1588 }
1589
1590 } else {
1591
1592 // Check consistency of coeff values from all waves
1593 bool hasRoot = false;
1594 beginId = starPtr->beginId;
1595 endId = starPtr->endId;
1596 for (iw = beginId; iw < endId; ++iw) {
1597 wavePtr = &basis().wave(iw);
1598 if (!(wavePtr->implicit)) {
1599 rank = dftMesh.rank(wavePtr->indicesDft);
1600 waveCoeff = std::complex<double>(in[rank][0], in[rank][1]);
1601 waveCoeff /= wavePtr->coeff;
1602 if (hasRoot) {
1603 diff = waveCoeff - rootCoeff;
1604 if (std::abs(diff) > uncancelledError) {
1605 uncancelledError = std::abs(diff);
1606 if ((!verbose) && (uncancelledError > epsilon)) {
1607 return false;
1608 }
1609 }
1610 } else {
1611 rootCoeff = waveCoeff;
1612 hasRoot = true;
1613 }
1614 }
1615 }
1616
1617 }
1618
1619 } // end loop over star index is
1620
1621 if ((cancelledError < epsilon) && (uncancelledError < epsilon)) {
1622 return true;
1623 } else if (verbose) {
1624 Log::file() << std::endl
1625 << "Maximum coefficient of a cancelled star: "
1626 << cancelledError << std::endl
1627 << "Maximum error of coefficient for uncancelled star: "
1628 << uncancelledError << std::endl;
1629 }
1630 return false;
1631 }
1632
1633 template <int D>
1634 void FieldIo<D>::checkWorkDft() const
1635 {
1636 if (!workDft_.isAllocated()) {
1637 workDft_.allocate(mesh().dimensions());
1638 } else {
1639 UTIL_CHECK(workDft_.meshDimensions() == fft().meshDimensions());
1640 }
1641 }
1642
1643} // namespace Pspc
1644} // namespace Pscf
1645#endif
Pscf::Basis::Star
A list of wavevectors that are related by space-group symmetries.
Definition: Basis.h:439
Pscf::Basis::Star::size
int size
Number of wavevectors in this star.
Definition: Basis.h:446
Pscf::Basis::Star::invertFlag
int invertFlag
Index for inversion symmetry of star.
Definition: Basis.h:485
Pscf::Basis::Star::waveBz
IntVec< D > waveBz
Integer indices indexBz of a characteristic wave of this star.
Definition: Basis.h:498
Pscf::Basis::Star::endId
int endId
Wave index of last wavevector in star.
Definition: Basis.h:456
Pscf::Basis::Star::basisId
int basisId
Index of basis function associated with this star.
Definition: Basis.h:517
Pscf::Basis::Star::beginId
int beginId
Wave index of first wavevector in star.
Definition: Basis.h:451
Pscf::Basis::Star::cancel
bool cancel
Is this star cancelled, i.e., associated with a zero function?
Definition: Basis.h:532
Pscf::Basis::Wave
Wavevector used to construct a basis function.
Definition: Basis.h:354
Pscf::Basis::Wave::implicit
bool implicit
Is this wave represented implicitly in DFT of real field?
Definition: Basis.h:405
Pscf::Basis::Wave::coeff
std::complex< double > coeff
Coefficient of wave within the associated star basis function.
Definition: Basis.h:361
Pscf::Basis::Wave::starId
int starId
Index of the star that contains this wavevector.
Definition: Basis.h:390
Pscf::Basis::Wave::indicesDft
IntVec< D > indicesDft
Integer indices of wave, on a discrete Fourier transform mesh.
Definition: Basis.h:370
Pscf::Basis
Symmetry-adapted Fourier basis for pseudo-spectral scft.
Definition: Basis.h:346
Pscf::IntVec
An IntVec<D, T> is a D-component vector of elements of integer type T.
Definition: IntVec.h:27
Pscf::MeshIterator
Iterator over points in a Mesh<D>.
Definition: MeshIterator.h:29
Pscf::MeshIterator::rank
int rank() const
Get the rank of current element.
Definition: MeshIterator.h:136
Pscf::MeshIterator::begin
void begin()
Set iterator to the first point in the mesh.
Definition: MeshIterator.tpp:61
Pscf::MeshIterator::atEnd
bool atEnd() const
Is this the end (i.e., one past the last point)?
Definition: MeshIterator.h:141
Pscf::Mesh
Description of a regular grid of points in a periodic domain.
Definition: Mesh.h:61
Pscf::Mesh::size
int size() const
Get total number of grid points.
Definition: Mesh.h:229
Pscf::Mesh::rank
int rank(IntVec< D > const &position) const
Get the rank of a grid point with specified position.
Definition: Mesh.tpp:72
Pscf::Pspc::FFT
Fourier transform wrapper for real data.
Definition: pspc/field/FFT.h:31
Pscf::Pspc::FieldIo
File input/output operations and format conversions for fields.
Definition: pspc/field/FieldIo.h:43
Pscf::Pspc::FieldIo::readFieldRGrid
void readFieldRGrid(std::istream &in, RField< D > &field, UnitCell< D > &unitCell) const
Read single RField (field on an r-space grid) from istream.
Definition: pspc/field/FieldIo.tpp:715
Pscf::Pspc::FieldIo::writeFieldsBasis
void writeFieldsBasis(std::ostream &out, DArray< DArray< double > > const &fields, UnitCell< D > const &unitCell) const
Write concentration or chemical potential field components to file.
Definition: pspc/field/FieldIo.tpp:430
Pscf::Pspc::FieldIo::convertBasisToKGrid
void convertBasisToKGrid(DArray< double > const &components, RFieldDft< D > &dft) const
Convert a field from symmetrized basis to Fourier transform (k-grid).
Definition: pspc/field/FieldIo.tpp:1137
Pscf::Pspc::FieldIo::readFieldsKGrid
void readFieldsKGrid(std::istream &in, DArray< RFieldDft< D > > &fields, UnitCell< D > &unitCell) const
Read array of RFieldDft objects (k-space fields) from file.
Definition: pspc/field/FieldIo.tpp:921
Pscf::Pspc::FieldIo::readFieldHeader
void readFieldHeader(std::istream &in, int &nMonomer, UnitCell< D > &unitCell) const
Reader header of field file (fortran pscf format)
Definition: pspc/field/FieldIo.tpp:1047
Pscf::Pspc::FieldIo::writeFieldHeader
void writeFieldHeader(std::ostream &out, int nMonomer, UnitCell< D > const &unitCell) const
Write header for field file (fortran pscf format)
Definition: pspc/field/FieldIo.tpp:1126
Pscf::Pspc::FieldIo::convertBasisToRGrid
void convertBasisToRGrid(DArray< double > const &in, RField< D > &out) const
Convert a field from symmetrized basis to spatial grid (r-grid).
Definition: pspc/field/FieldIo.tpp:1403
Pscf::Pspc::FieldIo::convertKGridToBasis
void convertKGridToBasis(RFieldDft< D > const &in, DArray< double > &out, bool checkSymmetry=true, double epsilon=1.0e-8) const
Convert a field from Fourier transform (kgrid) to symmetrized basis.
Definition: pspc/field/FieldIo.tpp:1240
Pscf::Pspc::FieldIo::writeFieldsKGrid
void writeFieldsKGrid(std::ostream &out, DArray< RFieldDft< D > > const &fields, UnitCell< D > const &unitCell) const
Write array of RFieldDft objects (k-space fields) to file.
Definition: pspc/field/FieldIo.tpp:997
Pscf::Pspc::FieldIo::writeFieldBasis
void writeFieldBasis(std::ostream &out, DArray< double > const &field, UnitCell< D > const &unitCell) const
Write single concentration or chemical potential field to output stream out.
Definition: pspc/field/FieldIo.tpp:106
Pscf::Pspc::FieldIo::associate
void associate(Mesh< D > const &mesh, FFT< D > const &fft, typename UnitCell< D >::LatticeSystem &lattice, std::string &groupName, SpaceGroup< D > &group, Basis< D > &basis, FileMaster const &fileMaster)
Get and store addresses of associated objects.
Definition: pspc/field/FieldIo.tpp:56
Pscf::Pspc::FieldIo::readFieldsRGrid
void readFieldsRGrid(std::istream &in, DArray< RField< D > > &fields, UnitCell< D > &unitCell) const
Read array of RField objects (fields on r-space grid) from istream.
Definition: pspc/field/FieldIo.tpp:478
Pscf::Pspc::FieldIo::convertRGridToBasis
void convertRGridToBasis(RField< D > const &in, DArray< double > &out, bool checkSymmetry=true, double epsilon=1.0e-8) const
Convert a field from spatial grid (r-grid) to symmetrized basis.
Definition: pspc/field/FieldIo.tpp:1428
Pscf::Pspc::FieldIo::writeFieldRGrid
void writeFieldRGrid(std::ostream &out, RField< D > const &field, UnitCell< D > const &unitCell, bool writeHeader=true) const
Write a single RField (field on an r-space grid) to ostream.
Definition: pspc/field/FieldIo.tpp:832
Pscf::Pspc::FieldIo::writeFieldsRGrid
void writeFieldsRGrid(std::ostream &out, DArray< RField< D > > const &fields, UnitCell< D > const &unitCell) const
Write array of RField objects (fields on r-space grid) to ostream.
Definition: pspc/field/FieldIo.tpp:614
Pscf::Pspc::FieldIo::readFieldBasis
void readFieldBasis(std::istream &in, DArray< double > &field, UnitCell< D > &unitCell) const
Read single concentration or chemical potential field from file.
Definition: pspc/field/FieldIo.tpp:74
Pscf::Pspc::FieldIo::convertKGridToRGrid
void convertKGridToRGrid(DArray< RFieldDft< D > > &in, DArray< RField< D > > &out) const
Convert fields from k-grid (DFT) to real space (r-grid) format.
Definition: pspc/field/FieldIo.tpp:1480
Pscf::Pspc::FieldIo::hasSymmetry
bool hasSymmetry(RField< D > const &in, double epsilon=1.0e-8, bool verbose=true) const
Check if an r-grid field has the declared space group symmetry.
Definition: pspc/field/FieldIo.tpp:1532
Pscf::Pspc::FieldIo::convertRGridToKGrid
void convertRGridToKGrid(DArray< RField< D > > const &in, DArray< RFieldDft< D > > &out) const
Convert fields from spatial grid (r-grid) to k-grid format.
Definition: pspc/field/FieldIo.tpp:1505
Pscf::Pspc::FieldIo::readFieldsBasis
void readFieldsBasis(std::istream &in, DArray< DArray< double > > &fields, UnitCell< D > &unitCell) const
Read concentration or chemical potential field components from file.
Definition: pspc/field/FieldIo.tpp:132
Pscf::Pspc::RFieldDft
Fourier transform of a real field on an FFT mesh.
Definition: RFieldDft.h:31
Pscf::Pspc::RFieldDft::dftDimensions
IntVec< D > const & dftDimensions() const
Return vector of dft (Fourier) grid dimensions by constant reference.
Definition: RFieldDft.h:144
Pscf::Pspc::RField
Field of real double precision values on an FFT mesh.
Definition: RField.h:29
Pscf::Pspc::RField::meshDimensions
const IntVec< D > & meshDimensions() const
Return mesh dimensions by constant reference.
Definition: RField.h:102
Pscf::Pspc::RField::allocate
void allocate(const IntVec< D > &meshDimensions)
Allocate the underlying C array for an FFT grid.
Definition: RField.tpp:97
Pscf::SpaceGroup
Crystallographic space group.
Definition: SpaceGroup.h:30
Pscf::UnitCellBase::nParameter
int nParameter() const
Get the number of parameters in the unit cell.
Definition: UnitCellBase.h:243
Pscf::UnitCellBase::isInitialized
bool isInitialized() const
Has this unit cell been initialized?
Definition: UnitCellBase.h:141
Pscf::UnitCell
Base template for UnitCell<D> classes, D=1, 2 or 3.
Definition: UnitCell.h:44
Pscf::Vec::negate
Vec< D, T > & negate(const Vec< D, T > &v)
Return negative of vector v.
Definition: Vec.h:520
Util::Array::capacity
int capacity() const
Return allocated size.
Definition: Array.h:159
Util::DArray
Dynamically allocatable contiguous array template.
Definition: DArray.h:32
Util::DArray::allocate
void allocate(int capacity)
Allocate the underlying C array.
Definition: DArray.h:199
Util::DArray::isAllocated
bool isAllocated() const
Return true if this DArray has been allocated, false otherwise.
Definition: DArray.h:247
Util::Dbl
Wrapper for a double precision number, for formatted ostream output.
Definition: Dbl.h:40
Util::FileMaster
A FileMaster manages input and output files for a simulation.
Definition: FileMaster.h:143
Util::Int
Wrapper for an int, for formatted ostream output.
Definition: Int.h:37
Util::Log::file
static std::ostream & file()
Get log ostream by reference.
Definition: Log.cpp:57
UTIL_CHECK
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition: global.h:68
UTIL_THROW
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
Definition: global.h:51
UTIL_ASSERT
#define UTIL_ASSERT(condition)
Assertion macro suitable for debugging serial or parallel code.
Definition: global.h:75
Pscf
C++ namespace for polymer self-consistent field theory (PSCF).
Definition: BlockDescriptor.cpp:11
Util
Utility classes for scientific computation.
Definition: accumulators.mod:1
Generated on Wed Sep 20 2023 00:52:46 for PSCF by doxygen 1.9.5