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>
18#include <pscf/inter/Interaction.h>
19#include <pscf/inter/Interaction.h>
20#include <pscf/floryHuggins/Clump.h>
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>
48 iteratorFactoryPtr_(0),
66 fieldIo_.associate(domain_, fileMaster_);
79 if (interactionPtr_) {
80 delete interactionPtr_;
85 if (iteratorFactoryPtr_) {
86 delete iteratorFactoryPtr_;
91 if (sweepFactoryPtr_) {
92 delete sweepFactoryPtr_;
114 while ((c = getopt(argc, argv,
"er:p:c:i:o:f")) != -1) {
136 Log::file() <<
"Unknown option -" << optopt << std::endl;
177 homogeneous_.initialize(
mixture());
189 iteratorPtr_ = iteratorFactoryPtr_->readObject(in, *
this,
192 std::string msg =
"Unrecognized Iterator subclass name ";
199 sweepFactoryPtr_->readObjectOptional(in, *
this,
className, isEnd);
202 sweepPtr_->setSystem(*
this);
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;
240 void System::allocateFields()
251 wFields_.allocate(nMonomer);
252 cFields_.allocate(nMonomer);
254 for (
int i = 0; i < nMonomer; ++i) {
270 std::string filename;
272 std::istream& inBuffer = in;
274 bool readNext =
true;
279 if (inBuffer.eof()) {
285 if (command ==
"FINISH") {
289 if (command ==
"READ_W") {
290 readEcho(inBuffer, filename);
291 fieldIo_.readFields(
wFields(), filename);
293 if (command ==
"COMPUTE") {
298 if (command ==
"ITERATE") {
300 bool isContinuation =
false;
301 int fail =
iterate(isContinuation);
306 if (command ==
"SWEEP") {
311 if (command ==
"WRITE_PARAM") {
312 readEcho(inBuffer, filename);
318 if (command ==
"WRITE_THERMO") {
319 readEcho(inBuffer, filename);
325 if (command ==
"COMPARE_HOMOGENEOUS") {
329 readEcho(inBuffer, filename);
331 Log::file() <<
"mode = " << mode << std::endl;
338 comparison.
output(mode, file);
341 if (command ==
"WRITE_W") {
342 readEcho(inBuffer, filename);
343 fieldIo_.writeFields(
wFields(), filename);
345 if (command ==
"WRITE_C") {
346 readEcho(inBuffer, filename);
347 fieldIo_.writeFields(
cFields(), filename);
349 if (command ==
"WRITE_BLOCK_C") {
350 readEcho(inBuffer, filename);
351 fieldIo_.writeBlockCFields(mixture_, filename);
353 if (command ==
"WRITE_Q_SLICE") {
354 int polymerId, blockId, directionId, segmentId;
355 readEcho(inBuffer, filename);
356 inBuffer >> polymerId;
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,
367 if (command ==
"WRITE_Q_TAIL") {
368 readEcho(inBuffer, filename);
369 int polymerId, blockId, directionId;
370 inBuffer >> polymerId;
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);
378 if (command ==
"WRITE_Q_VERTEX") {
379 int polymerId, vertexId;
380 inBuffer >> polymerId;
383 <<
Int(polymerId, 5) << std::endl;
384 inBuffer >> vertexId;
386 <<
Int(vertexId, 5) << std::endl;
387 inBuffer >> filename;
389 <<
Str(filename, 20) << std::endl;
390 fieldIo_.writeVertexQ(mixture_, polymerId, vertexId, filename);
392 if (command ==
"WRITE_Q") {
393 readEcho(inBuffer, filename);
394 int polymerId, blockId, directionId;
395 inBuffer >> polymerId;
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);
403 if (command ==
"WRITE_Q_ALL") {
404 readEcho(inBuffer, filename);
407 if (command ==
"REMESH_W") {
412 inBuffer >> filename;
413 Log::file() <<
"outfile = " <<
Str(filename, 20) << std::endl;
414 fieldIo_.remesh(
wFields(), nx, filename);
416 if (command ==
"EXTEND_W") {
421 inBuffer >> filename;
422 Log::file() <<
"outfile = " <<
Str(filename, 20) << std::endl;
423 fieldIo_.extend(
wFields(), m, filename);
425 Log::file() <<
" Error: Unknown command "
426 << command << std::endl;
447 void System::readEcho(std::istream& in, std::string&
string)
const
463 mixture_.compute(
wFields(), cFields_);
521 for (
int i = 0; i < np; ++i) {
523 phi = polymerPtr->
phi();
524 mu = polymerPtr->
mu();
526 length = polymerPtr->
length();
528 fHelmholtz_ += phi*( mu - 1.0 )/length;
537 for (
int i = 0; i < ns; ++i) {
539 phi = solventPtr->
phi();
540 mu = solventPtr->
mu();
542 size = solventPtr->
size();
544 fHelmholtz_ += phi*( mu - 1.0 )/size;
550 for (
int i = 0; i < nm; ++i) {
554 fIdeal_ = fHelmholtz_;
558 if (!f_.isAllocated()) f_.allocate(nx);
559 if (!c_.isAllocated()) c_.allocate(nm);
561 for (
int i = 0; i < nx; ++i) {
563 for (j = 0; j < nm; ++j) {
564 c_[j] = cFields_[j][i];
570 fInter_ = fHelmholtz_ - fIdeal_;
573 pressure_ = -fHelmholtz_;
577 for (
int i = 0; i < np; ++i) {
579 phi = polymerPtr->
phi();
580 mu = polymerPtr->
mu();
581 length = polymerPtr->
length();
583 pressure_ += phi*mu/length;
590 for (
int i = 0; i < ns; ++i) {
592 phi = solventPtr->
phi();
593 mu = solventPtr->
mu();
594 size = solventPtr->
size();
596 pressure_ += phi*mu/size;
606 out <<
"fHelmholtz " <<
Dbl(
fHelmholtz(), 18, 11) << std::endl;
607 out <<
"pressure " <<
Dbl(
pressure(), 18, 11) << std::endl;
609 out <<
"fIdeal " <<
Dbl(fIdeal_, 18, 11) << std::endl;
610 out <<
"fInter " <<
Dbl(fInter_, 18, 11) << std::endl;
616 out <<
"polymers:" << std::endl;
621 for (
int i = 0; i < np; ++i) {
623 <<
" " <<
Dbl(
mixture().polymer(i).phi(),18, 11)
624 <<
" " <<
Dbl(
mixture().polymer(i).mu(), 18, 11)
633 out <<
"solvents:" << std::endl;
638 for (
int i = 0; i < ns; ++i) {
640 <<
" " <<
Dbl(
mixture().solvent(i).phi(),18, 11)
641 <<
" " <<
Dbl(
mixture().solvent(i).mu(), 18, 11)
657 fieldIo_.writeFields(
wFields(), filename);
666 fieldIo_.writeFields(
cFields(), filename);
676 fieldIo_.writeBlockCFields(mixture_, filename);
683 int polymerId,
int blockId,
684 int directionId,
int segmentId)
689 Polymer const& polymer = mixture_.polymer(polymerId);
695 propagator = polymer.
propagator(blockId, directionId);
696 Field const& field = propagator.
q(segmentId);
697 fieldIo_.writeField(field, filename);
704 int polymerId,
int blockId,
int directionId)
709 Polymer const& polymer = mixture_.polymer(polymerId);
715 field = polymer.
propagator(blockId, directionId).tail();
716 fieldIo_.writeField(field, filename);
723 int polymerId,
int blockId,
int directionId)
728 Polymer const& polymer = mixture_.polymer(polymerId);
735 int nslice = propagator.
ns();
739 fileMaster_.openOutputFile(filename, file);
742 file <<
"ngrid" << std::endl
743 <<
" " << domain_.nx() << std::endl
744 <<
"nslice" << std::endl
745 <<
" " << nslice << std::endl;
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);
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) {
777 writeQ(filename, ip, ib,
id);
Flory-Huggins excess free energy model.
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.
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).
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.
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.
double pressure() const
Get precomputed pressure x monomer volume kT.
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.
CField & cField(int monomerId)
Get chemical potential field for a specific monomer type.
DArray< CField > & cFields()
Get array of all chemical potential fields.
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.
Domain & domain()
Get spatial domain (including grid info) by reference.
void writeQAll(std::string const &basename) const
Write all propagators of all blocks, each to a separate file.
FileMaster & fileMaster()
Get FileMaster by reference.
void setOptions(int argc, char **argv)
Process command line options.
WField & wField(int monomerId)
Get chemical potential field for a specific monomer type.
DArray< WField > & wFields()
Get array of all chemical potential fields.
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.
Mixture & mixture()
Get Mixture by reference.
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).
void allocate(int capacity)
Allocate the underlying C array.
Wrapper for a double precision number, for formatted ostream output.
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 ¶mFileName)
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.
static std::ostream & file()
Get log ostream by reference.
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.
#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.
SCFT with real 1D fields.
PSCF package top-level namespace.
void set(BracketPolicy::Type policy)
Set policy regarding use of bracket delimiters on arrays.