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