PSCF v1.3
r1d/System.cpp
1/*
2* PSCF - Polymer Self-Consistent Field
3*
4* Copyright 2015 - 2025, The Regents of the University of Minnesota
5* Distributed under the terms of the GNU General Public License.
6*/
7
8#include "System.h"
9
10#include <r1d/iterator/Iterator.h>
11#include <r1d/iterator/IteratorFactory.h>
12#include <r1d/sweep/Sweep.h>
13#include <r1d/sweep/SweepFactory.h>
14#include <r1d/iterator/NrIterator.h>
15#include <r1d/misc/HomogeneousComparison.h>
16#include <r1d/misc/FieldIo.h>
17
18#include <pscf/inter/Interaction.h>
19#include <pscf/inter/Interaction.h>
20#include <pscf/floryHuggins/Clump.h>
21
22#include <util/format/Str.h>
23#include <util/format/Int.h>
24#include <util/format/Dbl.h>
25#include <util/param/BracketPolicy.h>
26#include <util/misc/ioUtil.h>
27
28#include <string>
29#include <unistd.h>
30
31namespace Pscf {
32namespace R1d
33{
34
35 using namespace Util;
36
37 /*
38 * Constructor.
39 */
41 : mixture_(),
42 domain_(),
43 fileMaster_(),
44 fieldIo_(),
45 homogeneous_(),
46 interactionPtr_(0),
47 iteratorPtr_(0),
48 iteratorFactoryPtr_(0),
49 sweepPtr_(0),
50 sweepFactoryPtr_(0),
51 wFields_(),
52 cFields_(),
53 f_(),
54 c_(),
55 fHelmholtz_(0.0),
56 fIdeal_(0.0),
57 fInter_(0.0),
58 pressure_(0.0),
59 hasMixture_(0),
60 hasDomain_(0),
61 hasFields_(0),
62 hasSweep_(0)
63 {
64 setClassName("System");
65
66 fieldIo_.associate(domain_, fileMaster_);
67 interactionPtr_ = new Interaction();
68 iteratorFactoryPtr_ = new IteratorFactory(*this);
69 sweepFactoryPtr_ = new SweepFactory(*this);
70
71 BracketPolicy::set(BracketPolicy::Optional);
72 }
73
74 /*
75 * Destructor.
76 */
78 {
79 if (interactionPtr_) {
80 delete interactionPtr_;
81 }
82 if (iteratorPtr_) {
83 delete iteratorPtr_;
84 }
85 if (iteratorFactoryPtr_) {
86 delete iteratorFactoryPtr_;
87 }
88 if (sweepPtr_) {
89 delete sweepPtr_;
90 }
91 if (sweepFactoryPtr_) {
92 delete sweepFactoryPtr_;
93 }
94 }
95
96 /*
97 * Process command line options.
98 */
99 void System::setOptions(int argc, char **argv)
100 {
101 bool eflag = false; // echo
102 bool pFlag = false; // param file
103 bool cFlag = false; // command file
104 bool iFlag = false; // input prefix
105 bool oFlag = false; // output prefix
106 char* pArg = 0;
107 char* cArg = 0;
108 char* iArg = 0;
109 char* oArg = 0;
110
111 // Read program arguments
112 int c;
113 opterr = 0;
114 while ((c = getopt(argc, argv, "er:p:c:i:o:f")) != -1) {
115 switch (c) {
116 case 'e':
117 eflag = true;
118 break;
119 case 'p': // parameter file
120 pFlag = true;
121 pArg = optarg;
122 break;
123 case 'c': // command file
124 cFlag = true;
125 cArg = optarg;
126 break;
127 case 'i': // input prefix
128 iFlag = true;
129 iArg = optarg;
130 break;
131 case 'o': // output prefix
132 oFlag = true;
133 oArg = optarg;
134 break;
135 case '?':
136 Log::file() << "Unknown option -" << optopt << std::endl;
137 UTIL_THROW("Invalid command line option");
138 }
139 }
140
141 // Set flag to echo parameters as they are read.
142 if (eflag) {
144 }
145
146 // If option -p, set parameter file name
147 if (pFlag) {
148 fileMaster().setParamFileName(std::string(pArg));
149 }
150
151 // If option -c, set command file name
152 if (cFlag) {
153 fileMaster().setCommandFileName(std::string(cArg));
154 }
155
156 // If option -i, set path prefix for input files
157 if (iFlag) {
158 fileMaster().setInputPrefix(std::string(iArg));
159 }
160
161 // If option -o, set path prefix for output files
162 if (oFlag) {
163 fileMaster().setOutputPrefix(std::string(oArg));
164 }
165
166 }
167
168 /*
169 * Read parameters and initialize.
170 */
171 void System::readParameters(std::istream& in)
172 {
174 hasMixture_ = true;
175
176 // Initialize FloryHuggins::Mixture object
177 homogeneous_.initialize(mixture());
178
179 interaction().setNMonomer(mixture().nMonomer());
181
183 hasDomain_ = true;
184 allocateFields();
185
186 // Instantiate and initialize an Iterator
187 std::string className;
188 bool isEnd;
189 iteratorPtr_ = iteratorFactoryPtr_->readObject(in, *this,
190 className, isEnd);
191 if (!iteratorPtr_) {
192 std::string msg = "Unrecognized Iterator subclass name ";
193 msg += className;
194 UTIL_THROW(msg.c_str());
195 }
196
197 // Optionally instantiate and initialize a Sweep object
198 sweepPtr_ =
199 sweepFactoryPtr_->readObjectOptional(in, *this, className, isEnd);
200 if (sweepPtr_) {
201 hasSweep_ = true;
202 sweepPtr_->setSystem(*this);
203 } else {
204 hasSweep_ = false;
205 }
206 }
207
208 /*
209 * Read default parameter file.
210 */
211 void System::readParam(std::istream& in)
212 {
213 readBegin(in, className().c_str());
214 readParameters(in);
215 readEnd(in);
216 }
217
218 /*
219 * Read default parameter file.
220 */
222 { readParam(fileMaster().paramFile()); }
223
224 /*
225 * Write parameter file, omitting the sweep block.
226 */
227 void System::writeParamNoSweep(std::ostream& out) const
228 {
229 out << "System{" << std::endl;
230 mixture_.writeParam(out);
231 interactionPtr_->writeParam(out);
232 domain_.writeParam(out);
233 iteratorPtr_->writeParam(out);
234 out << "}" << std::endl;
235 }
236
237 /*
238 * Allocate memory for fields.
239 */
240 void System::allocateFields()
241 {
242 // Preconditions
243 UTIL_CHECK(hasMixture_);
244 UTIL_CHECK(hasDomain_);
245
246 // Allocate memory in mixture
248
249 // Allocate wFields and cFields
250 int nMonomer = mixture().nMonomer();
251 wFields_.allocate(nMonomer);
252 cFields_.allocate(nMonomer);
253 int nx = domain().nx();
254 for (int i = 0; i < nMonomer; ++i) {
255 wField(i).allocate(nx);
256 cField(i).allocate(nx);
257 }
258 hasFields_ = true;
259 }
260
261 /*
262 * Read and execute commands from a specified command file.
263 */
264 void System::readCommands(std::istream &in)
265 {
266 UTIL_CHECK(hasMixture_);
267 UTIL_CHECK(hasDomain_);
268
269 std::string command;
270 std::string filename;
271
272 std::istream& inBuffer = in;
273
274 bool readNext = true;
275 while (readNext) {
276
277 inBuffer >> command;
278
279 if (inBuffer.eof()) {
280 break;
281 } else {
282 Log::file() << command;
283 }
284
285 if (command == "FINISH") {
286 Log::file() << std::endl;
287 readNext = false;
288 } else
289 if (command == "READ_W") {
290 readEcho(inBuffer, filename);
291 fieldIo_.readFields(wFields(), filename);
292 } else
293 if (command == "COMPUTE") {
294 // Solve the modified diffusion equation, without iteration
295 Log::file() << std::endl;
296 compute();
297 } else
298 if (command == "ITERATE") {
299 // Attempt iteratively solve a single SCFT problem
300 bool isContinuation = false;
301 int fail = iterate(isContinuation);
302 if (fail) {
303 readNext = false;
304 }
305 } else
306 if (command == "SWEEP") {
307 // Attempt to solve a sequence of SCFT problems along a line
308 // through parameter space.
309 sweep();
310 } else
311 if (command == "WRITE_PARAM") {
312 readEcho(inBuffer, filename);
313 std::ofstream file;
314 fileMaster().openOutputFile(filename, file);
315 writeParam(file);
316 file.close();
317 } else
318 if (command == "WRITE_THERMO") {
319 readEcho(inBuffer, filename);
320 std::ofstream file;
321 fileMaster().openOutputFile(filename, file, std::ios_base::app);
322 writeThermo(file);
323 file.close();
324 } else
325 if (command == "COMPARE_HOMOGENEOUS") {
326 int mode;
327 inBuffer >> mode;
328
329 readEcho(inBuffer, filename);
330 Log::file() << std::endl;
331 Log::file() << "mode = " << mode << std::endl;
332
333 std::ofstream file;
334 fileMaster().openOutputFile(filename, file, std::ios_base::app);
335
336 HomogeneousComparison comparison(*this);
337 comparison.compute(mode);
338 comparison.output(mode, file);
339 file.close();
340 } else
341 if (command == "WRITE_W") {
342 readEcho(inBuffer, filename);
343 fieldIo_.writeFields(wFields(), filename);
344 } else
345 if (command == "WRITE_C") {
346 readEcho(inBuffer, filename);
347 fieldIo_.writeFields(cFields(), filename);
348 } else
349 if (command == "WRITE_BLOCK_C") {
350 readEcho(inBuffer, filename);
351 fieldIo_.writeBlockCFields(mixture_, filename);
352 } else
353 if (command == "WRITE_Q_SLICE") {
354 int polymerId, blockId, directionId, segmentId;
355 readEcho(inBuffer, filename);
356 inBuffer >> polymerId;
357 inBuffer >> blockId;
358 inBuffer >> directionId;
359 inBuffer >> segmentId;
360 Log::file() << Str("polymer ID ", 21) << polymerId << "\n"
361 << Str("block ID ", 21) << blockId << "\n"
362 << Str("direction ID ", 21) << directionId << "\n"
363 << Str("segment ID ", 21) << segmentId << std::endl;
364 writeQSlice(filename, polymerId, blockId, directionId,
365 segmentId);
366 } else
367 if (command == "WRITE_Q_TAIL") {
368 readEcho(inBuffer, filename);
369 int polymerId, blockId, directionId;
370 inBuffer >> polymerId;
371 inBuffer >> blockId;
372 inBuffer >> directionId;
373 Log::file() << Str("polymer ID ", 21) << polymerId << "\n"
374 << Str("block ID ", 21) << blockId << "\n"
375 << Str("direction ID ", 21) << directionId << "\n";
376 writeQTail(filename, polymerId, blockId, directionId);
377 } else
378 if (command == "WRITE_Q_VERTEX") {
379 int polymerId, vertexId;
380 inBuffer >> polymerId;
381 Log::file() << std::endl;
382 Log::file() << "polymerId = "
383 << Int(polymerId, 5) << std::endl;
384 inBuffer >> vertexId;
385 Log::file() << "vertexId = "
386 << Int(vertexId, 5) << std::endl;
387 inBuffer >> filename;
388 Log::file() << "outfile = "
389 << Str(filename, 20) << std::endl;
390 fieldIo_.writeVertexQ(mixture_, polymerId, vertexId, filename);
391 } else
392 if (command == "WRITE_Q") {
393 readEcho(inBuffer, filename);
394 int polymerId, blockId, directionId;
395 inBuffer >> polymerId;
396 inBuffer >> blockId;
397 inBuffer >> directionId;
398 Log::file() << Str("polymer ID ", 21) << polymerId << "\n"
399 << Str("block ID ", 21) << blockId << "\n"
400 << Str("direction ID ", 21) << directionId << "\n";
401 writeQ(filename, polymerId, blockId, directionId);
402 } else
403 if (command == "WRITE_Q_ALL") {
404 readEcho(inBuffer, filename);
405 writeQAll(filename);
406 } else
407 if (command == "REMESH_W") {
408 int nx;
409 inBuffer >> nx;
410 Log::file() << std::endl;
411 Log::file() << "nx = " << Int(nx, 20) << std::endl;
412 inBuffer >> filename;
413 Log::file() << "outfile = " << Str(filename, 20) << std::endl;
414 fieldIo_.remesh(wFields(), nx, filename);
415 } else
416 if (command == "EXTEND_W") {
417 int m;
418 inBuffer >> m;
419 Log::file() << std::endl;
420 Log::file() << "m = " << Int(m, 20) << std::endl;
421 inBuffer >> filename;
422 Log::file() << "outfile = " << Str(filename, 20) << std::endl;
423 fieldIo_.extend(wFields(), m, filename);
424 } else {
425 Log::file() << " Error: Unknown command "
426 << command << std::endl;
427 readNext = false;
428 }
429
430 }
431 }
432
433 /*
434 * Read and execute commands from the default command file.
435 */
437 {
438 if (fileMaster().commandFileName().empty()) {
439 UTIL_THROW("Empty command file name");
440 }
441 readCommands(fileMaster().commandFile());
442 }
443
444 /*
445 * Read a string (e.g., a filename) and echo to the log file.
446 */
447 void System::readEcho(std::istream& in, std::string& string) const
448 {
449 in >> string;
450 Log::file() << " " << Str(string, 20) << std::endl;
451 }
452
453 // Primary SCFT Computations
454
455 /*
456 * Solve MDE for current w-fields, without iteration.
457 */
459 {
460 //UTIL_CHECK(hasWFields_);
461
462 // Solve the modified diffusion equation (without iteration)
463 mixture_.compute(wFields(), cFields_);
464
465 //hasCFields_ = true;
466 }
467
468 /*
469 * Iteratively solve a SCFT problem for specified parameters.
470 */
471 int System::iterate(bool isContinuation)
472 {
473 // UTIL_CHECK(hasWFields_);
474 // hasCFields_ = false;
475
476 Log::file() << std::endl;
477
478 // Call iterator (return 0 for convergence, 1 for failure)
479 int error = iterator().solve(isContinuation);
480
481 //hasCFields_ = true;
482
483 // If converged, compute related properties
484 if (!error) {
487 }
488 return error;
489 }
490
491 /*
492 * Perform sweep along a line in parameter space.
493 */
495 {
496 //UTIL_CHECK(hasWFields_);
497 UTIL_CHECK(sweepPtr_);
498 Log::file() << std::endl;
499 Log::file() << std::endl;
500
501 sweepPtr_->sweep();
502 }
503
504 // Thermodynamic properties
505
506 /*
507 * Compute Helmholtz free energy and pressure
508 */
510 {
511 fHelmholtz_ = 0.0;
512 int np = mixture().nPolymer();
513 int ns = mixture().nSolvent();
514 int nm = mixture().nMonomer();
515 double phi, mu;
516
517 // Polymeric ideal gas contributions to fHelhmoltz_
518 if (np > 0) {
519 Polymer* polymerPtr;
520 double length;
521 for (int i = 0; i < np; ++i) {
522 polymerPtr = &mixture().polymer(i);
523 phi = polymerPtr->phi();
524 mu = polymerPtr->mu();
525 // Recall: mu = ln(phi/q)
526 length = polymerPtr->length();
527 if (phi > 1E-08) {
528 fHelmholtz_ += phi*( mu - 1.0 )/length;
529 }
530 }
531 }
532
533 // Solvent ideal gas contributions to fHelhmoltz_
534 if (ns > 0) {
535 Solvent* solventPtr;
536 double size;
537 for (int i = 0; i < ns; ++i) {
538 solventPtr = &mixture().solvent(i);
539 phi = solventPtr->phi();
540 mu = solventPtr->mu();
541 // Recall: mu = ln(phi/q)
542 size = solventPtr->size();
543 if (phi > 1E-08) {
544 fHelmholtz_ += phi*( mu - 1.0 )/size;
545 }
546 }
547 }
548
549 // Apply Legendre transform subtraction
550 for (int i = 0; i < nm; ++i) {
551 fHelmholtz_ -=
552 domain().innerProduct(wFields_[i], cFields_[i]);
553 }
554 fIdeal_ = fHelmholtz_;
555
556 // Add average interaction free energy density per monomer
557 int nx = domain().nx();
558 if (!f_.isAllocated()) f_.allocate(nx);
559 if (!c_.isAllocated()) c_.allocate(nm);
560 int j;
561 for (int i = 0; i < nx; ++i) {
562 // Store c_[j] = local concentration of species j
563 for (j = 0; j < nm; ++j) {
564 c_[j] = cFields_[j][i];
565 }
566 // Compute f_[i] = excess free eenrgy at grid point i
567 f_[i] = interaction().fHelmholtz(c_);
568 }
569 fHelmholtz_ += domain().spatialAverage(f_);
570 fInter_ = fHelmholtz_ - fIdeal_;
571
572 // Compute pressure
573 pressure_ = -fHelmholtz_;
574 if (np > 0) {
575 double length;
576 Polymer* polymerPtr;
577 for (int i = 0; i < np; ++i) {
578 polymerPtr = & mixture().polymer(i);
579 phi = polymerPtr->phi();
580 mu = polymerPtr->mu();
581 length = polymerPtr->length();
582 if (phi > 1E-08) {
583 pressure_ += phi*mu/length;
584 }
585 }
586 }
587 if (ns > 0) {
588 double size;
589 Solvent* solventPtr;
590 for (int i = 0; i < ns; ++i) {
591 solventPtr = & mixture().solvent(i);
592 phi = solventPtr->phi();
593 mu = solventPtr->mu();
594 size = solventPtr->size();
595 if (phi > 1E-08) {
596 pressure_ += phi*mu/size;
597 }
598 }
599 }
600
601 }
602
603 void System::writeThermo(std::ostream& out) const
604 {
605 out << std::endl;
606 out << "fHelmholtz " << Dbl(fHelmholtz(), 18, 11) << std::endl;
607 out << "pressure " << Dbl(pressure(), 18, 11) << std::endl;
608 out << std::endl;
609 out << "fIdeal " << Dbl(fIdeal_, 18, 11) << std::endl;
610 out << "fInter " << Dbl(fInter_, 18, 11) << std::endl;
611 out << std::endl;
612
613 // Polymers
614 int np = mixture().nPolymer();
615 if (np > 0) {
616 out << "polymers:" << std::endl;
617 out << " "
618 << " phi "
619 << " mu "
620 << std::endl;
621 for (int i = 0; i < np; ++i) {
622 out << Int(i, 5)
623 << " " << Dbl(mixture().polymer(i).phi(),18, 11)
624 << " " << Dbl(mixture().polymer(i).mu(), 18, 11)
625 << std::endl;
626 }
627 out << std::endl;
628 }
629
630 // Solvents
631 int ns = mixture().nSolvent();
632 if (ns > 0) {
633 out << "solvents:" << std::endl;
634 out << " "
635 << " phi "
636 << " mu "
637 << std::endl;
638 for (int i = 0; i < ns; ++i) {
639 out << Int(i, 5)
640 << " " << Dbl(mixture().solvent(i).phi(),18, 11)
641 << " " << Dbl(mixture().solvent(i).mu(), 18, 11)
642 << std::endl;
643 }
644 out << std::endl;
645 }
646
647 }
648
649 // Output operations (correspond to command file commands)
650
651 /*
652 * Write w-fields in symmetry-adapted basis format.
653 */
654 void System::writeW(std::string const & filename) const
655 {
656 //UTIL_CHECK(hasWFields_);
657 fieldIo_.writeFields(wFields(), filename);
658 }
659
660 /*
661 * Write all concentration fields in symmetry-adapted basis format.
662 */
663 void System::writeC(std::string const & filename) const
664 {
665 //UTIL_CHECK(hasCFields_);
666 fieldIo_.writeFields(cFields(), filename);
667 }
668
669 /*
670 * Write all concentration fields in real space (r-grid) format, for
671 * each block (or solvent) individually rather than for each species.
672 */
673 void System::writeBlockC(std::string const & filename) const
674 {
675 //UTIL_CHECK(hasCFields_);
676 fieldIo_.writeBlockCFields(mixture_, filename);
677 }
678
679 /*
680 * Write the last time slice of the propagator in r-grid format.
681 */
682 void System::writeQSlice(const std::string & filename,
683 int polymerId, int blockId,
684 int directionId, int segmentId)
685 const
686 {
687 UTIL_CHECK(polymerId >= 0);
688 UTIL_CHECK(polymerId < mixture_.nPolymer());
689 Polymer const& polymer = mixture_.polymer(polymerId);
690 UTIL_CHECK(blockId >= 0);
691 UTIL_CHECK(blockId < polymer.nBlock());
692 UTIL_CHECK(directionId >= 0);
693 UTIL_CHECK(directionId <= 1);
694 Propagator const &
695 propagator = polymer.propagator(blockId, directionId);
696 Field const& field = propagator.q(segmentId);
697 fieldIo_.writeField(field, filename);
698 }
699
700 /*
701 * Write the last time slice of the propagator in r-grid format.
702 */
703 void System::writeQTail(const std::string & filename,
704 int polymerId, int blockId, int directionId)
705 const
706 {
707 UTIL_CHECK(polymerId >= 0);
708 UTIL_CHECK(polymerId < mixture_.nPolymer());
709 Polymer const& polymer = mixture_.polymer(polymerId);
710 UTIL_CHECK(blockId >= 0);
711 UTIL_CHECK(blockId < polymer.nBlock());
712 UTIL_CHECK(directionId >= 0);
713 UTIL_CHECK(directionId <= 1);
714 Field const&
715 field = polymer.propagator(blockId, directionId).tail();
716 fieldIo_.writeField(field, filename);
717 }
718
719 /*
720 * Write the propagator for a block and direction.
721 */
722 void System::writeQ(const std::string & filename,
723 int polymerId, int blockId, int directionId)
724 const
725 {
726 UTIL_CHECK(polymerId >= 0);
727 UTIL_CHECK(polymerId < mixture_.nPolymer());
728 Polymer const& polymer = mixture_.polymer(polymerId);
729 UTIL_CHECK(blockId >= 0);
730 UTIL_CHECK(blockId < polymer.nBlock());
731 UTIL_CHECK(directionId >= 0);
732 UTIL_CHECK(directionId <= 1);
733 Propagator const& propagator
734 = polymer.propagator(blockId, directionId);
735 int nslice = propagator.ns();
736
737 // Open file
738 std::ofstream file;
739 fileMaster_.openOutputFile(filename, file);
740
741 // Write header
742 file << "ngrid" << std::endl
743 << " " << domain_.nx() << std::endl
744 << "nslice" << std::endl
745 << " " << nslice << std::endl;
746
747 // Write data
748 bool hasHeader = false;
749 for (int i = 0; i < nslice; ++i) {
750 file << "slice " << i << std::endl;
751 Field const& field = propagator.q(i);
752 fieldIo_.writeField(field, file, hasHeader);
753 }
754 file.close();
755 }
756
757 /*
758 * Write propagators for all blocks of all polymers to files.
759 */
760 void System::writeQAll(std::string const & basename) const
761 {
762 std::string filename;
763 int np, nb, ip, ib, id;
764 np = mixture_.nPolymer();
765 for (ip = 0; ip < np; ++ip) {
766 nb = mixture_.polymer(ip).nBlock();
767 for (ib = 0; ib < nb; ++ib) {
768 for (id = 0; id < 2; ++id) {
769 filename = basename;
770 filename += "_";
771 filename += toString(ip);
772 filename += "_";
773 filename += toString(ib);
774 filename += "_";
775 filename += toString(id);
776 //filename += ".q";
777 writeQ(filename, ip, ib, id);
778 }
779 }
780 }
781 }
782
783} // namespace R1d
784} // namespace Pscf
Flory-Huggins excess free energy model.
Definition Interaction.h:26
virtual double fHelmholtz(Array< double > const &c) const
Compute excess Helmholtz free energy per monomer.
void setNMonomer(int nMonomer)
Set the number of monomer types.
int nPolymer() const
Get number of polymer species.
int nMonomer() const
Get number of monomer types.
int nSolvent() const
Get number of solvent (point particle) species.
SolventT & solvent(int id)
Get a solvent solver object.
PolymerT & polymer(int id)
Get a polymer solver object by non-const reference.
double phi() const
Get the overall volume fraction for this species.
Definition Species.h:149
int nBlock() const
Number of blocks.
PropagatorT & propagator(int blockId, int directionId)
Get the propagator for a specific block and direction (non-const).
double mu() const
Get the chemical potential for this species (units kT=1).
Definition Species.h:155
double length() const
Sum of the lengths of all blocks in the polymer (thread model).
int nx() const
Get number of spatial grid points, including both endpoints.
double spatialAverage(Field const &f) const
Compute spatial average of a field.
double innerProduct(Field const &f, Field const &g) const
Compute inner product of two real fields.
Command to compute properties of homogeneous reference system.
void output(int mode, std::ostream &out)
Output comparison to a homogeneous reference system.
void compute(int mode)
Compute properties of a homogeneous reference system.
Factory for subclasses of Iterator.
virtual int solve(bool isContinuation=false)=0
Iterate to solution.
void setDomain(Domain const &domain)
Create an association with the domain and allocate memory.
Descriptor and solver for a block polymer species.
MDE solver for one-direction of one block.
QFieldT const & q(int i) const
Return q-field at specified step.
int ns() const
Number of values of s (or slices), including head and tail.
Solver and descriptor for a solvent species.
Default Factory for subclasses of Sweep.
DArray< double > Field
Generic Field type.
Definition r1d/System.h:70
void compute()
Solve the modified diffusion equation once, without iteration.
void readParam()
Read input parameters from default param file.
double fHelmholtz() const
Get precomputed Helmholtz free energy per monomer / kT.
Definition r1d/System.h:674
double pressure() const
Get precomputed pressure x monomer volume kT.
Definition r1d/System.h:680
void writeC(std::string const &filename) const
Write concentration fields in symmetrized basis format.
int iterate(bool isContinuation=false)
Iteratively solve a SCFT problem.
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.
Interaction & interaction()
Get interaction (i.e., excess free energy) by reference.
Definition r1d/System.h:574
CField & cField(int monomerId)
Get chemical potential field for a specific monomer type.
Definition r1d/System.h:662
DArray< CField > & cFields()
Get array of all chemical potential fields.
Definition r1d/System.h:649
void computeFreeEnergy()
Compute free energy density and pressure for current fields.
void writeW(std::string const &filename) const
Write chemical potential fields in symmetrized basis format.
void readCommands()
Read commands from default command file.
Iterator & iterator()
Get the Iterator by reference.
Definition r1d/System.h:605
Domain & domain()
Get spatial domain (including grid info) by reference.
Definition r1d/System.h:592
void writeQAll(std::string const &basename) const
Write all propagators of all blocks, each to a separate file.
FileMaster & fileMaster()
Get FileMaster by reference.
Definition r1d/System.h:614
void setOptions(int argc, char **argv)
Process command line options.
WField & wField(int monomerId)
Get chemical potential field for a specific monomer type.
Definition r1d/System.h:635
DArray< WField > & wFields()
Get array of all chemical potential fields.
Definition r1d/System.h:621
void writeQ(std::string const &filename, int polymerId, int blockId, int directionId) const
Write one propagator for one block, in r-grid format.
void writeQTail(std::string const &filename, int polymerId, int blockId, int directionId) const
Write the final slice of a propagator in r-grid format.
virtual void readParameters(std::istream &in)
Read input parameters (without opening and closing lines).
void writeBlockC(std::string const &filename) const
Write c-fields for all blocks and solvents in r-grid format.
void sweep()
Sweep in parameter space, solving an SCF problem at each point.
void writeThermo(std::ostream &out) const
Write thermodynamic properties to a file.
void writeParamNoSweep(std::ostream &out) const
Write parameter file to an ostream, omitting any Sweep block.
~System()
Destructor.
Mixture & mixture()
Get Mixture by reference.
Definition r1d/System.h:562
System()
Constructor.
double size() const
Get the size (number of monomers) in this solvent.
double phi() const
Get the overall volume fraction for this species.
Definition Species.h:149
double mu() const
Get the chemical potential for this species (units kT=1).
Definition Species.h:155
void allocate(int capacity)
Allocate the underlying C array.
Definition DArray.h:199
Wrapper for a double precision number, for formatted ostream output.
Definition Dbl.h:40
void openOutputFile(const std::string &filename, std::ofstream &out, std::ios_base::openmode mode=std::ios_base::out) const
Open an output file.
void setParamFileName(const std::string &paramFileName)
Set the parameter file name.
void setOutputPrefix(const std::string &outputPrefix)
Set the output file prefix string.
void setCommandFileName(const std::string &commandFileName)
Set the command file name.
void setInputPrefix(const std::string &inputPrefix)
Set the input file prefix string.
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
static void setEcho(bool echo=true)
Enable or disable echoing for all subclasses of ParamComponent.
Begin & readBegin(std::istream &in, const char *label, bool isRequired=true)
Add and read a class label and opening bracket.
void setClassName(const char *className)
Set class name string.
virtual void writeParam(std::ostream &out) const
Write all parameters to an output stream.
std::string className() const
Get class name string.
void readParamComposite(std::istream &in, ParamComposite &child, bool next=true)
Add and read a required child ParamComposite.
End & readEnd(std::istream &in)
Add and read the closing bracket.
Wrapper for a std::string, for formatted ostream output.
Definition Str.h:37
#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
std::string toString(int n)
Return string representation of an integer.
Definition ioUtil.cpp:52
SCFT with real 1D fields.
PSCF package top-level namespace.
Definition param_pc.dox:1
void set(BracketPolicy::Type policy)
Set policy regarding use of bracket delimiters on arrays.