rp/system/System.tpp

#ifndef RP_SYSTEM_TPP
#define RP_SYSTEM_TPP
 
/*
* PSCF - Polymer Self-Consistent Field
*
* Copyright 2015 - 2025, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/
 
#include "System.h"
 
#include <prdc/crystal/Basis.h>
#include <prdc/crystal/UnitCell.h>
 
#include <pscf/interaction/Interaction.h>
#include <pscf/mesh/Mesh.h>
#include <pscf/math/IntVec.h>
 
#include <util/containers/DArray.h>
#include <util/containers/FSArray.h>
#include <util/param/BracketPolicy.h>
#include <util/param/ParamComponent.h>
#include <util/signal/Signal.h>
#include <util/format/Str.h>
#include <util/format/Int.h>
#include <util/format/Dbl.h>
#include <util/misc/ioUtil.h>
#include <util/misc/FileMaster.h>
 
#include <string>
#include <unistd.h>
 
namespace Pscf {
namespace Rp {
 
   using namespace Util;
   using namespace Prdc;
 
   /*
   * Constructor.
   */
   template <int D, class T>
   System<D,T>::System(typename T::System& system)
    : systemPtr_(&system),
      cPtr_(nullptr),
      wPtr_(nullptr),
      hPtr_(nullptr),
      maskPtr_(nullptr),
      mixturePtr_(nullptr),
      mixtureModifierPtr_(nullptr),
      domainPtr_(nullptr),
      interactionPtr_(nullptr),
      scftPtr_(nullptr),
      environmentPtr_(nullptr),
      environmentFactoryPtr_(nullptr),
      iteratorPtr_(nullptr),
      iteratorFactoryPtr_(nullptr),
      sweepPtr_(nullptr),
      sweepFactoryPtr_(nullptr),
      simulatorPtr_(nullptr),
      simulatorFactoryPtr_(nullptr),
      fileMasterPtr_(nullptr),
      tmpUnitCellPtr_(nullptr),
      polymerModel_(PolymerModel::Thread),
      isAllocatedGrid_(false),
      isAllocatedBasis_(false),
      hasMixture_(false)
   {
      setClassName("System");  // Set block label in parameter file
      BracketPolicy::set(BracketPolicy::Optional);
 
      // Create field containers
      cPtr_ = new typename T::CFields();
      wPtr_ = new typename T::WFields();
      hPtr_ = new typename T::WFields();
      maskPtr_ = new typename T::Mask();
 
      // Create non-optional components
      mixturePtr_ = new typename T::Mixture();
      mixtureModifierPtr_ = new typename T::MixtureModifier();
      interactionPtr_ = new typename T::Interaction();
      domainPtr_ = new typename T::Domain();
      scftPtr_ = new typename T::ScftThermo(*systemPtr_);
      fileMasterPtr_ = new FileMaster();
      tmpUnitCellPtr_ = new UnitCell<D>();
 
      // Create Factories for optional components
      environmentFactoryPtr_ 
         = new typename T::EnvironmentFactory(*systemPtr_);
      iteratorFactoryPtr_ = new typename T::IteratorFactory(*systemPtr_);
      sweepFactoryPtr_ = new typename T::SweepFactory(*systemPtr_);
      simulatorFactoryPtr_ = new typename T::SimulatorFactory(*systemPtr_);
 
      // Create associations among class members
      mixtureModifier().associate(mixture_());
      domain_().setFileMaster(*fileMasterPtr_);
      w().setFieldIo(domain_().fieldIo());
      w().setReadUnitCell(domain_().unitCell());
      w().setWriteUnitCell(domain_().unitCell());
      h().setFieldIo(domain_().fieldIo());
      h().setReadUnitCell(*tmpUnitCellPtr_);
      h().setWriteUnitCell(domain_().unitCell());
      mask().setFieldIo(domain_().fieldIo());
      mask().setReadUnitCell(*tmpUnitCellPtr_);
      mask().setWriteUnitCell(domain_().unitCell());
      c_().setFieldIo(domain_().fieldIo());
      c_().setWriteUnitCell(domain_().unitCell());
 
      // Note the arguments of setReadUnitCell functions used above:
      // When w is read from a file in basis or r-grid format, the
      // parameters of the system unit cell, domain_().unitCell(), are
      // reset to those in the field file header. When imposed fields h()
      // and mask() are read from a file, however, unit cell parameters
      // from the field file header are read into a mutable workspace
      // object, *tmpUnitCellPtr_, and ultimately discarded.
 
      // Use of "signals" to maintain data relationships:
      // Signals are instances of Unit::Signal<void> that are used to
      // notify "observer" objects of modification of data owned by a
      // related "notifier" object. Each signal is owned by a notifier
      // object that maintains data that may be modified. Each signal
      // maintains a list of observers objects that should be notified
      // whenever the data owned by the notifier object changes. Each
      // observer is added by the Signal<void>::addObserver function
      // template, which takes two arguments: a reference to an observer
      // object (i.e., an instance of a class) and a pointer to a member
      // function of that class which will be invoked on that object when
      // the signal is triggered by modification of associated data.
 
      // Addition of observers to signals
 
      // Signal triggered by unit cell modification
      Signal<void>& cellSignal = domain_().unitCell().signal();
      cellSignal.addObserver(*this, &System<D,T>::clearUnitCellData);
 
      // Signal triggered by basis construction
      Signal<void>& basisSignal = domain_().basis().signal();
      basisSignal.addObserver(*this, &System<D,T>::allocateFieldsBasis);
 
      // Signal triggered by w-field modification
      w().signal().addObserver(*this, &System<D,T>::clearCFields);
 
      // Signal triggered by h-field modification
      h().signal().addObserver(*this, &System<D,T>::clearCFields);
 
      // Signal triggered by mask modification
      mask().signal().addObserver(*this, &System<D,T>::clearCFields);
   }
 
   /*
   * Destructor.
   */
   template <int D, class T>
   System<D,T>::~System()
   {
      // Delete optional elements
      if (environmentPtr_) {
         delete environmentPtr_;
      }
      if (iteratorPtr_) {
         delete iteratorPtr_;
      }
      if (sweepPtr_) {
         delete sweepPtr_;
      }
      if (simulatorPtr_) {
         delete simulatorPtr_;
      }
 
      // Delete factories
      delete environmentFactoryPtr_;
      delete iteratorFactoryPtr_;
      delete sweepFactoryPtr_;
      delete simulatorFactoryPtr_;
 
      // Delete other non-optional components
      delete mixturePtr_;
      delete mixtureModifierPtr_;
      delete interactionPtr_;
      delete domainPtr_;
      delete scftPtr_;
      delete fileMasterPtr_;
      delete tmpUnitCellPtr_;
 
      // Delete field containers
      delete maskPtr_;
      delete hPtr_;
      delete wPtr_;
      delete cPtr_;
   }
 
   // Lifetime (called in main program)
 
   /*
   * Process command line options.
   */
   template <int D, class T>
   void System<D,T>::setOptions(int argc, char **argv)
   {
      bool eFlag = false;  // echo
      bool dFlag = false;  // spatial dimension (1, 2, or 3)
      bool pFlag = false;  // param file
      bool cFlag = false;  // command file
      bool iFlag = false;  // input prefix
      bool oFlag = false;  // output prefix
      bool tFlag = false;  // thread count
      char* pArg = 0;
      char* cArg = 0;
      char* iArg = 0;
      char* oArg = 0;
      int dArg = 0;
      int tArg = 0;
 
      // Read program arguments
      int c;
      opterr = 0;
      while ((c = getopt(argc, argv, "ed:p:c:i:o:t:")) != -1) {
         switch (c) {
         case 'e':
            eFlag = true;
            break;
         case 'd':
            dFlag = true;
            dArg = atoi(optarg);
            break;
         case 'p': // parameter file
            pFlag = true;
            pArg  = optarg;
            break;
         case 'c': // command file
            cFlag = true;
            cArg  = optarg;
            break;
         case 'i': // input prefix
            iFlag = true;
            iArg  = optarg;
            break;
         case 'o': // output prefix
            oFlag = true;
            oArg  = optarg;
            break;
         case 't': // thread count, if set by user
            tFlag = true;
            tArg = atoi(optarg);
            break;
         case '?':
           Log::file() << "Unknown option -" << optopt << std::endl;
           UTIL_THROW("Invalid command line option");
         }
      }
 
      // Process program arguments
 
      // Check required option -d
      if (dFlag) {
         UTIL_CHECK(D == dArg);
      } else {
         UTIL_THROW("Missing required -d option");
      }
 
      // Check required option -p, set parameter file name
      if (pFlag) {
         fileMaster().setParamFileName(std::string(pArg));
      } else {
         UTIL_THROW("Missing required -p option - no parameter file");
      }
 
      // Check required option -c, set command file name
      if (cFlag) {
         fileMaster().setCommandFileName(std::string(cArg));
      } else {
         UTIL_THROW("Missing required -c option - no command file");
      }
 
      // If option -e, set to echo parameters as they are read
      if (eFlag) {
         Util::ParamComponent::setEcho(true);
      }
 
      // If option -i, set path prefix for input files
      if (iFlag) {
         fileMaster().setInputPrefix(std::string(iArg));
      }
 
      // If option -o, set path prefix for output files
      if (oFlag) {
         fileMaster().setOutputPrefix(std::string(oArg));
      }
 
      // If option -t, process the thread count
      if (tFlag) {
         if (tArg <= 0) {
            UTIL_THROW("Error: Non-positive thread count -t option");
         }
         setThreadCount(tArg);
      }
   }
 
   /*
   * Read parameters and initialize.
   */
   template <int D, class T>
   void System<D,T>::readParameters(std::istream& in)
   {
      // Optionally read polymerModel_ enum value, set the global value
      if (!PolymerModel::isLocked()) {
         polymerModel_ = PolymerModel::Thread;
         readOptional(in, "polymerModel", polymerModel_);
 
         // Set the global enumeration value
         PolymerModel::setModel(polymerModel_);
      }
 
      // Check number of times the global polymer model has been set.
      // Retain this value so we can check at the end of this function
      // that it was not set again within the function.
 
      int nSetPolymerModel = PolymerModel::nSet();
 
      // Read the Mixture{ ... } block
      readParamComposite(in, mixture_());
      hasMixture_ = true;
 
      int nm = mixture_().nMonomer();
      int np = mixture_().nPolymer();
      int ns = mixture_().nSolvent();
      UTIL_CHECK(nm > 0);
      UTIL_CHECK(np >= 0);
      UTIL_CHECK(ns >= 0);
      UTIL_CHECK(np + ns > 0);
 
      // Read the Interaction{ ... } block
      interaction().setNMonomer(nm);
      readParamComposite(in, interaction());
 
      // Read the Domain{ ... } block
      readParamComposite(in, domain_());
      UTIL_CHECK(domain_().mesh().size() > 0);
      UTIL_CHECK(domain_().unitCell().nParameter() > 0);
      UTIL_CHECK(domain_().unitCell().lattice() != UnitCell<D>::Null);
      domain_().fieldIo().setNMonomer(nm);
 
      // Setup the mixture
      mixture_().associate(domain_().mesh(), domain_().fft(),
                         domain_().unitCell(), domain_().waveList());
      mixture_().setFieldIo(domain_().fieldIo());
      mixture_().allocate();
 
      // Allocate memory for field containers in r-grid form
      allocateFieldsGrid();
 
      // Optionally read and construct an Environment
      std::string className;
      bool isEnd;
      environmentPtr_ =
         environmentFactoryPtr_->readObjectOptional(in, *this, 
                                                    className, isEnd);
      if (!environmentPtr_ && ParamComponent::echo()) {
         Log::file() << indent() << "  Environment{ [absent] }\n";
      }
      if (environmentPtr_) {
         environment().setNParams(domain_().unitCell().nParameter());
      }
 
      // Optionally read and construct an Iterator
      if (!isEnd) {
         iteratorPtr_ =
            iteratorFactoryPtr_->readObjectOptional(in, *this, 
                                                    className, isEnd);
         if (iteratorPtr_) {
            if (iteratorPtr_->isSymmetric()) {
               UTIL_CHECK(domain_().hasGroup()); 
            }
         } else {
            if (ParamComponent::echo()) {
               Log::file() << indent() << "  Iterator{ [absent] }\n";
            }
         }
      }
 
      // Optionally read and construct a Sweep (if an iterator exists)
      if (hasIterator() && !isEnd) {
         sweepPtr_ =
            sweepFactoryPtr_->readObjectOptional(in, *this,
                                                 className, isEnd);
         if (!sweepPtr_ && ParamComponent::echo()) {
            Log::file() << indent() << "  Sweep{ [absent] }\n";
         }
      }
 
      // Optionally read and construct a Simulator 
      if (!isEnd) {
         simulatorPtr_ =
            simulatorFactoryPtr_->readObjectOptional(in, *this,
                                                     className, isEnd);
         if (!simulatorPtr_ && ParamComponent::echo()) {
            Log::file() << indent() << "  Simulator{ [absent] }\n";
         }
      }
 
      // Check that the polymer model was not modified after initialization
      UTIL_CHECK(PolymerModel::nSet() == nSetPolymerModel);
 
      // Note: The main pscf_rpc or pscf_rpg program can irreversibly lock the
      // global polymer model after reading the parameter file, to also prohibit 
      // later changes. The model is not locked here to allow repeated calls to
      // this function during unit testing.
   }
 
   /*
   * Read parameter file (including open and closing brackets).
   */
   template <int D, class T>
   void System<D,T>::readParam(std::istream& in)
   {
      readBegin(in, className().c_str());
      readParameters(in);
      readEnd(in);
   }
 
   /*
   * Read default parameter file.
   */
   template <int D, class T>
   void System<D,T>::readParam()
   {  readParam(fileMaster().paramFile()); }
 
   /*
   * Read and execute commands from a specified command file.
   */
   template <int D, class T>
   void System<D,T>::readCommands(std::istream &in)
   {
      UTIL_CHECK(isAllocatedGrid_);
      std::string command, filename, inFileName, outFileName;
      typename T::FieldIo const & fieldIo = domain_().fieldIo();
 
      bool readNext = true;
      while (readNext) {
 
         in >> command;
 
         if (in.eof()) {
            break;
         } else {
            Log::file() << command << std::endl;
         }
 
         if (command == "FINISH") {
            Log::file() << std::endl;
            readNext = false;
         } else
         if (command == "READ_W_BASIS") {
            readEcho(in, filename);
            w().readBasis(filename);
            UTIL_CHECK(domain_().unitCell().isInitialized());
            UTIL_CHECK(domain_().basis().isInitialized());
            UTIL_CHECK(!domain_().waveList().hasKSq());
            UTIL_CHECK(isAllocatedBasis_);
            UTIL_CHECK(!c().hasData());
            UTIL_CHECK(!scft().hasData());
         } else
         if (command == "READ_W_RGRID") {
            readEcho(in, filename);
            w().readRGrid(filename);
            UTIL_CHECK(domain_().unitCell().isInitialized());
            UTIL_CHECK(!domain_().waveList().hasKSq());
            UTIL_CHECK(!c().hasData());
            UTIL_CHECK(!scft().hasData());
         } else
         if (command == "SET_UNIT_CELL") {
            UnitCell<D> unitCell;
            in >> unitCell;
            Log::file() << "   " << unitCell << std::endl;
            setUnitCell(unitCell);
            UTIL_CHECK(domain_().unitCell().isInitialized());
            UTIL_CHECK(!domain_().waveList().hasKSq());
            UTIL_CHECK(!c().hasData());
            UTIL_CHECK(!scft().hasData());
         } else
         if (command == "COMPUTE") {
            // Solve the modified diffusion equation, without iteration
            compute();
         } else
         if (command == "ITERATE") {
            // Attempt to iteratively solve a single SCFT problem
            bool isContinuation = false;
            int fail = iterate(isContinuation);
            if (fail) {
               readNext = false;
            }
         } else
         if (command == "SWEEP") {
            // Attempt to solve a sequence of SCFT problems along a path
            // through parameter space
            sweep();
         } else
         if (command == "COMPRESS") {
            // Impose incompressibility
            UTIL_CHECK(hasSimulator());
            simulator().compressor().compress();
         } else
         if (command == "SIMULATE") {
            // Perform a field theoretic simulation
            int nStep;
            in >> nStep;
            Log::file() << "   "  << nStep << "\n";
            simulate(nStep);
         } else
         if (command == "ANALYZE" || command == "ANALYZE_TRAJECTORY") {
            // Read and analyze a field trajectory file
            int min, max;
            in >> min;
            in >> max;
            Log::file() << "   "  << min ;
            Log::file() << "   "  << max << "\n";
            std::string classname;
            readEcho(in, classname);
            readEcho(in, filename);
            simulator().analyze(min, max, classname, filename);
         } else
         if (command == "WRITE_PARAM") {
            readEcho(in, filename);
            std::ofstream file;
            fileMaster().openOutputFile(filename, file);
            writeParamNoSweep(file);
            file.close();
         } else
         if (command == "WRITE_THERMO") {
            readEcho(in, filename);
            std::ofstream file;
            fileMaster().openOutputFile(filename, file,
                                       std::ios_base::app);
            scft().write(file);
            file.close();
         } else
         if (command == "WRITE_STRESS") {
            readEcho(in, filename);
            std::ofstream file;
            if (!mixture_().hasStress()) {
               computeStress();
            }
            fileMaster().openOutputFile(filename, file,
                                       std::ios_base::app);
            mixture_().writeStress(file);
            if (hasEnvironment()) {
               environment().writeStress(file);
            }
            file.close();
         } else
         if (command == "WRITE_W_BASIS") {
            readEcho(in, filename);
            w().writeBasis(filename);
         } else
         if (command == "WRITE_W_RGRID") {
            readEcho(in, filename);
            w().writeRGrid(filename);
         } else
         if (command == "WRITE_C_BASIS") {
            readEcho(in, filename);
            c().writeBasis(filename);
         } else
         if (command == "WRITE_C_RGRID") {
            readEcho(in, filename);
            c().writeRGrid(filename);
         } else
         if (command == "WRITE_BLOCK_C_RGRID") {
            readEcho(in, filename);
            mixture_().writeBlockCRGrid(filename);
         } else
         if (command == "WRITE_Q_SLICE") {
            int polymerId, blockId, directionId, segmentId;
            readEcho(in, filename);
            in >> polymerId;
            in >> blockId;
            in >> directionId;
            in >> segmentId;
            Log::file() << Str("polymer ID  ", 21) << polymerId << "\n"
                        << Str("block ID  ", 21) << blockId << "\n"
                        << Str("direction ID  ", 21) << directionId
                        << "\n"
                        << Str("segment ID  ", 21) << segmentId
                        << std::endl;
            mixture_().writeQSlice(filename, polymerId, blockId,
                                 directionId, segmentId);
         } else
         if (command == "WRITE_Q_TAIL") {
            readEcho(in, filename);
            int polymerId, blockId, directionId;
            in >> polymerId;
            in >> blockId;
            in >> directionId;
            Log::file() << Str("polymer ID  ", 21) << polymerId << "\n"
                        << Str("block ID  ", 21) << blockId << "\n"
                        << Str("direction ID  ", 21) << directionId
                        << "\n";
            mixture_().writeQTail(filename, polymerId, blockId, directionId);
         } else
         if (command == "WRITE_Q") {
            readEcho(in, filename);
            int polymerId, blockId, directionId;
            in >> polymerId;
            in >> blockId;
            in >> directionId;
            Log::file() << Str("polymer ID  ", 21) << polymerId << "\n"
                        << Str("block ID  ", 21) << blockId << "\n"
                        << Str("direction ID  ", 21) << directionId
                        << "\n";
            mixture_().writeQ(filename, polymerId, blockId, directionId);
         } else
         if (command == "WRITE_Q_ALL") {
            readEcho(in, filename);
            mixture_().writeQAll(filename);
         } else
         if (command == "WRITE_STARS") {
            readEcho(in, filename);
            domain_().writeStars(filename);
         } else
         if (command == "WRITE_WAVES") {
            readEcho(in, filename);
            domain_().writeWaves(filename);
         } else
         if (command == "WRITE_GROUP") {
            readEcho(in, filename);
            domain_().writeGroup(filename);
         } else
         if (command == "BASIS_TO_RGRID") {
            readEcho(in, inFileName);
            readEcho(in, outFileName);
            fieldIo.convertBasisToRGrid(inFileName, outFileName);
         } else
         if (command == "RGRID_TO_BASIS") {
            readEcho(in, inFileName);
            readEcho(in, outFileName);
            fieldIo.convertRGridToBasis(inFileName, outFileName);
         } else
         if (command == "KGRID_TO_RGRID") {
            readEcho(in, inFileName);
            readEcho(in, outFileName);
            fieldIo.convertKGridToRGrid(inFileName, outFileName);
         } else
         if (command == "RGRID_TO_KGRID") {
            readEcho(in, inFileName);
            readEcho(in, outFileName);
            fieldIo.convertRGridToKGrid(inFileName, outFileName);
         } else
         if (command == "BASIS_TO_KGRID") {
            readEcho(in, inFileName);
            readEcho(in, outFileName);
            fieldIo.convertBasisToKGrid(inFileName, outFileName);
         } else
         if (command == "KGRID_TO_BASIS") {
            readEcho(in, inFileName);
            readEcho(in, outFileName);
            fieldIo.convertKGridToBasis(inFileName, outFileName);
         } else
         if (command == "CHECK_RGRID_SYMMETRY") {
            double epsilon;
            readEcho(in, inFileName);
            readEcho(in, epsilon);
            bool hasSymmetry;
            hasSymmetry = fieldIo.hasSymmetry(inFileName, epsilon);
            if (hasSymmetry) {
               Log::file() << std::endl
                   << "Symmetry of r-grid file matches this space group"
                   << std::endl << std::endl;
            } else {
               Log::file() << std::endl
                   << "Symmetry of r-grid file does not match this\n"
                   << "space group to within error threshold of "
                   << Dbl(epsilon) << " ." << std::endl << std::endl;
            }
         } else
         if (command == "COMPARE_BASIS") {
            std::string filecompare1, filecompare2;
            readEcho(in, filecompare1);
            readEcho(in, filecompare2);
            fieldIo.compareFieldsBasis(filecompare1, filecompare2);
         } else
         if (command == "COMPARE_RGRID") {
            std::string filecompare1, filecompare2;
            readEcho(in, filecompare1);
            readEcho(in, filecompare2);
            fieldIo.compareFieldsRGrid(filecompare1, filecompare2);
         } else
         if (command == "SCALE_BASIS") {
            double factor;
            readEcho(in, inFileName);
            readEcho(in, outFileName);
            readEcho(in, factor);
            fieldIo.scaleFieldsBasis(inFileName, outFileName, factor);
         } else
         if (command == "SCALE_RGRID") {
            double factor;
            readEcho(in, inFileName);
            readEcho(in, outFileName);
            readEcho(in, factor);
            fieldIo.scaleFieldsRGrid(inFileName, outFileName, factor);
         } else
         if (command == "EXPAND_RGRID_DIMENSION") {
            readEcho(in, inFileName);
            readEcho(in, outFileName);
 
            // Read expanded spatial dimension d
            int d;
            in >> d;
            UTIL_CHECK(d > D);
            Log::file() << Str("Expand dimension to: ") << d << "\n";
 
            // Read numbers of grid points along added dimensions
            DArray<int> newGridDimensions;
            newGridDimensions.allocate(d-D);
            for (int i = 0; i < d-D; i++) {
               in >> newGridDimensions[i];
            }
            Log::file()
               << Str("Numbers of grid points in added dimensions :  ");
            for (int i = 0; i < d-D; i++) {
               Log::file() << newGridDimensions[i] << " ";
            }
            Log::file() << "\n";
 
            fieldIo.expandRGridDimension(inFileName, outFileName,
                                         d, newGridDimensions);
 
         } else
         if (command == "REPLICATE_UNIT_CELL") {
            readEcho(in, inFileName);
            readEcho(in, outFileName);
 
            // Read numbers of replicas along each direction
            IntVec<D> replicas;
            for (int i = 0; i < D; i++){
               in >> replicas[i];
            }
            for (int i = 0; i < D; i++){
               Log::file() << "Replicate unit cell in direction "
                           << i << " : ";
               Log::file() << replicas[i] << " times ";
               Log::file() << "\n";
            }
            fieldIo.replicateUnitCell(inFileName, outFileName, replicas);
 
         } else
         if (command == "ESTIMATE_W_BASIS") {
            readEcho(in, inFileName);
            readEcho(in, outFileName);
            fieldIo.estimateWBasis(inFileName, outFileName,
                                   interaction().chi());
         } else
         if (command == "READ_H_BASIS") {
            readEcho(in, filename);
            h().readBasis(filename);
            UTIL_CHECK(!c().hasData());
         } else
         if (command == "READ_H_RGRID") {
            readEcho(in, filename);
            h().readRGrid(filename);
            UTIL_CHECK(!c().hasData());
         } else
         if (command == "WRITE_H_BASIS") {
            readEcho(in, filename);
            h().writeBasis(filename);
         } else
         if (command == "WRITE_H_RGRID") {
            readEcho(in, filename);
            h().writeRGrid(filename);
         } else
         if (command == "READ_MASK_BASIS") {
            readEcho(in, filename);
            mask().readBasis(filename);
         } else
         if (command == "READ_MASK_RGRID") {
            readEcho(in, filename);
            mask().readRGrid(filename);
         } else
         if (command == "WRITE_MASK_BASIS") {
            readEcho(in, filename);
            mask().writeBasis(filename);
         } else
         if (command == "WRITE_MASK_RGRID") {
            readEcho(in, filename);
            mask().writeRGrid(filename);
         } else
         if (command == "WRITE_TIMERS") {
            readEcho(in, filename);
            std::ofstream file;
            fileMaster().openOutputFile(filename, file);
            writeTimers(file);
            file.close();
         } else
         if (command == "CLEAR_TIMERS") {
            clearTimers();
         } else {
            Log::file() << "Error: Unknown command  "
                        << command << std::endl;
            readNext = false;
         }
      }
   }
 
   /*
   * Read and execute commands from the default command file.
   */
   template <int D, class T>
   void System<D,T>::readCommands()
   {
      if (fileMaster().commandFileName().empty()) {
         UTIL_THROW("Empty command file name");
      }
      readCommands(fileMaster().commandFile());
   }
 
   // Primary Field Theory Computations
 
   /*
   * Solve MDE for current w fields, without iteration.
   */
   template <int D, class T>
   void System<D,T>::compute(bool needStress)
   {
      // Preconditions
      UTIL_CHECK(isAllocatedGrid_);
      UTIL_CHECK(w().isAllocatedRGrid());
      UTIL_CHECK(c().isAllocatedRGrid());
      UTIL_CHECK(w().hasData());
      clearCFields();
 
      // Make sure Environment is up-to-date
      if (hasEnvironment()) {
         if (environment().needsUpdate()) {
            environment().generate();
         }
      }
 
      // Solve the modified diffusion equation (without iteration)
      mixture_().compute(w().rgrid(), c_().rgrid(), mask().phiTot());
      mixture_().setIsSymmetric(w().isSymmetric());
      c_().setHasData(true);
      scft().clear();
 
      // If w fields are symmetric, compute basis components for c fields
      if (w().isSymmetric()) {
         UTIL_CHECK(c().isAllocatedBasis());
         domain_().fieldIo().convertRGridToBasis(c().rgrid(), c_().basis(),
                                               false);
         c_().setIsSymmetric(true);
      }
 
      // Compute stress if needed
      if (needStress) {
         computeStress();
      }
   }
 
   /*
   * Compute SCFT stress for current fields.
   */
   template <int D, class T>
   void System<D,T>::computeStress()
   {
      // Compute and store standard Mixture stress
      if (!mixture_().hasStress()) {
         mixture_().computeStress(mask().phiTot());
      }
 
      // If necessary, compute and store Environment stress
      if (hasEnvironment()) {
         if (!environment().hasStress()) {
            environment().computeStress(mixture_(),
                                        iterator().flexibleParams());
         }
      }
   }
 
   /*
   * Iteratively solve a SCFT problem.
   */
   template <int D, class T>
   int System<D,T>::iterate(bool isContinuation)
   {
      // Preconditions
      UTIL_CHECK(hasIterator());
      UTIL_CHECK(w().hasData());
      if (iterator().isSymmetric()) {
         UTIL_CHECK(isAllocatedBasis_);
         UTIL_CHECK(w().isSymmetric());
      }
      clearCFields();
 
      Log::file() << std::endl;
      Log::file() << std::endl;
 
      // Make sure Environment is up-to-date
      if (hasEnvironment()) {
         if (environment().needsUpdate()) {
            environment().generate();
         }
      }
 
      // Call iterator (return 0 for convergence, 1 for failure)
      int error = iterator().solve(isContinuation);
      UTIL_CHECK(c().hasData());
 
      // If converged, compute related thermodynamic properties
      if (!error) {
         scft().compute();
         scft().write(Log::file());
         if (!iterator().isFlexible()) {
            if (!mixture_().hasStress()) {
               computeStress();
            }
            mixture_().writeStress(Log::file());
            if (hasEnvironment()) {
               environment().writeStress(Log::file());
            }
         }
      }
 
      return error;
   }
 
   /*
   * Perform an SCFT sweep along a path in parameter space.
   */
   template <int D, class T>
   void System<D,T>::sweep()
   {
      // Preconditions
      UTIL_CHECK(hasIterator());
      UTIL_CHECK(hasSweep());
      UTIL_CHECK(w().hasData());
      if (iterator().isSymmetric()) {
         UTIL_CHECK(isAllocatedBasis_);
         UTIL_CHECK(w().isSymmetric());
      }
 
      Log::file() << std::endl;
      Log::file() << std::endl;
 
      // Perform sweep
      sweepPtr_->sweep();
   }
 
   /*
   * Perform a field theoretic simulation of nStep steps.
   */
   template <int D, class T>
   void System<D,T>::simulate(int nStep)
   {
      UTIL_CHECK(nStep > 0);
      UTIL_CHECK(hasSimulator());
      clearCFields();
 
      simulator().simulate(nStep);
   }
 
   /*
   * Mark c-fields and free energy as outdated/invalid.
   */
   template <int D, class T>
   void System<D,T>::clearCFields()
   {
      c_().setHasData(false);
      scft().clear();
   }
 
   // Unit Cell Modifiers
 
   /*
   * Set the system unit cell.
   */
   template <int D, class T>
   void System<D,T>::setUnitCell(UnitCell<D> const & unitCell)
   {
      // Preconditions
      UTIL_CHECK(domain_().lattice() != UnitCell<D>::Null);
      UTIL_CHECK(domain_().lattice() == unitCell.lattice());
 
      // Set system unit cell, using UnitCell<D> assignment operator
      domain_().unitCell() = unitCell;
 
      // If necessary, make basis
      if (domain_().hasGroup() && !domain_().basis().isInitialized()) {
         domain_().makeBasis();
      }
 
      // Postconditions
      UTIL_CHECK(domain_().unitCell().isInitialized());
      UTIL_CHECK(domain_().unitCell().lattice() == domain_().lattice());
      if (domain_().hasGroup()) {
         UTIL_CHECK(domain_().basis().isInitialized());
         UTIL_CHECK(isAllocatedBasis_);
      }
      UTIL_CHECK(!domain_().waveList().hasKSq());
   }
 
   /*
   * Set parameters of the system unit cell.
   */
   template <int D, class T>
   void System<D,T>::setUnitCell(FSArray<double, 6> const & parameters)
   {
      // Precondition
      UTIL_CHECK(domain_().lattice() != UnitCell<D>::Null);
 
      // Set system unit cell
      if (domain_().unitCell().lattice() == UnitCell<D>::Null) {
         domain_().unitCell().set(domain_().lattice(), parameters);
      } else {
         domain_().unitCell().setParameters(parameters);
      }
 
      // If necessary, make basis
      if (domain_().hasGroup() && !domain_().basis().isInitialized()) {
         domain_().makeBasis();
      }
 
      // Postconditions
      UTIL_CHECK(domain_().unitCell().isInitialized());
      UTIL_CHECK(domain_().unitCell().lattice() == domain_().lattice());
      if (domain_().hasGroup()) {
         UTIL_CHECK(domain_().basis().isInitialized());
         UTIL_CHECK(isAllocatedBasis_);
      }
      UTIL_CHECK(!domain_().waveList().hasKSq());
   }
 
   /*
   * Notify System members that unit cell parameters have been modified.
   */
   template <int D, class T>
   void System<D,T>::clearUnitCellData()
   {
      clearCFields();
      mixture_().clearUnitCellData();
      domain_().waveList().clearUnitCellData();
      if (hasEnvironment()) {
         environment().reset();
      }
   }
 
   // Property Output
 
   /*
   * Write parameter file for SCFT, omitting any sweep block.
   */
   template <int D, class T>
   void System<D,T>::writeParamNoSweep(std::ostream& out) const
   {
      out << "System{" << std::endl;
      mixture().writeParam(out);
      interaction().writeParam(out);
      domain().writeParam(out);
      if (hasEnvironment()) {
         environment().writeParam(out);
      }
      if (hasIterator()) {
         iterator().writeParam(out);
      }
      out << "}" << std::endl;
   }
 
   // Timer Operations
 
   /*
   * Write timer values to output stream (computational cost).
   */
   template <int D, class T>
   void System<D,T>::writeTimers(std::ostream& out) const
   {
      if (hasIterator()) {
         iterator().outputTimers(Log::file());
         iterator().outputTimers(out);
      }
      if (hasSimulator()){
         simulator().outputTimers(Log::file());
         simulator().outputTimers(out);
      }
   }
 
   /*
   * Clear state of all timers.
   */
   template <int D, class T>
   void System<D,T>::clearTimers()
   {
      if (hasIterator()) {
         iterator().clearTimers();
      }
      if (hasSimulator()){
         simulator().clearTimers();
      }
   }
 
   // Private member functions
 
   /*
   * Allocate memory for fields in grid format.
   */
   template <int D, class T>
   void System<D,T>::allocateFieldsGrid()
   {
      // Preconditions
      UTIL_CHECK(hasMixture_);
      int nMonomer = mixture_().nMonomer();
      UTIL_CHECK(nMonomer > 0);
      UTIL_CHECK(domain_().mesh().size() > 0);
      UTIL_CHECK(!isAllocatedGrid_);
 
      // Alias for mesh dimensions
      IntVec<D> const & dimensions = domain_().mesh().dimensions();
 
      // Allocate w (chemical potential) fields
      w().setNMonomer(nMonomer);
      w().allocateRGrid(dimensions);
 
      // Allocate c (monomer concentration) fields
      c_().setNMonomer(nMonomer);
      c_().allocateRGrid(dimensions);
 
      // Set nMonomer for external field container
      h().setNMonomer(nMonomer);
 
      // If hasEnvironment(), allocate mask and h fields
      if (hasEnvironment()) {
         if (environment().generatesMask()) {
            mask().allocateRGrid(dimensions);
         }
         if (environment().generatesExternalFields()) {
            h().allocateRGrid(dimensions);
         }
      }
 
      isAllocatedGrid_ = true;
   }
 
   /*
   * Allocate memory for fields in basis format.
   */
   template <int D, class T>
   void System<D,T>::allocateFieldsBasis()
   {
      // Preconditions and constants
      UTIL_CHECK(hasMixture_);
      const int nMonomer = mixture_().nMonomer();
      UTIL_CHECK(nMonomer > 0);
      UTIL_CHECK(domain_().mesh().size() > 0);
      UTIL_CHECK(domain_().hasGroup());
      UTIL_CHECK(domain_().basis().isInitialized());
      const int nBasis = domain_().basis().nBasis();
      UTIL_CHECK(nBasis > 0);
      UTIL_CHECK(isAllocatedGrid_);
      UTIL_CHECK(!isAllocatedBasis_);
 
      // Allocate basis fields in w and c field containers
      w().allocateBasis(nBasis);
      c_().allocateBasis(nBasis);
 
      // If hasEnvironment(), allocate mask and h fields
      if (hasEnvironment()) {
         if (environment().generatesMask()) {
            mask().allocateBasis(nBasis);
         }
         if (environment().generatesExternalFields()) {
            h().allocateBasis(nBasis);
         }
      }
 
      isAllocatedBasis_ = true;
   }
 
   /*
   * Read a filename string and echo to log file (used in readCommands).
   */
   template <int D, class T>
   void System<D,T>::readEcho(std::istream& in, std::string& string) const
   {
      in >> string;
      if (in.fail()) {
          UTIL_THROW("Unable to read a string parameter.");
      }
      Log::file() << " " << Str(string, 20) << std::endl;
   }
 
   /*
   * Read floating point number, echo to log file (used in readCommands).
   */
   template <int D, class T>
   void System<D,T>::readEcho(std::istream& in, double& value) const
   {
      in >> value;
      if (in.fail()) {
          UTIL_THROW("Unable to read floating point parameter.");
      }
      Log::file() << " " << Dbl(value, 20) << std::endl;
   }
 
} // namespace Rp
} // namespace Pscf
#endif