13#include <rpc/fts/simulator/Simulator.h>
14#include <rpc/fts/simulator/SimulatorFactory.h>
15#include <rpc/fts/compressor/Compressor.h>
16#include <rpc/scft/sweep/Sweep.h>
17#include <rpc/scft/sweep/SweepFactory.h>
18#include <rpc/scft/iterator/Iterator.h>
19#include <rpc/scft/iterator/IteratorFactory.h>
20#include <rpc/solvers/Polymer.h>
21#include <rpc/solvers/Solvent.h>
24#include <prdc/cpu/RField.h>
25#include <prdc/cpu/RFieldComparison.h>
26#include <prdc/crystal/BFieldComparison.h>
28#include <pscf/inter/Interaction.h>
29#include <pscf/math/IntVec.h>
30#include <pscf/homogeneous/Clump.h>
32#include <util/containers/FSArray.h>
33#include <util/param/BracketPolicy.h>
34#include <util/param/ParamComponent.h>
35#include <util/format/Str.h>
36#include <util/format/Int.h>
37#include <util/format/Dbl.h>
38#include <util/misc/ioUtil.h>
59 interactionPtr_(nullptr),
60 iteratorPtr_(nullptr),
61 iteratorFactoryPtr_(nullptr),
63 sweepFactoryPtr_(nullptr),
64 simulatorPtr_(nullptr),
65 simulatorFactoryPtr_(nullptr),
75 isAllocatedGrid_(false),
76 isAllocatedBasis_(false),
82 domain_.setFileMaster(fileMaster_);
83 w_.setFieldIo(
domain().fieldIo());
84 h_.setFieldIo(
domain().fieldIo());
85 mask_.setFieldIo(
domain().fieldIo());
101 if (interactionPtr_) {
102 delete interactionPtr_;
107 if (iteratorFactoryPtr_) {
108 delete iteratorFactoryPtr_;
113 if (sweepFactoryPtr_) {
114 delete sweepFactoryPtr_;
117 delete simulatorPtr_;
119 if (simulatorFactoryPtr_) {
120 delete simulatorFactoryPtr_;
147 while ((c = getopt(argc, argv,
"ed:p:c:i:o:t:")) != -1) {
177 Log::file() <<
"Unknown option -" << optopt << std::endl;
196 fileMaster_.setParamFileName(std::string(pArg));
198 UTIL_THROW(
"Missing required -p option - no parameter file");
203 fileMaster_.setCommandFileName(std::string(cArg));
205 UTIL_THROW(
"Missing required -c option - no command file");
210 fileMaster_.setInputPrefix(std::string(iArg));
215 fileMaster_.setOutputPrefix(std::string(oArg));
220 UTIL_THROW(
"Error: Non-positive thread count -t option");
233 readParamComposite(in, mixture_);
236 int nm = mixture().nMonomer();
237 int np = mixture().nPolymer();
238 int ns = mixture().nSolvent();
245 interaction().setNMonomer(nm);
246 readParamComposite(in, interaction());
249 readParamComposite(in, domain_);
251 UTIL_CHECK(domain().unitCell().nParameter() > 0);
255 mixture_.associate(domain().mesh(), domain().fft(),
256 domain().unitCell());
258 mixture_.clearUnitCellData();
261 allocateFieldsGrid();
264 std::string className;
267 iteratorFactoryPtr_->readObjectOptional(in, *
this, className,
270 Log::file() << indent() <<
" Iterator{ [absent] }\n";
274 if (iteratorPtr_ && !isEnd) {
276 sweepFactoryPtr_->readObjectOptional(in, *
this, className,
279 Log::file() << indent() <<
" Sweep{ [absent] }\n";
286 simulatorFactoryPtr_->readObjectOptional(in, *
this,
289 Log::file() << indent() <<
" Simulator{ [absent] }\n";
295 homogeneous_.setNMolecule(np+ns);
296 homogeneous_.setNMonomer(nm);
307 readBegin(in, className().c_str());
317 { readParam(fileMaster_.paramFile()); }
326 std::string command, filename, inFileName, outFileName;
328 bool readNext =
true;
339 if (command ==
"FINISH") {
343 if (command ==
"READ_W_BASIS") {
344 readEcho(in, filename);
345 readWBasis(filename);
347 if (command ==
"READ_W_RGRID") {
348 readEcho(in, filename);
349 readWRGrid(filename);
351 if (command ==
"ESTIMATE_W_FROM_C") {
352 readEcho(in, inFileName);
353 estimateWfromC(inFileName);
355 if (command ==
"SET_UNIT_CELL") {
358 Log::file() <<
" " << unitCell << std::endl;
359 setUnitCell(unitCell);
361 if (command ==
"COMPUTE") {
365 if (command ==
"ITERATE") {
367 bool isContinuation =
false;
368 int fail = iterate(isContinuation);
373 if (command ==
"SWEEP") {
378 if (command ==
"COMPRESS") {
381 simulator().compressor().compress();
383 if (command ==
"SIMULATE") {
390 if (command ==
"ANALYZE" || command ==
"ANALYZE_TRAJECTORY") {
397 std::string classname;
398 readEcho(in, classname);
399 readEcho(in, filename);
400 simulator().analyze(min, max, classname, filename);
402 if (command ==
"WRITE_TIMERS") {
403 readEcho(in, filename);
405 fileMaster_.openOutputFile(filename, file);
409 if (command ==
"CLEAR_TIMERS") {
412 if (command ==
"WRITE_PARAM") {
413 readEcho(in, filename);
415 fileMaster_.openOutputFile(filename, file);
416 writeParamNoSweep(file);
419 if (command ==
"WRITE_THERMO") {
420 readEcho(in, filename);
422 fileMaster_.openOutputFile(filename, file,
427 if (command ==
"WRITE_STRESS") {
428 readEcho(in, filename);
430 fileMaster_.openOutputFile(filename, file,
435 if (command ==
"WRITE_W_BASIS") {
436 readEcho(in, filename);
437 writeWBasis(filename);
439 if (command ==
"WRITE_W_RGRID") {
440 readEcho(in, filename);
441 writeWRGrid(filename);
443 if (command ==
"WRITE_C_BASIS") {
444 readEcho(in, filename);
445 writeCBasis(filename);
447 if (command ==
"WRITE_C_RGRID") {
448 readEcho(in, filename);
449 writeCRGrid(filename);
451 if (command ==
"WRITE_BLOCK_C_RGRID") {
452 readEcho(in, filename);
453 writeBlockCRGrid(filename);
455 if (command ==
"WRITE_Q_SLICE") {
456 int polymerId, blockId, directionId, segmentId;
457 readEcho(in, filename);
462 Log::file() <<
Str(
"polymer ID ", 21) << polymerId <<
"\n"
463 <<
Str(
"block ID ", 21) << blockId <<
"\n"
464 <<
Str(
"direction ID ", 21) << directionId <<
"\n"
465 <<
Str(
"segment ID ", 21) << segmentId
467 writeQSlice(filename, polymerId, blockId, directionId,
470 if (command ==
"WRITE_Q_TAIL") {
471 readEcho(in, filename);
472 int polymerId, blockId, directionId;
476 Log::file() <<
Str(
"polymer ID ", 21) << polymerId <<
"\n"
477 <<
Str(
"block ID ", 21) << blockId <<
"\n"
478 <<
Str(
"direction ID ", 21) << directionId
480 writeQTail(filename, polymerId, blockId, directionId);
482 if (command ==
"WRITE_Q") {
483 readEcho(in, filename);
484 int polymerId, blockId, directionId;
488 Log::file() <<
Str(
"polymer ID ", 21) << polymerId <<
"\n"
489 <<
Str(
"block ID ", 21) << blockId <<
"\n"
490 <<
Str(
"direction ID ", 21) << directionId
492 writeQ(filename, polymerId, blockId, directionId);
494 if (command ==
"WRITE_Q_ALL") {
495 readEcho(in, filename);
498 if (command ==
"WRITE_STARS") {
499 readEcho(in, filename);
500 writeStars(filename);
502 if (command ==
"WRITE_WAVES") {
503 readEcho(in, filename);
504 writeWaves(filename);
506 if (command ==
"WRITE_GROUP") {
507 readEcho(in, filename);
508 writeGroup(filename);
510 if (command ==
"BASIS_TO_RGRID") {
511 readEcho(in, inFileName);
512 readEcho(in, outFileName);
513 basisToRGrid(inFileName, outFileName);
515 if (command ==
"RGRID_TO_BASIS") {
516 readEcho(in, inFileName);
517 readEcho(in, outFileName);
518 rGridToBasis(inFileName, outFileName);
520 if (command ==
"KGRID_TO_RGRID") {
521 readEcho(in, inFileName);
522 readEcho(in, outFileName);
523 kGridToRGrid(inFileName, outFileName);
525 if (command ==
"RGRID_TO_KGRID") {
526 readEcho(in, inFileName);
527 readEcho(in, outFileName);
528 rGridToKGrid(inFileName, outFileName);
530 if (command ==
"BASIS_TO_KGRID") {
531 readEcho(in, inFileName);
532 readEcho(in, outFileName);
533 basisToKGrid(inFileName, outFileName);
535 if (command ==
"KGRID_TO_BASIS") {
536 readEcho(in, inFileName);
537 readEcho(in, outFileName);
538 kGridToBasis(inFileName, outFileName);
540 if (command ==
"CHECK_RGRID_SYMMETRY") {
542 readEcho(in, inFileName);
543 readEcho(in, epsilon);
545 hasSymmetry = checkRGridFieldSymmetry(inFileName, epsilon);
548 <<
"Symmetry of r-grid file matches this space group."
549 << std::endl << std::endl;
552 <<
"Symmetry of r-grid file does not match this\n"
553 <<
"space group to within error threshold of "
554 <<
Dbl(epsilon) <<
" ." << std::endl << std::endl;
557 if (command ==
"COMPARE_BASIS") {
560 std::string filecompare1, filecompare2;
561 readEcho(in, filecompare1);
562 readEcho(in, filecompare2);
566 domain().fieldIo().readFieldsBasis(filecompare1, Bfield1,
568 domain().fieldIo().readFieldsBasis(filecompare2, Bfield2,
573 compare(Bfield1, Bfield2);
576 if (command ==
"COMPARE_RGRID") {
578 std::string filecompare1, filecompare2;
579 readEcho(in, filecompare1);
580 readEcho(in, filecompare2);
584 domain().fieldIo().readFieldsRGrid(filecompare1, Rfield1,
586 domain().fieldIo().readFieldsRGrid(filecompare2, Rfield2,
591 compare(Rfield1, Rfield2);
594 if (command ==
"SCALE_BASIS") {
596 readEcho(in, inFileName);
597 readEcho(in, outFileName);
598 readEcho(in, factor);
599 scaleFieldsBasis(inFileName, outFileName, factor);
601 if (command ==
"SCALE_RGRID") {
603 readEcho(in, inFileName);
604 readEcho(in, outFileName);
605 readEcho(in, factor);
606 scaleFieldsRGrid(inFileName, outFileName, factor);
608 if (command ==
"EXPAND_RGRID_DIMENSION") {
609 readEcho(in, inFileName);
610 readEcho(in, outFileName);
616 Log::file() <<
Str(
"Expand dimension to: ") << d <<
"\n";
621 for (
int i = 0; i < d-D; i++) {
622 in >> newGridDimensions[i];
625 <<
Str(
"Numbers of grid points in added dimensions : ");
626 for (
int i = 0; i < d-D; i++) {
627 Log::file() << newGridDimensions[i] <<
" ";
632 d, newGridDimensions);
635 if (command ==
"REPLICATE_UNIT_CELL") {
636 readEcho(in, inFileName);
637 readEcho(in, outFileName);
641 for (
int i = 0; i < D; i++){
644 for (
int i = 0; i < D; i++){
645 Log::file() <<
"Replicate unit cell in direction "
653 if (command ==
"READ_H_BASIS") {
654 readEcho(in, filename);
655 if (!h_.isAllocatedBasis()) {
656 h_.allocateBasis(domain().basis().nBasis());
658 if (!h_.isAllocatedRGrid()) {
659 h_.allocateRGrid(domain().mesh().dimensions());
662 h_.readBasis(filename, tmpUnitCell);
664 if (command ==
"READ_H_RGRID") {
665 readEcho(in, filename);
666 if (!h_.isAllocatedRGrid()) {
667 h_.allocateRGrid(domain().mesh().dimensions());
669 if (iterator().isSymmetric() && !h_.isAllocatedBasis()) {
670 mask_.allocateBasis(domain().basis().nBasis());
673 h_.readRGrid(filename, tmpUnitCell);
675 if (command ==
"WRITE_H_BASIS") {
676 readEcho(in, filename);
679 domain().fieldIo().writeFieldsBasis(filename, h_.basis(),
680 domain().unitCell());
682 if (command ==
"WRITE_H_RGRID") {
683 readEcho(in, filename);
685 domain().fieldIo().writeFieldsRGrid(filename, h_.rgrid(),
686 domain().unitCell());
688 if (command ==
"READ_MASK_BASIS") {
689 readEcho(in, filename);
691 if (!mask_.isAllocatedBasis()) {
692 mask_.allocateBasis(domain().basis().nBasis());
694 if (!mask_.isAllocatedRGrid()) {
695 mask_.allocateRGrid(domain().mesh().dimensions());
698 mask_.readBasis(filename, tmpUnitCell);
700 if (command ==
"READ_MASK_RGRID") {
701 readEcho(in, filename);
702 if (!mask_.isAllocatedRGrid()) {
703 mask_.allocateRGrid(domain().mesh().dimensions());
705 if (iterator().isSymmetric() && !mask_.isAllocatedBasis()) {
707 mask_.allocateBasis(domain().basis().nBasis());
710 mask_.readRGrid(filename, tmpUnitCell);
712 if (command ==
"WRITE_MASK_BASIS") {
713 readEcho(in, filename);
716 domain().fieldIo().writeFieldBasis(filename, mask_.basis(),
717 domain().unitCell());
719 if (command ==
"WRITE_MASK_RGRID") {
720 readEcho(in, filename);
722 domain().fieldIo().writeFieldRGrid(filename, mask_.rgrid(),
724 mask_.isSymmetric());
727 << command << std::endl;
739 if (fileMaster_.commandFileName().empty()) {
742 readCommands(fileMaster_.commandFile());
756 if (!domain().unitCell().isInitialized()) {
759 UTIL_CHECK(domain().unitCell().isInitialized());
762 if (!isAllocatedBasis_) {
763 allocateFieldsBasis();
767 w_.readBasis(filename, domain_.unitCell());
770 mixture_.clearUnitCellData();
773 hasFreeEnergy_ =
false;
786 if (!domain().unitCell().isInitialized()) {
789 UTIL_CHECK(domain().unitCell().isInitialized());
790 if (domain().hasGroup() && !isAllocatedBasis_) {
793 allocateFieldsBasis();
797 w_.readRGrid(filename, domain_.unitCell());
800 mixture_.clearUnitCellData();
803 hasFreeEnergy_ =
false;
815 const int nm = mixture().nMonomer();
820 if (!domain().unitCell().isInitialized()) {
823 UTIL_CHECK(domain().unitCell().isInitialized());
825 const int nb = domain().basis().nBasis();
827 if (!isAllocatedBasis_) {
828 allocateFieldsBasis();
832 domain().fieldIo().readFieldsBasis(filename, tmpFieldsBasis_,
836 mixture_.clearUnitCellData();
844 for (i = 0; i < nb; ++i) {
845 for (j = 0; j < nm; ++j) {
847 for (k = 0; k < nm; ++k) {
848 wtmp[j] += interaction().chi(j,k)*tmpFieldsBasis_[k][i];
851 for (j = 0; j < nm; ++j) {
852 tmpFieldsBasis_[j][i] = wtmp[j];
857 w_.setBasis(tmpFieldsBasis_);
860 hasFreeEnergy_ =
false;
875 hasFreeEnergy_ =
false;
887 hasFreeEnergy_ =
false;
898 domain_.setUnitCell(unitCell);
900 mixture_.clearUnitCellData();
901 if (domain().hasGroup() && !isAllocatedBasis_) {
902 allocateFieldsBasis();
914 domain_.setUnitCell(lattice, parameters);
916 mixture_.clearUnitCellData();
917 if (domain().hasGroup() && !isAllocatedBasis_) {
918 allocateFieldsBasis();
928 domain_.setUnitCell(parameters);
930 mixture_.clearUnitCellData();
931 if (domain().hasGroup() && !isAllocatedBasis_) {
932 allocateFieldsBasis();
949 mixture_.compute(w_.rgrid(), c_.rgrid(), mask_.phiTot());
951 hasFreeEnergy_ =
false;
954 if (w_.isSymmetric()) {
956 domain().fieldIo().convertRGridToBasis(c_.rgrid(), c_.basis(),
962 mixture_.computeStress(mask().phiTot());
975 if (iterator().isSymmetric()) {
979 hasFreeEnergy_ =
false;
985 int error = iterator().solve(isContinuation);
991 if (!iterator().isFlexible()) {
992 mixture_.computeStress(mask().phiTot());
995 if (!iterator().isFlexible()) {
1026 hasCFields_ =
false;
1027 hasFreeEnergy_ =
false;
1029 simulator().simulate(nStep);
1038 std::string
const & outFileName,
1046 domain().fieldIo().readFieldsRGrid(inFileName,
1051 domain().fieldIo().expandRGridDimension(outFileName,
1054 d, newGridDimensions);
1062 std::string
const & outFileName,
1067 domain().fieldIo().readFieldsRGrid(inFileName, tmpFieldsRGrid_,
1071 domain().fieldIo().replicateUnitCell(outFileName,
1085 if (hasFreeEnergy_)
return;
1090 int nm = mixture().nMonomer();
1091 int np = mixture().nPolymer();
1092 int ns = mixture().nSolvent();
1106 for (
int i = 0; i < np; ++i) {
1107 polymerPtr = &mixture_.polymer(i);
1108 phi = polymerPtr->
phi();
1109 mu = polymerPtr->
mu();
1110 length = polymerPtr->
length();
1112 if (phi > 1.0E-08) {
1113 fIdeal_ += phi*( mu - 1.0 )/length;
1122 for (
int i = 0; i < ns; ++i) {
1123 solventPtr = &mixture_.solvent(i);
1124 phi = solventPtr->
phi();
1125 mu = solventPtr->
mu();
1126 size = solventPtr->
size();
1127 if (phi > 1.0E-08) {
1128 fIdeal_ += phi*( mu - 1.0 )/size;
1143 const double phiTot = mask().phiTot();
1147 if (w_.isSymmetric()) {
1149 const int nBasis = domain_.basis().nBasis();
1150 for (
int i = 0; i < nm; ++i) {
1151 for (
int k = 0; k < nBasis; ++k) {
1152 temp -= w_.basis(i)[k] * c_.basis(i)[k];
1157 const int meshSize = domain().mesh().size();
1158 for (
int i = 0; i < nm; ++i) {
1159 for (
int k = 0; k < meshSize; ++k) {
1160 temp -= w_.rgrid(i)[k] * c_.rgrid(i)[k];
1163 temp /= double(meshSize);
1167 fHelmholtz_ += fIdeal_;
1170 if (hasExternalFields()) {
1171 if (w_.isSymmetric()) {
1173 const int nBasis = domain_.basis().nBasis();
1174 for (
int i = 0; i < nm; ++i) {
1175 for (
int k = 0; k < nBasis; ++k) {
1176 fExt_ += h_.basis(i)[k] * c_.basis(i)[k];
1181 const int meshSize = domain().mesh().size();
1182 for (
int i = 0; i < nm; ++i) {
1183 for (
int k = 0; k < meshSize; ++k) {
1184 fExt_ += h_.rgrid(i)[k] * c_.rgrid(i)[k];
1187 fExt_ /= double(meshSize);
1190 fHelmholtz_ += fExt_;
1194 if (w_.isSymmetric()) {
1195 const int nBasis = domain_.basis().nBasis();
1196 for (
int i = 0; i < nm; ++i) {
1197 for (
int j = i; j < nm; ++j) {
1198 const double chi = interaction().chi(i,j);
1199 if (std::abs(chi) > 1.0E-9) {
1201 for (
int k = 0; k < nBasis; ++k) {
1202 temp += c_.basis(i)[k] * c_.basis(j)[k];
1205 fInter_ += 0.5*chi*temp;
1207 fInter_ += chi*temp;
1213 const int meshSize = domain().mesh().size();
1214 for (
int i = 0; i < nm; ++i) {
1215 for (
int j = i; j < nm; ++j) {
1216 const double chi = interaction().chi(i,j);
1217 if (std::abs(chi) > 1.0E-9) {
1219 for (
int k = 0; k < meshSize; ++k) {
1220 temp += c_.rgrid(i)[k] * c_.rgrid(j)[k];
1223 fInter_ += 0.5*chi*temp;
1225 fInter_ += chi*temp;
1230 fInter_ /= double(meshSize);
1233 fHelmholtz_ += fInter_;
1236 pressure_ = -fHelmholtz_;
1242 for (
int i = 0; i < np; ++i) {
1243 polymerPtr = &mixture_.polymer(i);
1244 phi = polymerPtr->
phi();
1245 mu = polymerPtr->
mu();
1246 length = polymerPtr->
length();
1247 if (phi > 1.0E-08) {
1248 pressure_ += mu * phi / length;
1257 for (
int i = 0; i < ns; ++i) {
1258 solventPtr = &mixture_.solvent(i);
1259 phi = solventPtr->
phi();
1260 mu = solventPtr->
mu();
1261 size = solventPtr->
size();
1262 if (phi > 1.0E-08) {
1263 pressure_ += mu * phi / size;
1268 hasFreeEnergy_ =
true;
1281 iterator().outputTimers(out);
1283 if (hasSimulator()){
1285 simulator().outputTimers(out);
1296 iterator().clearTimers();
1298 if (hasSimulator()){
1299 simulator().clearTimers();
1309 out <<
"System{" << std::endl;
1310 mixture_.writeParam(out);
1311 interaction().writeParam(out);
1312 domain_.writeParam(out);
1314 iterator().writeParam(out);
1316 out <<
"}" << std::endl;
1325 if (!hasFreeEnergy_) {
1326 computeFreeEnergy();
1330 out <<
"fHelmholtz " <<
Dbl(fHelmholtz(), 18, 11) << std::endl;
1331 out <<
"pressure " <<
Dbl(pressure(), 18, 11) << std::endl;
1333 out <<
"fIdeal " <<
Dbl(fIdeal_, 18, 11) << std::endl;
1334 out <<
"fInter " <<
Dbl(fInter_, 18, 11) << std::endl;
1335 if (hasExternalFields()) {
1336 out <<
"fExt " <<
Dbl(fExt_, 18, 11) << std::endl;
1340 int np = mixture().nPolymer();
1341 int ns = mixture().nSolvent();
1344 out <<
"polymers:" << std::endl;
1349 for (
int i = 0; i < np; ++i) {
1351 <<
" " <<
Dbl(mixture_.polymer(i).phi(),18, 11)
1352 <<
" " <<
Dbl(mixture_.polymer(i).mu(), 18, 11)
1359 out <<
"solvents:" << std::endl;
1364 for (
int i = 0; i < ns; ++i) {
1366 <<
" " <<
Dbl(mixture_.solvent(i).phi(),18, 11)
1367 <<
" " <<
Dbl(mixture_.solvent(i).mu(), 18, 11)
1373 out <<
"cellParams:" << std::endl;
1374 for (
int i = 0; i < domain().unitCell().nParameter(); ++i) {
1377 <<
Dbl(domain().unitCell().parameter(i), 18, 11)
1389 out <<
"stress:" << std::endl;
1390 for (
int i = 0; i < domain().unitCell().nParameter(); ++i) {
1393 <<
Dbl(mixture_.stress(i), 18, 11)
1409 domain_.fieldIo().writeFieldsBasis(filename, w_.basis(),
1410 domain().unitCell());
1421 domain_.fieldIo().writeFieldsRGrid(filename, w_.rgrid(),
1422 domain().unitCell(),
1436 domain_.fieldIo().writeFieldsBasis(filename, c_.basis(),
1437 domain().unitCell());
1448 domain_.fieldIo().writeFieldsRGrid(filename, c_.rgrid(),
1449 domain().unitCell(),
1468 mixture_.createBlockCRGrid(blockCFields);
1471 domain().fieldIo().writeFieldsRGrid(filename, blockCFields,
1472 domain().unitCell(),
1481 int polymerId,
int blockId,
1482 int directionId,
int segmentId)
1486 UTIL_CHECK(polymerId < mixture().nPolymer());
1487 Polymer<D> const& polymer = mixture_.polymer(polymerId);
1493 propagator = polymer.
propagator(blockId, directionId);
1494 RField<D> const& field = propagator.
q(segmentId);
1495 domain().fieldIo().writeFieldRGrid(filename, field,
1496 domain().unitCell(),
1505 int polymerId,
int blockId,
int directionId)
1509 UTIL_CHECK(polymerId < mixture().nPolymer());
1510 Polymer<D> const& polymer = mixture_.polymer(polymerId);
1517 domain().fieldIo().writeFieldRGrid(filename, field,
1518 domain().unitCell(),
1527 int polymerId,
int blockId,
int directionId)
1531 UTIL_CHECK(polymerId < mixture().nPolymer());
1532 Polymer<D> const& polymer = mixture_.polymer(polymerId);
1539 int ns = propagator.
ns();
1543 fileMaster_.openOutputFile(filename, file);
1546 domain().fieldIo().writeFieldHeader(file, 1, domain().unitCell(),
1548 file <<
"mesh " << std::endl
1549 <<
" " << domain().mesh().dimensions() << std::endl
1550 <<
"nslice" << std::endl
1551 <<
" " << ns << std::endl;
1554 bool hasHeader =
false;
1555 for (
int i = 0; i < ns; ++i) {
1556 file <<
"slice " << i << std::endl;
1557 domain().fieldIo().writeFieldRGrid(file, propagator.
q(i),
1558 domain().unitCell(),
1570 std::string filename;
1571 int np, nb, ip, ib, id;
1572 np = mixture().nPolymer();
1573 for (ip = 0; ip < np; ++ip) {
1574 nb = mixture_.polymer(ip).nBlock();
1575 for (ib = 0; ib < nb; ++ib) {
1576 for (
id = 0;
id < 2; ++id) {
1577 filename = basename;
1585 writeQ(filename, ip, ib,
id);
1599 fileMaster_.openOutputFile(filename, file);
1600 bool isSymmetric =
true;
1601 domain().fieldIo().writeFieldHeader(file, mixture().nMonomer(),
1602 domain().unitCell(),
1604 domain_.basis().outputStars(file);
1617 fileMaster_.openOutputFile(filename, file);
1618 bool isSymmetric =
true;
1619 domain().fieldIo().writeFieldHeader(file, mixture().nMonomer(),
1620 domain().unitCell(),
1622 domain_.basis().outputWaves(file);
1633 Pscf::Prdc::writeGroup(filename, domain_.group());
1643 std::string
const & outFileName)
1649 if (!isAllocatedBasis_) {
1651 allocateFieldsBasis();
1656 FieldIo<D> const & fieldIo = domain().fieldIo();
1668 std::string
const & outFileName)
1674 if (!isAllocatedBasis_) {
1676 allocateFieldsBasis();
1681 FieldIo<D> const & fieldIo = domain().fieldIo();
1693 std::string
const & outFileName)
1697 if (domain_.hasGroup() && !isAllocatedBasis_) {
1699 allocateFieldsBasis();
1704 FieldIo<D> const & fieldIo = domain().fieldIo();
1706 for (
int i = 0; i < mixture().nMonomer(); ++i) {
1707 domain().fft().inverseTransformUnsafe(tmpFieldsKGrid_[i],
1708 tmpFieldsRGrid_[i]);
1719 std::string
const & outFileName)
1723 if (domain_.hasGroup() && !isAllocatedBasis_) {
1725 allocateFieldsBasis();
1730 FieldIo<D> const & fieldIo = domain().fieldIo();
1733 for (
int i = 0; i < mixture().nMonomer(); ++i) {
1734 domain().fft().forwardTransform(tmpFieldsRGrid_[i],
1735 tmpFieldsKGrid_[i]);
1746 std::string
const & outFileName)
1752 if (!isAllocatedBasis_) {
1754 allocateFieldsBasis();
1759 domain_.fieldIo().readFieldsKGrid(inFileName, tmpFieldsKGrid_,
1761 domain_.fieldIo().convertKGridToBasis(tmpFieldsKGrid_,
1763 domain_.fieldIo().writeFieldsBasis(outFileName,
1764 tmpFieldsBasis_, tmpUnitCell);
1772 std::string
const & outFileName)
1778 if (!isAllocatedBasis_) {
1780 allocateFieldsBasis();
1785 domain_.fieldIo().readFieldsBasis(inFileName,
1786 tmpFieldsBasis_, tmpUnitCell);
1787 domain_.fieldIo().convertBasisToKGrid(tmpFieldsBasis_,
1789 domain_.fieldIo().writeFieldsKGrid(outFileName,
1790 tmpFieldsKGrid_, tmpUnitCell);
1805 if (!isAllocatedBasis_) {
1807 allocateFieldsBasis();
1812 domain_.fieldIo().readFieldsRGrid(inFileName,
1813 tmpFieldsRGrid_, tmpUnitCell);
1816 for (
int i = 0; i < mixture().nMonomer(); ++i) {
1818 symmetric = domain_.fieldIo().hasSymmetry(tmpFieldsRGrid_[i],
1833 std::string
const & outFileName,
1838 if (!isAllocatedBasis_) {
1840 allocateFieldsBasis();
1844 FieldIo<D> const & fieldIo = domain().fieldIo();
1856 std::string
const & outFileName,
1857 double factor)
const
1860 FieldIo<D> const & fieldIo = domain().fieldIo();
1878 comparison.
compare(field1,field2);
1880 Log::file() <<
"\n Basis expansion field comparison results"
1882 Log::file() <<
" Maximum Absolute Difference: "
1883 << comparison.
maxDiff() << std::endl;
1884 Log::file() <<
" Root-Mean-Square Difference: "
1885 << comparison.
rmsDiff() <<
"\n" << std::endl;
1896 comparison.
compare(field1, field2);
1898 Log::file() <<
"\n Real-space field comparison results"
1900 Log::file() <<
" Maximum Absolute Difference: "
1901 << comparison.
maxDiff() << std::endl;
1902 Log::file() <<
" Root-Mean-Square Difference: "
1903 << comparison.
rmsDiff() <<
"\n" << std::endl;
1916 int nMonomer = mixture().nMonomer();
1922 IntVec<D> const & dimensions = domain().mesh().dimensions();
1925 w_.setNMonomer(nMonomer);
1926 w_.allocateRGrid(dimensions);
1929 c_.setNMonomer(nMonomer);
1930 c_.allocateRGrid(dimensions);
1932 h_.setNMonomer(nMonomer);
1935 tmpFieldsRGrid_.allocate(nMonomer);
1936 tmpFieldsKGrid_.allocate(nMonomer);
1937 for (
int i = 0; i < nMonomer; ++i) {
1938 tmpFieldsRGrid_[i].allocate(dimensions);
1939 tmpFieldsKGrid_[i].allocate(dimensions);
1942 isAllocatedGrid_ =
true;
1949 void System<D>::allocateFieldsBasis()
1953 const int nMonomer = mixture().nMonomer();
1959 const int nBasis = domain_.basis().nBasis();
1962 w_.allocateBasis(nBasis);
1963 c_.allocateBasis(nBasis);
1966 tmpFieldsBasis_.allocate(nMonomer);
1967 for (
int i = 0; i < nMonomer; ++i) {
1968 tmpFieldsBasis_[i].allocate(nBasis);
1971 isAllocatedBasis_ =
true;
1978 void System<D>::readFieldHeader(std::string filename)
1985 fileMaster_.openInputFile(filename, file);
1990 domain_.fieldIo().readFieldHeader(file, nMonomer,
1991 domain_.unitCell(), isSymmetric);
1999 UTIL_CHECK(mixture().nMonomer() == nMonomer);
2000 UTIL_CHECK(domain().unitCell().nParameter() > 0);
2002 UTIL_CHECK(domain().unitCell().isInitialized());
2003 if (domain_.hasGroup() && isSymmetric) {
2014 void System<D>::readEcho(std::istream& in, std::string&
string)
const
2018 UTIL_THROW(
"Unable to read string parameter.");
2027 void System<D>::readEcho(std::istream& in,
double& value)
const
2031 UTIL_THROW(
"Unable to read floating point parameter.");
2041 void System<D>::initHomogeneous()
2045 int nm = mixture().nMonomer();
2046 int np = mixture().nPolymer();
2047 int ns = mixture().nSolvent();
2048 UTIL_CHECK(homogeneous_.nMolecule() == np + ns);
2066 for (i = 0; i < np; ++i) {
2069 for (j = 0; j < nm; ++j) {
2074 nb = mixture_.polymer(i).nBlock();
2075 for (k = 0; k < nb; ++k) {
2076 Block<D>& block = mixture_.polymer(i).block(k);
2077 j = block.monomerId();
2078 s[j] += block.length();
2083 for (j = 0; j < nm; ++j) {
2084 if (s[j] > 1.0E-8) {
2088 homogeneous_.molecule(i).setNClump(nc);
2092 for (j = 0; j < nm; ++j) {
2093 if (s[j] > 1.0E-8) {
2094 homogeneous_.molecule(i).clump(k).setMonomerId(j);
2095 homogeneous_.molecule(i).clump(k).setSize(s[j]);
2099 homogeneous_.molecule(i).computeSize();
2108 for (
int is = 0; is < ns; ++is) {
2110 monomerId = mixture_.solvent(is).monomerId();
2111 size = mixture_.solvent(is).size();
2112 homogeneous_.molecule(i).setNClump(1);
2113 homogeneous_.molecule(i).clump(0).setMonomerId(monomerId);
2114 homogeneous_.molecule(i).clump(0).setSize(size);
2115 homogeneous_.molecule(i).computeSize();
double rmsDiff() const
Return the precomputed root-mean-squared difference.
double compare(FT const &a, FT const &b)
Compare individual fields.
double maxDiff() const
Return the precomputed maximum element-by-element difference.
An IntVec<D, T> is a D-component vector of elements of integer type T.
Flory-Huggins excess free energy model.
int nBlock() const
Number of blocks.
Propagator & propagator(int blockId, int directionId)
Get propagator for a specific block and direction.
double length() const
Sum of the lengths of all blocks in the polymer.
Comparator for fields in symmetry-adapted basis format.
Comparator for fields in real-space (r-grid) format.
Field of real double precision values on an FFT mesh.
void scaleFieldsRGrid(DArray< RFRT > &fields, double factor) const
Scale an array of r-grid fields by a real scalar.
void convertBasisToRGrid(DArray< double > const &in, RFRT &out) const
Convert a field from symmetrized basis to spatial grid (r-grid).
void convertRGridToBasis(RFRT const &in, DArray< double > &out, bool checkSymmetry=true, double epsilon=1.0e-8) const
Convert a field from spatial grid (r-grid) to symmetrized basis.
void readFieldsBasis(std::istream &in, DArray< DArray< double > > &fields, UnitCell< D > &unitCell) const
Read concentration or chemical potential fields from file.
void scaleFieldsBasis(DArray< DArray< double > > &fields, double factor) const
Scale an array of fields in basis format by a real scalar.
void writeFieldsBasis(std::ostream &out, DArray< DArray< double > > const &fields, UnitCell< D > const &unitCell) const
Write concentration or chemical potential field components to file.
Base template for UnitCell<D> classes, D=1, 2 or 3.
File input/output operations and format conversions for fields.
void writeFieldsKGrid(std::ostream &out, DArray< RFieldDft< D > > const &fields, UnitCell< D > const &unitCell, bool isSymmetric=true) const override
Write array of RFieldDft objects (k-space fields) to file.
void readFieldsRGrid(std::istream &in, DArray< RField< D > > &fields, UnitCell< D > &unitCell) const override
Read array of RField objects (r-grid fields) from a stream.
void readFieldsKGrid(std::istream &in, DArray< RFieldDft< D > > &fields, UnitCell< D > &unitCell) const override
Read array of RFieldDft objects (k-space fields) from a stream.
void writeFieldsRGrid(std::ostream &out, DArray< RField< D > > const &fields, UnitCell< D > const &unitCell, bool writeHeader=true, bool isSymmetric=true, bool writeMeshSize=true) const override
Write array of RField objects (fields on r-space grid) to a stream.
Factory for subclasses of Iterator.
Descriptor and solver for one polymer species.
MDE solver for one direction of one block.
const QField & q(int i) const
Return q-field at specified step.
int ns() const
Number of values of s (or slices), including head and tail.
const QField & tail() const
Return q-field at the end of the block.
Factory for subclasses of Simulator.
Solver and descriptor for a solvent species.
Default Factory for subclasses of Sweep.
Main class for SCFT or PS-FTS simulation of one system.
void expandRGridDimension(const std::string &inFileName, const std::string &outFileName, int d, DArray< int > newGridDimensions)
Expand the number of spatial dimensions of an r-grid field.
int iterate(bool isContinuation=false)
Iteratively solve a SCFT problem.
void writeGroup(std::string const &filename) const
Output all elements of the space group.
void scaleFieldsBasis(const std::string &inFileName, const std::string &outFileName, double factor)
Multiply all components of an array of basis fields by a scalar.
void writeQTail(std::string const &filename, int polymerId, int blockId, int directionId) const
Write the final slice of a propagator in r-grid format.
void readWBasis(const std::string &filename)
Read chemical potential fields in symmetry adapted basis format.
void writeWBasis(const std::string &filename) const
Write chemical potential fields in symmetrized basis format.
void compute(bool needStress=false)
Solve the modified diffusion equation once, without iteration.
void setWRGrid(DArray< RField< D > > const &fields)
Set new w fields, in real-space (r-grid) format.
void writeWRGrid(const std::string &filename) const
Write chemical potential fields in real space grid (r-grid) format.
void setOptions(int argc, char **argv)
Process command line options.
void readCommands()
Read and process commands from the default command file.
void simulate(int nStep)
Perform a field theoretic simulation (PS-FTS).
void writeQAll(std::string const &basename)
Write all propagators of all blocks, each to a separate file.
bool checkRGridFieldSymmetry(const std::string &inFileName, double epsilon=1.0E-8)
Check if r-grid fields have the declared space group symmetry.
void replicateUnitCell(const std::string &inFileName, const std::string &outFileName, IntVec< D > const &replicas)
Replicate the crystal unit cell to create a larger cell.
void writeTimers(std::ostream &out)
Write timer information to an output stream.
void compare(const DArray< DArray< double > > field1, const DArray< DArray< double > > field2)
Compare arrays of fields in basis format, output a report.
void writeCBasis(const std::string &filename) const
Write concentration fields in symmetrized basis format.
void writeThermo(std::ostream &out)
Write thermodynamic properties to a file.
void setUnitCell(UnitCell< D > const &unitCell)
Set parameters of the associated unit cell.
void rGridToKGrid(const std::string &inFileName, const std::string &outFileName)
Convert fields from real-space (r-grid) to Fourier (k-grid) format.
void basisToKGrid(const std::string &inFileName, const std::string &outFileName)
Convert fields from symmetrized basis to Fourier (k-grid) format.
void computeFreeEnergy()
Compute free energy density and pressure for current fields.
void setWBasis(DArray< DArray< double > > const &fields)
Set chemical potential fields, in symmetry-adapted basis format.
void clearTimers()
Clear timers.
void scaleFieldsRGrid(const std::string &inFileName, const std::string &outFileName, double factor) const
Multiply all elements of an array of r-grid fields by a scalar.
void writeQ(std::string const &filename, int polymerId, int blockId, int directionId) const
Write one propagator for one block, in r-grid format.
void writeBlockCRGrid(const std::string &filename) const
Write c-fields for all blocks and solvents in r-grid format.
virtual void readParameters(std::istream &in)
Read body of parameter block (without opening and closing lines).
void writeParamNoSweep(std::ostream &out) const
Write parameter file to an ostream, omitting any sweep block.
void sweep()
Sweep in parameter space, solving an SCF problem at each point.
void kGridToBasis(const std::string &inFileName, const std::string &outFileName)
Convert fields from Fourier (k-grid) to symmetrized basis format.
void readWRGrid(const std::string &filename)
Read chemical potential fields in real space grid (r-grid) format.
void kGridToRGrid(const std::string &inFileName, const std::string &outFileName)
Convert fields from Fourier (k-grid) to real-space (r-grid) format.
void estimateWfromC(const std::string &filename)
Construct trial w-fields from c-fields in symmetry-adapted form.
void writeCRGrid(const std::string &filename) const
Write concentration fields in real space grid (r-grid) format.
Domain< D > const & domain() const
Get Domain by const reference.
void writeQSlice(std::string const &filename, int polymerId, int blockId, int directionId, int segmentId) const
Write slice of a propagator at fixed s in r-grid format.
void basisToRGrid(const std::string &inFileName, const std::string &outFileName)
Convert a field from symmetrized basis format to r-grid format.
void writeStars(std::string const &filename) const
Output information about stars and symmetrized basis functions.
void readParam()
Read input parameters from default param file.
void writeStress(std::ostream &out)
Write stress properties to a file.
void rGridToBasis(const std::string &inFileName, const std::string &outFileName)
Convert a field from real-space grid to symmetrized basis format.
void writeWaves(std::string const &filename) const
Output information about waves.
double size() const
Get the size (number of monomers) in this solvent.
double phi() const
Get the overall volume fraction for this species.
double mu() const
Get the chemical potential for this species (units kT=1).
Dynamically allocatable contiguous array template.
void allocate(int capacity)
Allocate the underlying C array.
Wrapper for a double precision number, for formatted ostream output.
A fixed capacity (static) contiguous array with a variable logical size.
Wrapper for an int, for formatted ostream output.
static std::ostream & file()
Get log ostream by reference.
static bool echo()
Get echo parameter.
static void setEcho(bool echo=true)
Enable or disable echoing for all subclasses of ParamComponent.
void setClassName(const char *className)
Set class name string.
Wrapper for a std::string, for formatted ostream output.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
std::string toString(int n)
Return string representation of an integer.
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 replicateUnitCell(std::ostream &out, DArray< AT > const &fields, IntVec< D > const &meshDimensions, UnitCell< D > const &unitCell, IntVec< D > const &replicas)
Write r-grid fields in a replicated unit cell to std::ostream.
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 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.
Fields and FFTs for periodic boundary conditions (CPU)
Periodic fields and crystallography.
PSCF package top-level namespace.
void set(BracketPolicy::Type policy)
Set policy regarding use of bracket delimiters on arrays.
Utility classes for scientific computation.