PSCF v1.3
fieldIoUtil.tpp
1#ifndef PRDC_FIELD_IO_UTIL_TPP
2#define PRDC_FIELD_IO_UTIL_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 <prdc/crystal/Basis.h>
12#include <prdc/crystal/UnitCell.h>
13#include <prdc/crystal/unitCellHeader.h>
14#include <prdc/crystal/replicateUnitCell.h>
15#include <prdc/crystal/shiftToMinimum.h>
16
17#include <pscf/mesh/Mesh.h>
18#include <pscf/mesh/MeshIterator.h>
19#include <pscf/mesh/MeshIteratorFortran.h>
20#include <pscf/math/IntVec.h>
21#include <pscf/math/complex.h>
22
23#include <util/containers/DArray.h>
24#include <util/misc/FileMaster.h>
25#include <util/misc/Log.h>
26#include <util/format/Int.h>
27#include <util/format/Dbl.h>
28#include <util/math/Constants.h>
29
30#include <string>
31
32namespace Pscf {
33namespace Prdc {
34
35 // Field inspection and allocation functions
36
37 /*
38 * Check allocation of a single field of type FT, allocate if needed.
39 */
40 template <int D, class FT>
41 void checkAllocateField(FT& field,
42 IntVec<D> const& dimensions)
43 {
44 if (field.isAllocated()) {
45 UTIL_CHECK(field.meshDimensions() == dimensions);
46 } else {
47 field.allocate(dimensions);
48 }
49 }
50
51 /*
52 * Check allocation of an array of fields of type FT, allocate if needed.
53 */
54 template <int D, class FT>
56 int nMonomer,
57 IntVec<D> const& dimensions)
58 {
59 if (fields.isAllocated()) {
60 int nMonomerFields = fields.capacity();
61 UTIL_CHECK(nMonomerFields > 0)
62 UTIL_CHECK(nMonomerFields == nMonomer)
63 for (int i = 0; i < nMonomer; ++i) {
64 UTIL_CHECK(fields[i].isAllocated());
65 UTIL_CHECK(fields[i].meshDimensions() == dimensions);
66 }
67 } else {
68 fields.allocate(nMonomer);
69 for (int i = 0; i < nMonomer; ++i) {
70 fields[i].allocate(dimensions);
71 }
72 }
73 }
74
75 /*
76 * Inspect dimensions of a DArray of fields, each of type FT.
77 */
78 template <int D, class FT>
79 void inspectFields(DArray<FT> const& fields,
80 int & nMonomer,
81 IntVec<D> & dimensions)
82 {
83 UTIL_CHECK(fields.isAllocated());
84 nMonomer = fields.capacity();
85 UTIL_CHECK(nMonomer > 0);
86
87 dimensions = fields[0].meshDimensions();
88 for (int i = 0; i < nMonomer; ++i) {
89 UTIL_CHECK(fields[i].isAllocated());
90 UTIL_CHECK(fields[i].meshDimensions() == dimensions);
91 }
92 }
93
94 /*
95 * Check allocation of an array of arrays of type AT, allocate if needed.
96 */
97 template <class AT>
99 int nMonomer,
100 int capacity)
101 {
102 if (arrays.isAllocated()) {
103 int nMonomerArrays = arrays.capacity();
104 UTIL_CHECK(nMonomerArrays > 0)
105 UTIL_CHECK(nMonomerArrays == nMonomer)
106 for (int i = 0; i < nMonomer; ++i) {
107 UTIL_CHECK(arrays[i].isAllocated());
108 UTIL_CHECK(arrays[i].capacity() == capacity);
109 }
110 } else {
111 arrays.allocate(nMonomer);
112 for (int i = 0; i < nMonomer; ++i) {
113 arrays[i].allocate(capacity);
114 }
115 }
116 }
117
118 /*
119 * Inspect dimensions of an array of arrays of type AT.
120 */
121 template <class AT>
122 void inspectArrays(DArray<AT> const& arrays,
123 int & nMonomer,
124 int & capacity)
125 {
126 UTIL_CHECK(arrays.isAllocated());
127 nMonomer = arrays.capacity();
128 UTIL_CHECK(nMonomer > 0);
129
130 capacity = arrays[0].capacity();
131 UTIL_CHECK(capacity > 0);
132 for (int i = 0; i < nMonomer; ++i) {
133 UTIL_CHECK(arrays[i].isAllocated());
134 UTIL_CHECK(arrays[i].capacity() == capacity);
135 }
136 }
137
138 // Field File Header (common to all formats)
139
140 /*
141 * Read common part of field header (fortran PSCF format).
142 */
143 template <int D>
144 void readFieldHeader(std::istream& in,
145 int& ver1, int& ver2,
146 UnitCell<D>& cell,
147 std::string& groupName,
148 int& nMonomer)
149 {
150 std::string label;
151
152 // Read format v1 v2
153 in >> label;
154 UTIL_CHECK(label == "format");
155 in >> ver1;
156 in >> ver2;
157
158 // Read dimension of space
159 in >> label;
160 UTIL_CHECK(label == "dim");
161 int dim;
162 in >> dim;
163 UTIL_CHECK(dim == D);
164
165 // Read unit cell lattice type, number of parameters, parameters
166 // Function declared in UnitCell.h and defined in UnitCell.tpp
167 readUnitCellHeader(in, cell);
168
169 // Optionally read groupName
170 in >> label;
171 groupName = "";
172 if (label == "group_name") {
173 in >> groupName;
174 in >> label;
175 }
176
177 // Read nMonomer value
178 UTIL_CHECK(label == "N_monomer");
179 in >> nMonomer;
180 UTIL_CHECK(nMonomer > 0);
181 }
182
183 /*
184 * Write common part of field header (fortran PSCF format).
185 */
186 template <int D>
187 void writeFieldHeader(std::ostream &out,
188 int ver1, int ver2,
189 UnitCell<D> const & cell,
190 std::string const & groupName,
191 int nMonomer)
192 {
193 // Write file format id and dimension of space D
194 out << "format " << Int(ver1,3) << " " << Int(ver2,3) << std::endl;
195 out << "dim" << std::endl
196 << " " << D << std::endl;
197
198 // Write unit cell lattice type, number of parameters, parameters
199 // Function declared in UnitCell.h and defined in UnitCell.tpp
200 writeUnitCellHeader(out, cell);
201
202 // Write group_name if not empty
203 if (groupName != "") {
204 out << "group_name" << std::endl
205 << " " << groupName << std::endl;
206 }
207
208 // Write nMonomer, number of monomer types
209 if (nMonomer != 0) {
210 UTIL_CHECK(nMonomer > 0);
211 out << "N_monomer" << std::endl
212 << " " << nMonomer << std::endl;
213 }
214
215 // Note: The option to not write nMonomer when the value is zero
216 // is designed to allow this function to be used in contexts in
217 // which the value of N_monomer is irrelevant and possibly unknown,
218 // such as files that contain information about stars in a basis
219 // or group operations.
220
221 }
222
223
224 // RGrid File Io Templates
225
226 template <int D>
227 void readMeshDimensions(std::istream& in,
228 IntVec<D> const& meshDimensions)
229 {
230 // Read and check input stream mesh dimensions
231 std::string label;
232 in >> label;
233 UTIL_ASSERT(in.good());
234 if (label != "mesh" && label != "ngrid") {
235 std::string msg = "\n";
236 msg += "Error reading field file:\n";
237 msg += "Expected mesh or ngrid, but found [";
238 msg += label;
239 msg += "]";
240 UTIL_THROW(msg.c_str());
241 }
242 IntVec<D> meshDimensionsIn;
243 in >> meshDimensionsIn;
244 UTIL_CHECK(in.good());
245 if (meshDimensionsIn != meshDimensions) {
246 Log::file()
247 << "Inconsistent mesh dimensions:\n"
248 << "meshDimensions (expected) = " << meshDimensions << "\n"
249 << "meshDimensions (from file) = " << meshDimensionsIn << "\n";
250 UTIL_THROW("Unexpected mesh dimensions in field file header");
251 }
252 }
253
254 template <int D>
255 void writeMeshDimensions(std::ostream &out,
256 IntVec<D> const& meshDimensions)
257 {
258 out << "mesh " << std::endl
259 << " " << meshDimensions << std::endl;
260 }
261
262 template <int D, class ART>
263 void readRGridData(std::istream& in,
264 DArray<ART>& fields,
265 int nMonomer,
266 IntVec<D> const& dimensions)
267 {
268 UTIL_CHECK(nMonomer > 0);
269 UTIL_CHECK(fields.capacity() == nMonomer);
270
271 MeshIteratorFortran<D> iter(dimensions);
272 int rank;
273 for (iter.begin(); !iter.atEnd(); ++iter) {
274 rank = iter.rank();
275 for (int k = 0; k < nMonomer; ++k) {
276 in >> fields[k][rank];
277 }
278 UTIL_ASSERT(in.good());
279 }
280 UTIL_CHECK(in.good());
281 }
282
283 template <int D, class ART>
284 void readRGridData(std::istream& in,
285 ART& field,
286 IntVec<D> const& dimensions)
287 {
288 MeshIteratorFortran<D> iter(dimensions);
289 int rank;
290 for (iter.begin(); !iter.atEnd(); ++iter) {
291 rank = iter.rank();
292 in >> field[rank];
293 UTIL_ASSERT(in.good());
294 }
295 UTIL_CHECK(in.good());
296 }
297
298 template <int D, class ART>
299 void writeRGridData(std::ostream& out,
300 DArray<ART> const& fields,
301 int nMonomer,
302 IntVec<D> const& dimensions)
303 {
304 UTIL_CHECK(nMonomer > 0);
305 UTIL_CHECK(nMonomer == fields.capacity());
306
307 MeshIteratorFortran<D> iter(dimensions);
308 int rank, j;
309 for (iter.begin(); !iter.atEnd(); ++iter) {
310 rank = iter.rank();
311 for (j = 0; j < nMonomer; ++j) {
312 out << " " << Dbl(fields[j][rank], 21, 13);
313 }
314 out << std::endl;
315 }
316
317 }
318
319 template <int D, class ART>
320 void writeRGridData(std::ostream& out,
321 ART const& field,
322 IntVec<D> const& dimensions)
323 {
324 MeshIteratorFortran<D> iter(dimensions);
325 int rank;
326 for (iter.begin(); !iter.atEnd(); ++iter) {
327 rank = iter.rank();
328 out << " " << Dbl(field[rank], 21, 13);
329 out << std::endl;
330 }
331 }
332
333 // KGrid file IO templates
334
335 template <int D, class ACT>
336 void readKGridData(std::istream& in,
337 DArray<ACT> & fields,
338 int nMonomer,
339 IntVec<D> const& dftDimensions)
340 {
341 typedef typename ACT::Complex CT;
342 typedef typename ACT::Real RT;
343
344 RT x, y;
345 MeshIterator<D> iter(dftDimensions);
346 int rank, i, j, idum;
347 i = 0;
348 for (iter.begin(); !iter.atEnd(); ++iter) {
349 rank = iter.rank();
350 in >> idum;
351 UTIL_CHECK(i == idum);
352 UTIL_CHECK(i == rank);
353 for (j = 0; j < nMonomer; ++j) {
354 in >> x;
355 in >> y;
356 assign<CT, RT>(fields[j][rank], x, y);
357 }
358 UTIL_ASSERT(in.good());
359 ++i;
360 }
361 UTIL_CHECK(in.good());
362 }
363
364 template <int D, class ACT>
365 void readKGridData(std::istream& in,
366 ACT& field,
367 IntVec<D> const& dftDimensions)
368 {
369 typedef typename ACT::Complex CT;
370 typedef typename ACT::Real RT;
371
372 RT x, y;
373 MeshIterator<D> iter(dftDimensions);
374 int rank, idum;
375 int i = 0;
376 for (iter.begin(); !iter.atEnd(); ++iter) {
377 rank = iter.rank();
378 in >> idum;
379 UTIL_CHECK(i == idum);
380 UTIL_CHECK(i == rank);
381 in >> x;
382 in >> y;
383 UTIL_ASSERT(in.good());
384 assign<CT, RT>(field[rank], x, y);
385 ++i;
386 }
387 UTIL_CHECK(in.good());
388 }
389
390 template <int D, class ACT>
391 void writeKGridData(std::ostream& out,
392 DArray<ACT> const& fields,
393 int nMonomer,
394 IntVec<D> const& dftDimensions)
395 {
396 UTIL_CHECK(nMonomer > 0);
397 UTIL_CHECK(nMonomer == fields.capacity());
398
399 typedef typename ACT::Complex CT;
400 typedef typename ACT::Real RT;
401
402 RT x, y;
403 MeshIterator<D> iter(dftDimensions);
404 int rank;
405 int i = 0;
406 for (iter.begin(); !iter.atEnd(); ++iter) {
407 rank = iter.rank();
408 UTIL_CHECK(i == rank);
409 out << Int(rank, 5);
410 for (int j = 0; j < nMonomer; ++j) {
411 x = real<CT, RT>(fields[j][rank]);
412 y = imag<CT, RT>(fields[j][rank]);
413 out << " "
414 << Dbl(x, 21, 13)
415 << Dbl(y, 21, 13);
416 }
417 out << std::endl;
418 ++i;
419 }
420
421 }
422
423 template <int D, class ACT>
424 void writeKGridData(std::ostream& out,
425 ACT const& field,
426 IntVec<D> const& dftDimensions)
427 {
428 typedef typename ACT::Complex CT;
429 typedef typename ACT::Real RT;
430
431 RT x, y;
432 MeshIterator<D> iter(dftDimensions);
433 int rank, i;
434 i = 0;
435 for (iter.begin(); !iter.atEnd(); ++iter) {
436 rank = iter.rank();
437 UTIL_CHECK(i == rank);
438 x = real<CT, RT>(field[rank]);
439 y = imag<CT, RT>(field[rank]);
440 out << Int(rank, 5);
441 out << " "
442 << Dbl(x, 21, 13)
443 << Dbl(y, 21, 13);
444 out << std::endl;
445 ++i;
446 }
447 }
448
449 // Functions for files in symmetry-adapted basis format
450
451 /*
452 * Read the data section for an array of fields in basis format.
453 */
454 template <int D>
455 void readBasisData(std::istream& in,
456 DArray< DArray<double> >& fields,
457 UnitCell<D> const& unitCell,
458 Mesh<D> const& mesh,
459 Basis<D> const& basis,
460 int nBasisIn)
461 {
462 UTIL_CHECK(basis.isInitialized());
463 UTIL_CHECK(fields.isAllocated());
464
465 int nMonomer = fields.capacity();
466 UTIL_CHECK(nMonomer > 0);
467
468 // Initialize all field array elements to zero
469 int fieldCapacity = fields[0].capacity();
470 UTIL_CHECK(fieldCapacity > 0);
471 for (int i = 0; i < nMonomer; ++i) {
472 UTIL_CHECK(fields[i].capacity() == fieldCapacity);
473 for (int j = 0; j < fieldCapacity; ++j) {
474 fields[i][j] = 0.0;
475 }
476 }
477
478 // Reset nBasis = min(nBasisIn, fieldCapacity)
479 int nBasis = nBasisIn;
480 if (fieldCapacity < nBasisIn) {
481 nBasis = fieldCapacity;
482 }
483
484 // Allocate temp arrays used to read in components
485 DArray<double> temp, temp2;
486 temp.allocate(nMonomer);
487 temp2.allocate(nMonomer);
488
489 typename Basis<D>::Star const * starPtr;
490 typename Basis<D>::Star const * starPtr2;
491 IntVec<D> waveIn, waveIn2;
492 int sizeIn, sizeIn2;
493 int starId, starId2;
494 int basisId, basisId2;
495 int waveId, waveId2;
496
497 std::complex<double> coeff, phasor;
498 IntVec<D> waveBz, waveDft;
499 int nReversedPair = 0;
500 bool waveExists, sizeMatches;
501
502 // Loop over stars in input file to read field components
503 int i = 0;
504 while (i < nBasis) {
505
506 // Read next line of data
507 for (int j = 0; j < nMonomer; ++j) {
508 in >> temp[j]; // field components
509 }
510 in >> waveIn; // wave of star
511 in >> sizeIn; // # of waves in star
512 UTIL_CHECK(in.good());
513 ++i;
514
515 sizeMatches = false;
516 waveExists = false;
517
518 // Check if wavevector waveIn is a minimum image
519 waveBz = shiftToMinimum(waveIn, mesh.dimensions(), unitCell);
520 waveExists = (waveIn == waveBz);
521
522 if (waveExists) {
523
524 // Find the star containing waveIn
525 waveDft = waveIn;
526 mesh.shift(waveDft);
527 waveId = basis.waveId(waveDft);
528 starId = basis.wave(waveId).starId;
529 starPtr = &basis.star(starId);
530 UTIL_CHECK(!(starPtr->cancel));
531 basisId = starPtr->basisId;
532
533 if (starPtr->size == sizeIn) {
534 sizeMatches = true;
535 } else {
536 Log::file()
537 << "Warning: Inconsistent star size (line ignored)\n"
538 << "wave from file = " << waveIn << "\n"
539 << "size from file = " << sizeIn << "\n"
540 << "size of star = " << starPtr->size
541 << "\n";
542 sizeMatches = false;
543 }
544
545 }
546
547 if (waveExists && sizeMatches) { // Attempt to process wave
548
549 if (starPtr->invertFlag == 0) {
550
551 if (starPtr->waveBz == waveIn) {
552
553 // Copy components of closed star to fields array
554 for (int j = 0; j < nMonomer; ++j) {
555 fields[j][basisId] = temp[j];
556 }
557
558 } else {
559
560 Log::file()
561 << "Inconsistent wave of closed star on input\n"
562 << "wave from file = " << waveIn << "\n"
563 << "starId of wave = " << starId << "\n"
564 << "waveBz of star = " << starPtr->waveBz
565 << "\n";
566
567 }
568
569 } else {
570
571 // Read the next line
572 for (int j = 0; j < nMonomer; ++j) {
573 in >> temp2[j]; // components of field
574 }
575 in >> waveIn2; // wave of star
576 in >> sizeIn2; // # of wavevectors in star
577 UTIL_CHECK(in.good());
578 ++i;
579
580 // Check that waveIn2 is also in the 1st BZ
581 waveBz =
582 shiftToMinimum(waveIn2, mesh.dimensions(), unitCell);
583 UTIL_CHECK(waveIn2 == waveBz);
584
585 // Identify the star containing waveIn2
586 waveDft = waveIn2;
587 mesh.shift(waveDft);
588 waveId2 = basis.waveId(waveDft);
589 starId2 = basis.wave(waveId2).starId;
590 starPtr2 = &basis.star(starId2);
591 UTIL_CHECK(!(starPtr2->cancel));
592 basisId2 = starPtr2->basisId;
593 UTIL_CHECK(starPtr2->size == sizeIn2);
594 UTIL_CHECK(sizeIn == sizeIn2);
595
596 if (starPtr->invertFlag == 1) {
597
598 // This is a pair of open stars written in the same
599 // order as in this basis. Check preconditions:
600 UTIL_CHECK(starPtr2->invertFlag == -1);
601 UTIL_CHECK(starId2 = starId + 1);
602 UTIL_CHECK(basisId2 = basisId + 1);
603 UTIL_CHECK(starPtr->waveBz == waveIn);
604 UTIL_CHECK(starPtr2->waveBz == waveIn2);
605
606 // Copy components for both stars into fields array
607 for (int j = 0; j < nMonomer; ++j) {
608 fields[j][basisId] = temp[j];
609 fields[j][basisId2] = temp2[j];
610 }
611
612 } else
613 if (starPtr->invertFlag == -1) {
614
615 // This is a pair of open stars written in opposite
616 // order from in this basis. Check preconditions:
617 UTIL_CHECK(starPtr2->invertFlag == 1);
618 UTIL_CHECK(starId == starId2 + 1);
619 UTIL_CHECK(basisId == basisId2 + 1);
620 UTIL_CHECK(waveId == starPtr->beginId);
621
622 // Check that waveIn2 is negation of waveIn
623 IntVec<D> nVec;
624 nVec.negate(waveIn);
625 nVec =
626 shiftToMinimum(nVec, mesh.dimensions(), unitCell);
627 UTIL_CHECK(waveIn2 == nVec);
628
629 /*
630 * Consider two related stars, C and D, that are listed
631 * in the order (C,D) in the basis used in this code (the
632 * reading program), but that were listed in the opposite
633 * order (D,C) in the program that wrote the file (the
634 * writing program). In the basis of the reading program,
635 * star C has star index starId2, while star D has index
636 * starId = starid2 + 1.
637 *
638 * Let f(r) and f^{*}(r) denote the basis functions used
639 * by the reading program for stars C and D, respectively.
640 * Let u(r) and u^{*}(r) denote the corresponding basis
641 * functions used by the writing program for stars C
642 * and D. Let exp(i phi) denote the unit magnitude
643 * coefficient (i.e., phasor) within f(r) of the wave
644 * with wave index waveId2, which was the characteristic
645 * wave for star C in the writing program. The
646 * coefficient of this wave within the basis function
647 * u(r) used by the writing program must instead be real
648 * and positive. This implies that
649 * u(r) = exp(-i phi) f(r).
650 *
651 * Focus on the contribution to the field for a specific
652 * monomer type j. Let a and b denote the desired
653 * coefficients of stars C and D in the reading program,
654 * for which the total contribution of both stars to the
655 * field is:
656 *
657 * (a - ib) f(r) + (a + ib) f^{*}(r)
658 *
659 * Let A = temp[j] and B = temp2[j] denote the
660 * coefficients read from file in order (A,B). Noting
661 * that the stars were listed in the order (D,C) in the
662 * basis used by the writing program, the contribution
663 * of both stars must be (A-iB)u^{*}(r)+(A+iB)u(r), or:
664 *
665 * (A+iB) exp(-i phi)f(r) + (A-iB) exp(i phi) f^{*}(r)
666 *
667 * Comparing coefficients of f^{*}(r), we find that
668 *
669 * (a + ib) = (A - iB) exp(i phi)
670 *
671 * This equality is implemented below, where the
672 * variable "phasor" is set equal to exp(i phi).
673 */
674 phasor = basis.wave(waveId2).coeff;
675 phasor = phasor/std::abs(phasor);
676 for (int j = 0; j < nMonomer; ++j) {
677 coeff = std::complex<double>(temp[j],-temp2[j]);
678 coeff *= phasor;
679 fields[j][basisId2] = real(coeff);
680 fields[j][basisId ] = imag(coeff);
681 }
682
683 // Increment count of number of reversed open pairs
684 ++nReversedPair;
685
686 } else {
687 UTIL_THROW("Invalid starInvert value");
688 }
689
690 } // if (wavePtr->invertFlag == 0) ... else ...
691
692 } // if (waveExists && sizeMatches)
693
694 } // end while (i < nBasis)
695
696 if (nReversedPair > 0) {
697 Log::file() << "\n";
698 Log::file() << nReversedPair << " reversed pairs of open stars"
699 << " detected in readFieldsBasis\n";
700 }
701
702 }
703
704 /*
705 * Read the data section for a single field in basis format.
706 */
707 template <int D>
708 void readBasisData(std::istream& in,
709 DArray<double>& field,
710 UnitCell<D> const& unitCell,
711 Mesh<D> const& mesh,
712 Basis<D> const& basis,
713 int nBasisIn)
714 {
715 UTIL_CHECK(field.isAllocated());
716
717 // Local array of fields containing only one field
718 DArray< DArray<double> > fields;
719 fields.allocate(1);
720 fields[0].allocate(field.capacity());
721
722 // Call overloaded function that takes an array of fields
723 readBasisData(in, fields, unitCell, mesh, basis, nBasisIn);
724
725 // Copy data from local array fields to function parameter field
726 field = fields[0];
727 }
728
729 /*
730 * Write an array of fields in basis format to an output stream.
731 */
732 template <int D>
733 void writeBasisData(std::ostream &out,
734 DArray< DArray<double> > const & fields,
735 Basis<D> const & basis)
736 {
737 // Preconditions, set nMonomer and fieldCapacity
738 int nMonomer = fields.capacity();
739 UTIL_CHECK(nMonomer > 0);
740 int fieldCapacity = fields[0].capacity();
741 for (int j = 0; j < nMonomer; ++j) {
742 UTIL_CHECK(fieldCapacity == fields[j].capacity());
743 }
744 UTIL_CHECK(basis.isInitialized());
745
746 // Write fields
747 int ib = 0;
748 int nStar = basis.nStar();
749 for (int i = 0; i < nStar; ++i) {
750 if (ib >= fieldCapacity) break;
751 if (!basis.star(i).cancel) {
752 for (int j = 0; j < nMonomer; ++j) {
753 out << Dbl(fields[j][ib], 20, 10);
754 }
755 out << " ";
756 for (int j = 0; j < D; ++j) {
757 out << Int(basis.star(i).waveBz[j], 5);
758 }
759 out << Int(basis.star(i).size, 5) << std::endl;
760 ++ib;
761 }
762 }
763
764 }
765
766 /*
767 * Write data section of a single field in basis format.
768 */
769 template <int D>
770 void writeBasisData(std::ostream &out,
771 DArray<double> const & field,
772 Basis<D> const & basis)
773 {
774 UTIL_CHECK(field.isAllocated());
775
776 // Local array of fields containing only one field
777 DArray< DArray<double> > fields;
778 fields.allocate(1);
779 fields[0].allocate(field.capacity());
780 fields[0] = field;
781
782 // Call overloaded function that takes and array of fields
783 writeBasisData(out, fields, basis);
784 }
785
786 /*
787 * Convert an array of fields from basis to k-grid format.
788 */
789 template <int D, class ACT>
790 void convertBasisToKGrid(DArray<double> const & in,
791 ACT& out,
792 Basis<D> const& basis,
793 IntVec<D> const& dftDimensions)
794 {
795 UTIL_CHECK(basis.isInitialized());
796
797 typedef typename ACT::Complex CT;
798 typedef typename ACT::Real RT;
799
800 // Create Mesh<D> with dimensions of DFT Fourier grid.
801 Mesh<D> dftMesh(dftDimensions);
802
803 typename Basis<D>::Star const* starPtr; // pointer to current star
804 typename Basis<D>::Wave const* wavePtr; // pointer to current wave
805 std::complex<double> component; // coefficient for star
806 std::complex<double> coeff; // coefficient for wave
807 IntVec<D> indices; // dft grid indices of wave
808 int rank; // dft grid rank of wave
809 int is; // star index
810 int ib; // basis index
811 int iw; // wave index
812
813 // Initialize all dft coponents to zero
814 for (rank = 0; rank < dftMesh.size(); ++rank) {
815 assign<CT, RT>(out[rank], 0.0, 0.0);
816 }
817
818 // Loop over stars, skipping cancelled stars
819 is = 0;
820 while (is < basis.nStar()) {
821 starPtr = &(basis.star(is));
822
823 if (starPtr->cancel) {
824 ++is;
825 continue;
826 }
827
828 // Set basisId for uncancelled star
829 ib = starPtr->basisId;
830
831 if (starPtr->invertFlag == 0) {
832
833 // Make complex coefficient for star basis function
834 component = std::complex<double>(in[ib], 0.0);
835
836 // Loop over waves in closed star
837 for (iw = starPtr->beginId; iw < starPtr->endId; ++iw) {
838 wavePtr = &basis.wave(iw);
839 if (!wavePtr->implicit) {
840 coeff = component*(wavePtr->coeff);
841 indices = wavePtr->indicesDft;
842 rank = dftMesh.rank(indices);
843 //out[rank][0] = coeff.real();
844 //out[rank][1] = coeff.imag();
845 assign<CT, RT>(out[rank], coeff);
846 }
847 }
848 ++is;
849
850 } else
851 if (starPtr->invertFlag == 1) {
852
853 // Loop over waves in first star
854 component = std::complex<double>(in[ib], -in[ib+1]);
855 component /= sqrt(2.0);
856 starPtr = &(basis.star(is));
857 for (iw = starPtr->beginId; iw < starPtr->endId; ++iw) {
858 wavePtr = &basis.wave(iw);
859 if (!(wavePtr->implicit)) {
860 coeff = component*(wavePtr->coeff);
861 indices = wavePtr->indicesDft;
862 rank = dftMesh.rank(indices);
863 //out[rank][0] = coeff.real();
864 //out[rank][1] = coeff.imag();
865 assign<CT, RT>(out[rank], coeff);
866 }
867 }
868
869 // Loop over waves in second star
870 starPtr = &(basis.star(is+1));
871 UTIL_CHECK(starPtr->invertFlag == -1);
872 component = std::complex<double>(in[ib], +in[ib+1]);
873 component /= sqrt(2.0);
874 for (iw = starPtr->beginId; iw < starPtr->endId; ++iw) {
875 wavePtr = &basis.wave(iw);
876 if (!(wavePtr->implicit)) {
877 coeff = component*(wavePtr->coeff);
878 indices = wavePtr->indicesDft;
879 rank = dftMesh.rank(indices);
880 //out[rank][0] = coeff.real();
881 //out[rank][1] = coeff.imag();
882 assign<CT, RT>(out[rank], coeff);
883 }
884 }
885
886 // Increment is by 2 (two stars were processed)
887 is += 2;
888
889 } else {
890
891 UTIL_THROW("Invalid invertFlag value");
892
893 }
894
895 }
896
897 }
898
899 template <int D, class ACT>
900 void convertKGridToBasis(ACT const& in,
901 DArray<double>& out,
902 Basis<D> const& basis,
903 IntVec<D> const& dftDimensions,
904 bool checkSymmetry,
905 double epsilon)
906 {
907 UTIL_CHECK(basis.isInitialized());
908
909 typedef typename ACT::Complex CT;
910 typedef typename ACT::Real RT;
911
912 // Check if input field in k-grid format has specified symmetry
913 if (checkSymmetry) {
914 bool symmetric;
915 symmetric = hasSymmetry(in, basis, dftDimensions, epsilon, true);
916 if (!symmetric) {
917 Log::file() << std::endl
918 << "WARNING: Conversion of asymmetric field to"
919 << "symmetry-adapted basis in Prdc::convertKgridToBasis."
920 << std::endl << std::endl;
921 }
922 }
923
924 // Create Mesh<D> with dimensions of DFT Fourier grid.
925 Mesh<D> dftMesh(dftDimensions);
926
927 typename Basis<D>::Star const* starPtr; // pointer to current star
928 typename Basis<D>::Wave const* wavePtr; // pointer to current wave
929 std::complex<double> component; // coefficient for star
930 int rank; // dft grid rank of wave
931 int is; // star index
932 int ib; // basis index
933 int iw; // wave index
934
935 // Initialize all components to zero
936 for (is = 0; is < basis.nBasis(); ++is) {
937 out[is] = 0.0;
938 }
939
940 // Loop over stars
941 is = 0;
942 while (is < basis.nStar()) {
943 starPtr = &(basis.star(is));
944
945 if (starPtr->cancel) {
946 ++is;
947 continue;
948 }
949
950 // Set basis id for uncancelled star
951 ib = starPtr->basisId;
952
953 if (starPtr->invertFlag == 0) {
954
955 // Choose a wave in the star that is not implicit
956 int beginId = starPtr->beginId;
957 int endId = starPtr->endId;
958 iw = 0;
959 bool isImplicit = true;
960 while (isImplicit) {
961 wavePtr = &basis.wave(beginId + iw);
962 if (!wavePtr->implicit) {
963 isImplicit = false;
964 } else {
965 UTIL_CHECK(beginId + iw < endId - 1 - iw);
966 wavePtr = &basis.wave(endId - 1 - iw);
967 if (!wavePtr->implicit) {
968 isImplicit = false;
969 }
970 }
971 ++iw;
972 }
973 UTIL_CHECK(wavePtr->starId == is);
974
975 // Compute component value
976 rank = dftMesh.rank(wavePtr->indicesDft);
977 //component = std::complex<double>(in[rank][0], in[rank][1]);
978 assign<CT,RT>(component, in[rank]);
979 component /= wavePtr->coeff;
980 out[ib] = component.real();
981 ++is;
982
983 } else
984 if (starPtr->invertFlag == 1) {
985
986 // Identify a characteristic wave that is not implicit:
987 // Either the first wave of the 1st star or its inverse
988 // in the second star
989 wavePtr = &basis.wave(starPtr->beginId);
990 UTIL_CHECK(wavePtr->starId == is);
991 if (wavePtr->implicit) {
992 iw = wavePtr->inverseId;
993 starPtr = &(basis.star(is+1));
994 UTIL_CHECK(starPtr->invertFlag == -1);
995 wavePtr = &basis.wave(iw);
996 UTIL_CHECK(!(wavePtr->implicit));
997 UTIL_CHECK(wavePtr->starId == is+1);
998 }
999 rank = dftMesh.rank(wavePtr->indicesDft);
1000 //component = std::complex<double>(in[rank][0], in[rank][1]);
1001 assign<CT, RT>(component, in[rank]);
1002 UTIL_CHECK(std::abs(wavePtr->coeff) > 1.0E-8);
1003 component /= wavePtr->coeff;
1004 component *= sqrt(2.0);
1005
1006 // Compute basis function coefficient values
1007 if (starPtr->invertFlag == 1) {
1008 out[ib] = component.real();
1009 out[ib+1] = -component.imag();
1010 } else {
1011 out[ib] = component.real();
1012 out[ib+1] = component.imag();
1013 }
1014
1015 is += 2;
1016 } else {
1017 UTIL_THROW("Invalid invertFlag value");
1018 }
1019
1020 } // loop over star index is
1021 }
1022
1023 /*
1024 * Test if an K-grid array has the declared space group symmetry.
1025 * Return true if symmetric, false otherwise. Print error values
1026 * if verbose == true and hasSymmetry == false.
1027 */
1028 template <int D, class ACT>
1029 bool hasSymmetry(ACT const & in,
1030 Basis<D> const& basis,
1031 IntVec<D> const& dftDimensions,
1032 double epsilon,
1033 bool verbose)
1034 {
1035 UTIL_CHECK(basis.isInitialized());
1036
1037 typedef typename ACT::Complex CT;
1038 typedef typename ACT::Real RT;
1039
1040 typename Basis<D>::Star const* starPtr; // pointer to current star
1041 typename Basis<D>::Wave const* wavePtr; // pointer to current wave
1042 std::complex<double> waveCoeff; // coefficient from wave
1043 std::complex<double> rootCoeff; // coefficient from root
1044 std::complex<double> diff; // coefficient difference
1045 int is; // star index
1046 int iw; // wave index
1047 int beginId, endId; // star begin, end ids
1048 int rank; // dft grid rank of wave
1049
1050 double cancelledError(0.0); // max error from cancelled stars
1051 double uncancelledError(0.0); // max error from uncancelled stars
1052
1053 // Create Mesh<D> with dimensions of DFT Fourier grid.
1054 Mesh<D> dftMesh(dftDimensions);
1055
1056 // Loop over all stars
1057 for (is = 0; is < basis.nStar(); ++is) {
1058 starPtr = &(basis.star(is));
1059
1060 if (starPtr->cancel) {
1061
1062 // Check that coefficients are zero for all waves in star
1063 beginId = starPtr->beginId;
1064 endId = starPtr->endId;
1065 for (iw = beginId; iw < endId; ++iw) {
1066 wavePtr = &basis.wave(iw);
1067 if (!wavePtr->implicit) {
1068 rank = dftMesh.rank(wavePtr->indicesDft);
1069 //waveCoeff
1070 // = std::complex<double>(in[rank][0], in[rank][1]);
1071 assign<CT, RT>(waveCoeff, in[rank]);
1072 if (std::abs(waveCoeff) > cancelledError) {
1073 cancelledError = std::abs(waveCoeff);
1074 if ((!verbose) && (cancelledError > epsilon)) {
1075 return false;
1076 }
1077 }
1078 }
1079 }
1080
1081 } else {
1082
1083 // Check consistency of coeff values from all waves
1084 bool hasRoot = false;
1085 beginId = starPtr->beginId;
1086 endId = starPtr->endId;
1087 for (iw = beginId; iw < endId; ++iw) {
1088 wavePtr = &basis.wave(iw);
1089 if (!(wavePtr->implicit)) {
1090 rank = dftMesh.rank(wavePtr->indicesDft);
1091 //waveCoeff
1092 // = std::complex<double>(in[rank][0], in[rank][1]);
1093 assign<CT, RT>(waveCoeff, in[rank]);
1094 waveCoeff /= wavePtr->coeff;
1095 if (hasRoot) {
1096 diff = waveCoeff - rootCoeff;
1097 if (std::abs(diff) > uncancelledError) {
1098 uncancelledError = std::abs(diff);
1099 if ((!verbose) && (uncancelledError > epsilon)) {
1100 return false;
1101 }
1102 }
1103 } else {
1104 rootCoeff = waveCoeff;
1105 hasRoot = true;
1106 }
1107 }
1108 }
1109
1110 }
1111
1112 } // end loop over star index is
1113
1114 if ((cancelledError < epsilon) && (uncancelledError < epsilon)) {
1115 return true;
1116 } else if (verbose) {
1117 Log::file()
1118 << std::endl
1119 << "Maximum coefficient of a cancelled star: "
1120 << cancelledError << std::endl
1121 << "Maximum error of coefficient for an uncancelled star: "
1122 << uncancelledError << std::endl;
1123 }
1124 return false;
1125 }
1126
1127 // Field file manipulations
1128
1129 template <int D, class ART>
1130 void replicateUnitCell(std::ostream &out,
1131 DArray<ART> const & fields,
1132 IntVec<D> const & meshDimensions,
1133 UnitCell<D> const & unitCell,
1134 IntVec<D> const & replicas)
1135 {
1136 // Obtain number of monomer types
1137 int nMonomer = fields.capacity();
1138 UTIL_CHECK(nMonomer > 0);
1139
1140 // Compute properties of mesh for replicated fields
1141 IntVec<D> repDimensions; // Dimensions
1142 int repSize = 1; // Total number of grid points
1143 for (int i = 0; i < D; ++i) {
1144 UTIL_CHECK(replicas[i] > 0);
1145 UTIL_CHECK(meshDimensions[i] > 0);
1146 repDimensions[i] = replicas[i] * meshDimensions[i];
1147 repSize *= repDimensions[i];
1148 }
1149
1150 // Allocate arrays for replicated fields
1151 DArray< DArray<double> > repFields;
1152 repFields.allocate(nMonomer);
1153 for (int i = 0; i < nMonomer; ++i) {
1154 repFields[i].allocate(repSize);
1155 }
1156
1157 // Replicate data
1158 Mesh<D> mesh(meshDimensions);
1159 Mesh<D> cellMesh(replicas);
1160 Mesh<D> repMesh(repDimensions);
1161 MeshIterator<D> iter(meshDimensions);
1162 MeshIterator<D> cellIter(replicas);
1163 IntVec<D> position;
1164 IntVec<D> cellPosition;
1165 IntVec<D> repPosition;
1166 double value;
1167 int repRank;
1168 for (int i = 0; i < nMonomer; ++i) {
1169 for (iter.begin(); !iter.atEnd(); ++iter) {
1170 position = iter.position();
1171 value = fields[i][iter.rank()];
1172 for (cellIter.begin(); !cellIter.atEnd(); ++cellIter) {
1173 cellPosition = cellIter.position();
1174 for (int j=0; j < D; ++j) {
1175 repPosition[j] = position[j]
1176 + meshDimensions[j]*cellPosition[j];
1177 }
1178 repRank = repMesh.rank(repPosition);
1179 repFields[i][repRank] = value;
1180 }
1181 }
1182 }
1183
1184 // Create a replicated UnitCell<D> object
1185 // Use function declared in prdc/crystal/replicateUnitCell.h
1186 UnitCell<D> cell;
1187 replicateUnitCell(replicas, unitCell, cell);
1188
1189 // Write header
1190 int v1 = 1;
1191 int v2 = 0;
1192 std::string gName = "";
1193 Pscf::Prdc::writeFieldHeader(out, v1, v2, cell, gName, nMonomer);
1194 writeMeshDimensions(out, repDimensions);
1195
1196 // Write field data
1197 writeRGridData(out, repFields, nMonomer, repDimensions);
1198 }
1199
1200 template <int D, class ART>
1201 void
1202 expandRGridDimension(std::ostream &out,
1203 DArray<ART> const & fields,
1204 IntVec<D> const & meshDimensions,
1205 UnitCell<D> const & unitCell,
1206 int d,
1207 DArray<int> newGridDimensions)
1208 {
1209 UTIL_THROW("Unimplemented base template");
1210 }
1211
1212 template <class ART>
1213 void
1214 expandRGridDimension(std::ostream &out,
1215 DArray<ART> const & fields,
1216 IntVec<1> const & meshDimensions,
1217 UnitCell<1> const & unitCell,
1218 int d,
1219 DArray<int> newGridDimensions)
1220
1221 {
1222 // Check validity of expanded dimension d and newGridDimensions
1223 UTIL_CHECK(d > 1);
1224 UTIL_CHECK(d <= 3);
1225 UTIL_CHECK(newGridDimensions.capacity() == (d - 1));
1226
1227
1228 // Obtain number of monomer types
1229 int nMonomer = fields.capacity();
1230 UTIL_CHECK(nMonomer > 0);
1231
1232 // Set up necessary objects
1233 FSArray<double, 6> cellParameters;
1234 cellParameters.append(unitCell.parameter(0));
1235 int v1 = 1;
1236 int v2 = 0;
1237 std::string gName = "";
1238
1239 if (d == 2) {
1240
1241 // 1D expanded to 2D
1242
1243 // Set dimensions
1244 IntVec<2> dimensions;
1245 dimensions[0] = meshDimensions[0];
1246 dimensions[1] = newGridDimensions[0];
1247
1248 // Assign unit cell
1249 UnitCell<2> cell;
1250 if (dimensions[0] == dimensions[1]) {
1251 cell.set(UnitCell<2>::Square, cellParameters);
1252 } else {
1253 cellParameters.append((double)dimensions[1]/dimensions[0]
1254 * cellParameters[0]);
1255 cell.set(UnitCell<2>::Rectangular, cellParameters);
1256 }
1257
1258 // Write Header
1259 Pscf::Prdc::writeFieldHeader(out, v1, v2, cell, gName, nMonomer);
1260 out << "mesh " << std::endl
1261 << " " << dimensions << std::endl;
1262
1263 } else if (d == 3) {
1264
1265 // 1D expanded to 3D
1266
1267 // Set dimensions
1268 IntVec<3> dimensions;
1269 dimensions[0] = meshDimensions[0];
1270 dimensions[1] = newGridDimensions[0];
1271 dimensions[2] = newGridDimensions[1];
1272
1273 // Assign unit cell
1274 UnitCell<3> cell;
1275 if (dimensions[2] == dimensions[1]) {
1276 if (dimensions[1] == dimensions[0]) {
1277 cell.set(UnitCell<3>::Cubic, cellParameters);
1278 } else {
1279 cellParameters.append((double)dimensions[1]/dimensions[0]
1280 * cellParameters[0]);
1281 cellParameters.append((double)dimensions[2]/dimensions[0]
1282 * cellParameters[0]);
1283 cell.set(UnitCell<3>::Orthorhombic, cellParameters);
1284 }
1285 } else {
1286 if (dimensions[1] == dimensions[0]) {
1287 cellParameters.append((double)dimensions[2]/dimensions[0]
1288 * cellParameters[0]);
1289 cell.set(UnitCell<3>::Tetragonal, cellParameters);
1290 } else {
1291 cellParameters.append((double)dimensions[1]/dimensions[0]
1292 * cellParameters[0]);
1293 cellParameters.append((double)dimensions[2]/dimensions[0]
1294 * cellParameters[0]);
1295 cell.set(UnitCell<3>::Orthorhombic, cellParameters);
1296 }
1297 }
1298
1299 // Write Header
1300 Pscf::Prdc::writeFieldHeader(out, v1, v2, cell, gName, nMonomer);
1301 out << "mesh " << std::endl
1302 << " " << dimensions << std::endl;
1303
1304 } else {
1305
1306 UTIL_THROW("Invalid d value");
1307
1308 }
1309
1310 // Write expanded fields
1311 int nReplica = newGridDimensions[0];
1312 if (d == 3) {
1313 nReplica *= newGridDimensions[1];
1314 }
1315 MeshIteratorFortran<1> iter(meshDimensions);
1316 int rank;
1317 for (int i = 0; i < nReplica; ++i) {
1318 for (iter.begin(); !iter.atEnd(); ++iter) {
1319 rank = iter.rank();
1320 for (int j = 0; j < nMonomer; ++j) {
1321 out << " " << Dbl(fields[j][rank], 21, 13);
1322 }
1323 out << std::endl;
1324 }
1325 }
1326
1327 }
1328
1329 template <class ART>
1330 void expandRGridDimension(std::ostream &out,
1331 DArray<ART> const & fields,
1332 IntVec<2> const & meshDimensions,
1333 UnitCell<2> const & unitCell,
1334 int d,
1335 DArray<int> newGridDimensions)
1336 {
1337 // 2D expanded to 3D
1338
1339 // Check validity of expanded dimension d and newGridDimensions
1340 UTIL_CHECK(d == 3);
1341 UTIL_CHECK(newGridDimensions.capacity() == 1);
1342
1343 // Obtain number of monomer types
1344 int nMonomer = fields.capacity();
1345 UTIL_CHECK(nMonomer > 0);
1346
1347 // Set dimensions of mesh for replicated fields
1348 IntVec<3> dimensions;
1349 dimensions[0] = meshDimensions[0];
1350 dimensions[1] = meshDimensions[1];
1351 dimensions[2] = newGridDimensions[0];
1352
1353 // Set unit cell for replicated fields
1354 UnitCell<3> cell;
1355 FSArray<double, 6> cellParameters;
1356 cellParameters.append(unitCell.parameter(0));
1357 if (unitCell.lattice() == UnitCell<2>::Square) {
1358 if (newGridDimensions[0] == meshDimensions[0]){
1359 cell.set(UnitCell<3>::Cubic, cellParameters);
1360 } else {
1361 cellParameters.append((double)dimensions[2]/dimensions[0]
1362 * cellParameters[0]);
1363 cell.set(UnitCell<3>::Tetragonal, cellParameters);
1364 }
1365 } else if (unitCell.lattice() == UnitCell<2>::Rectangular) {
1366 cellParameters.append(unitCell.parameter(1));
1367 cellParameters.append((double)dimensions[2]/dimensions[0]
1368 * cellParameters[0]);
1369 cell.set(UnitCell<3>::Orthorhombic, cellParameters);
1370 } else if (unitCell.lattice() == UnitCell<2>::Hexagonal){
1371 cellParameters.append((double)dimensions[2]/dimensions[0]
1372 * cellParameters[0]);
1373 cell.set(UnitCell<3>::Hexagonal, cellParameters);
1374 } else if (unitCell.lattice() == UnitCell<2>::Rhombic) {
1375 cellParameters.append(unitCell.parameter(0));
1376 cellParameters.append((double)dimensions[2]/dimensions[0]
1377 * cellParameters[0]);
1378 cellParameters.append(Constants::Pi / 2);
1379 cellParameters.append(0.0);
1380 cellParameters.append(unitCell.parameter(1));
1381 cell.set(UnitCell<3>::Triclinic, cellParameters);
1382 } else if (unitCell.lattice() == UnitCell<2>::Oblique) {
1383 cellParameters.append(unitCell.parameter(1));
1384 cellParameters.append((double)dimensions[2]/dimensions[0]
1385 * cellParameters[0]);
1386 cellParameters.append(Constants::Pi / 2);
1387 cellParameters.append(0.0);
1388 cellParameters.append(unitCell.parameter(2));
1389 cell.set(UnitCell<3>::Triclinic, cellParameters);
1390 } else {
1391 UTIL_THROW("Unrecognized 2D lattice system.");
1392 }
1393
1394 // Write Header
1395 int v1 = 1;
1396 int v2 = 0;
1397 std::string gName = "";
1398 Pscf::Prdc::writeFieldHeader(out, v1, v2, cell, gName, nMonomer);
1399 out << "mesh " << std::endl
1400 << " " << dimensions << std::endl;
1401
1402 // Write expanded fields
1403 int nReplica = newGridDimensions[0];
1404 MeshIteratorFortran<2> iter(meshDimensions);
1405 int rank;
1406 for (int i = 0; i < nReplica; ++i) {
1407 for (iter.begin(); !iter.atEnd(); ++iter) {
1408 rank = iter.rank();
1409 for (int j = 0; j < nMonomer; ++j) {
1410 out << " " << Dbl(fields[j][rank], 21, 13);
1411 }
1412 out << std::endl;
1413 }
1414 }
1415
1416 }
1417
1418 template <class ART>
1419 void expandRGridDimension(std::ostream &out,
1420 DArray<ART> const & fields,
1421 IntVec<3> const & meshDimensions,
1422 UnitCell<3> const & unitCell,
1423 int d,
1424 DArray<int> newGridDimensions)
1425
1426 { UTIL_THROW("expandRGridDimension is invalid when D = 3."); }
1427
1428} // namespace Prdc
1429} // namespace Pscf
1430#endif
An IntVec<D, T> is a D-component vector of elements of integer type T.
Definition IntVec.h:27
Iterator over points in a mesh in "Fortran" order.
Iterator over points in a Mesh<D>.
Description of a regular grid of points in a periodic domain.
Definition Mesh.h:61
IntVec< D > dimensions() const
Get an IntVec<D> of the grid dimensions.
Definition Mesh.h:217
int shift(int &coordinate, int i) const
Shift a periodic coordinate into range.
Definition Mesh.tpp:131
A list of wavevectors that are related by space-group symmetries.
Definition Basis.h:482
int basisId
Index of basis function associated with this star.
Definition Basis.h:563
IntVec< D > waveBz
Integer indices indicesBz of a characteristic wave of this star.
Definition Basis.h:544
int invertFlag
Index for inversion symmetry of star.
Definition Basis.h:531
int size
Number of wavevectors in this star.
Definition Basis.h:489
int beginId
Wave index of first wavevector in star.
Definition Basis.h:494
bool cancel
Is this star cancelled, i.e., associated with a zero function?
Definition Basis.h:578
Wavevector used to construct a basis function.
Definition Basis.h:391
std::complex< double > coeff
Coefficient of wave within the associated star basis function.
Definition Basis.h:398
int starId
Index of the star that contains this wavevector.
Definition Basis.h:427
Symmetry-adapted Fourier basis for pseudo-spectral SCFT.
Definition Basis.h:383
Wave const & wave(int id) const
Get a specific Wave, access by integer index.
Definition Basis.h:876
int waveId(IntVec< D > vector) const
Get the integer index of a wave, as required by wave(int id).
Definition Basis.h:890
int nStar() const
Total number of stars.
Definition Basis.h:871
bool isInitialized() const
Returns true iff this basis is fully initialized.
Definition Basis.h:898
Star const & star(int id) const
Get a Star, accessed by integer star index.
Definition Basis.h:881
Base template for UnitCell<D> classes, D=1, 2 or 3.
Definition UnitCell.h:56
Vec< D, T > & negate(const Vec< D, T > &v)
Return negative of vector v.
Definition Vec.h:520
int capacity() const
Return allocated size.
Definition Array.h:159
static const double Pi
Trigonometric constant Pi.
Definition Constants.h:35
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
Wrapper for a double precision number, for formatted ostream output.
Definition Dbl.h:40
A fixed capacity (static) contiguous array with a variable logical size.
Definition FSArray.h:38
void append(Data const &data)
Append data to the end of the array.
Definition FSArray.h:258
Wrapper for an int, for formatted ostream output.
Definition Int.h:37
static std::ostream & file()
Get log ostream by reference.
Definition Log.cpp:59
#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 readUnitCellHeader(std::istream &in, UnitCell< D > &cell)
Read UnitCell<D> from a field file header (fortran PSCF format).
void replicateUnitCell(IntVec< 1 > const &replicas, UnitCell< 1 > const &cellIn, UnitCell< 1 > &cellOut)
Create a replicated UnitCell<1>.
void writeFieldHeader(std::ostream &out, int ver1, int ver2, UnitCell< D > const &cell, std::string const &groupName, int nMonomer)
Write common part of field header (fortran PSCF format).
void writeUnitCellHeader(std::ostream &out, UnitCell< D > const &cell)
Write UnitCell<D> to a field file header (fortran PSCF format).
void readFieldHeader(std::istream &in, int &ver1, int &ver2, UnitCell< D > &cell, std::string &groupName, int &nMonomer)
Read common part of field header (fortran PSCF format).
void writeKGridData(std::ostream &out, DArray< AT > const &fields, int nMonomer, IntVec< D > const &dftDimensions)
Write data for array of k-grid fields, with no header section.
bool hasSymmetry(AT const &in, Basis< D > const &basis, IntVec< D > const &dftDimensions, double epsilon=1.0e-8, bool verbose=true)
Check if a k-grid field has the declared space group symmetry.
void inspectArrays(DArray< AT > const &arrays, int &nMonomer, int &capacity)
Inspect dimensions of a DArray of 1D arrays, each of type AT.
void inspectFields(DArray< FT > const &fields, int &nMonomer, IntVec< D > &dimensions)
Inspect dimensions of a DArray of fields, each of type FT.
void checkAllocateArrays(DArray< AT > &arrays, int nMonomer, int capacity)
Check allocation of a DArray of 1D arrays, allocate if necessary.
void readKGridData(std::istream &in, DArray< AT > &fields, int nMonomer, IntVec< D > const &dftDimensions)
Read data for array of k-grid fields, with no header section.
void convertBasisToKGrid(DArray< double > const &components, AT &dft, Basis< D > const &basis, IntVec< D > const &dftDimensions)
Convert a real field from symmetrized basis to Fourier grid.
void writeMeshDimensions(std::ostream &out, IntVec< D > const &meshDimensions)
Write mesh dimensions to a field file header.
void expandRGridDimension(std::ostream &out, DArray< AT > const &fields, IntVec< D > const &meshDimensions, UnitCell< D > const &unitCell, int d, DArray< int > newGridDimensions)
Expand the dimensionality of space from D to d.
void writeRGridData(std::ostream &out, DArray< AT > const &fields, int nMonomer, IntVec< D > const &dimensions)
Write data for array of r-grid fields, with no header section.
void convertKGridToBasis(AT const &in, DArray< double > &out, Basis< D > const &basis, IntVec< D > const &dftDimensions, bool checkSymmetry=true, double epsilon=1.0e-8)
Convert a real field from Fourier grid to symmetrized basis.
void readRGridData(std::istream &in, DArray< AT > &fields, int nMonomer, IntVec< D > const &dimensions)
Read data for array of r-grid fields, with no header section.
void readMeshDimensions(std::istream &in, IntVec< D > const &meshDimensions)
Read mesh dimensions from a field file header.
void writeBasisData(std::ostream &out, DArray< DArray< double > > const &fields, Basis< D > const &basis)
Write data section of an array of fields in basis format.
void readBasisData(std::istream &in, DArray< DArray< double > > &fields, UnitCell< D > const &unitCell, Mesh< D > const &mesh, Basis< D > const &basis, int nStarIn)
Read an array of fields in basis format, without a header.
void checkAllocateFields(DArray< FT > &fields, int nMonomer, IntVec< D > const &dimensions)
Check allocation of an array of fields, allocate if necessary.
void checkAllocateField(FT &field, IntVec< D > const &dimensions)
Check allocation of a single field, allocate if necessary.
RT imag(CT const &a)
Return the imaginary part of a complex number.
void assign(CT &z, RT const &a, RT const &b)
Create a complex number from real and imaginary parts, z = a + ib.
RT real(CT const &a)
Return the real part of a complex number.
Periodic fields and crystallography.
Definition CField.cpp:11
PSCF package top-level namespace.
Definition param_pc.dox:1