rp/fts/simulator/Simulator.h

#ifndef RP_SIMULATOR_H
#define RP_SIMULATOR_H
 
/*
* 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 <util/param/ParamComposite.h>     // base class
#include <util/containers/DArray.h>        // member (template)
#include <util/containers/DMatrix.h>       // member (template)
#include <iostream>
 
// Forward declaration
namespace Util { class Random; }
 
namespace Pscf {
namespace Rp {
 
   using namespace Util;

   template <int D, class T>
   class Simulator : public ParamComposite
   {
 
   public:

      using RFieldT = typename T::RField;
 
      // Protected constructor and destructor (see below).


      virtual void readParameters(std::istream &in);

      void allocate();


      virtual void simulate(int nStep);

      virtual void analyze(int min, int max,
                           std::string classname,
                           std::string filename);

      void clearData();


      virtual void outputTimers(std::ostream& out) const;

      virtual void outputMdeCounter(std::ostream& out) const;

      virtual void clearTimers();

      long iStep();

      long iTotalStep();


      void analyzeChi();

      DArray<double> const & chiEvals() const;

      double chiEval(int a) const;

      DMatrix<double> const & chiEvecs() const;

      double chiEvecs(int a, int i) const;

      DArray<double> const & sc() const;

      double sc(int a) const;


      void computeHamiltonian();

      double hamiltonian() const;

      double idealHamiltonian() const;

      double fieldHamiltonian() const;

      double perturbationHamiltonian() const;

      bool hasHamiltonian() const;


      void computeWc();

      DArray<RFieldT> const & wc() const;

      RFieldT const & wc(int a) const;

      bool hasWc() const;


      void computeCc();

      DArray<RFieldT> const & cc() const;

      RFieldT const & cc(int a) const;

      bool hasCc() const;


      void computeDc();

      DArray<RFieldT> const & dc() const;

      RFieldT const & dc(int i) const;

      bool hasDc() const;


      void saveState();

      void restoreState();

      void clearState();


      typename T::System& system();

      Random& random();

      typename T::VecRandom& vecRandom();

      bool hasCompressor() const;

      typename T::Compressor const & compressor() const;

      typename T::Compressor& compressor();

      bool hasPerturbation() const;

      typename T::Perturbation const & perturbation() const;

      typename T::Perturbation& perturbation();

      bool hasRamp() const;

      typename T::Ramp const & ramp() const;

      typename T::Ramp& ramp();

 
   protected:
 
      // Protected member functions

      Simulator(typename T::System& system, 
                typename T::Simulator& simulator);

      ~Simulator();
 
      // Prohibit copying and assignment.
      Simulator(Simulator<D,T> const &) = delete;
      Simulator<D,T>& operator = (Simulator<D,T> const &) = delete;

      void readRandomSeed(std::istream& in);

      virtual void initializeVecRandom()
      {};

      typename T::CompressorFactory& compressorFactory();

      void readCompressor(std::istream& in, bool& isEnd);

      typename T::PerturbationFactory& perturbationFactory();

      void readPerturbation(std::istream& in, bool& isEnd);

      void setPerturbation(typename T::Perturbation* ptr);

      typename T::RampFactory& rampFactory();

      void readRamp(std::istream& in, bool& isEnd);

      void setRamp(typename T::Ramp* ptr);

      typename T::SimState& state();
 
      // Protected data members

      DArray<RFieldT> wc_;

      DArray<RFieldT> cc_;

      DArray<RFieldT> dc_;

      mutable typename T::SimState state_;

      double hamiltonian_;

      double idealHamiltonian_;

      double fieldHamiltonian_;

      double perturbationHamiltonian_;

      long iStep_;

      long iTotalStep_;

      long seed_;
 
      // Boolean status flags

      bool hasHamiltonian_;

      bool hasWc_;

      bool hasCc_;

      bool hasDc_;
 
   private:

      DMatrix<double> chiP_;

      DMatrix<double> chiEvecs_;

      DArray<double>  chiEvals_;

      DArray<double>  sc_;

      mutable RFieldT tmpField_;
 
      // Pointers to associated objects

      typename T::System* systemPtr_;

      typename T::Simulator* simulatorPtr_;

      Random* randomPtr_;

      typename T::VecRandom* vecRandomPtr_;

      typename T::CompressorFactory* compressorFactoryPtr_;

      typename T::Compressor* compressorPtr_;

      typename T::PerturbationFactory* perturbationFactoryPtr_;

      typename T::Perturbation* perturbationPtr_;

      typename T::RampFactory* rampFactoryPtr_;

      typename T::Ramp* rampPtr_;

      bool isAllocated_;
 
   };
 
   // Inline functions
 
   // Access to associated objects via pointers
 
   // Get the parent System by reference.
   template <int D, class T> inline 
   typename T::System& Simulator<D,T>::system()
   {
      UTIL_ASSERT(systemPtr_);
      return *systemPtr_;
   }
 
   // Get the scalar random number generator by reference.
   template <int D, class T> inline 
   Random& Simulator<D,T>::random()
   {
      UTIL_ASSERT(randomPtr_);
      return *randomPtr_;
   }
 
   // Get the vector random number generator by reference.
   template <int D, class T> inline 
   typename T::VecRandom& Simulator<D,T>::vecRandom()
   {
      UTIL_ASSERT(vecRandomPtr_);
      return *vecRandomPtr_;
   }
 
   // Does this Simulator have a Compressor?
   template <int D, class T> inline 
   bool Simulator<D,T>::hasCompressor() const
   {  return (bool)compressorPtr_; }
 
   // Get the Compressor by non-const reference.
   template <int D, class T> inline 
   typename T::Compressor& Simulator<D,T>::compressor()
   {
      UTIL_CHECK(compressorPtr_);
      return *compressorPtr_;
   }
 
   // Get the Compressor by const reference.
   template <int D, class T> inline 
   typename T::Compressor const & Simulator<D,T>::compressor() const
   {
      UTIL_CHECK(compressorPtr_);
      return *compressorPtr_;
   }
 
   // Does this Simulator have an associated Perturbation?
   template <int D, class T> inline 
   bool Simulator<D,T>::hasPerturbation() const
   {  return (bool)perturbationPtr_; }
 
   // Get a Perturbation by const reference.
   template <int D, class T> inline 
   typename T::Perturbation const & Simulator<D,T>::perturbation() const
   {
      UTIL_CHECK(perturbationPtr_);
      return *perturbationPtr_;
   }
 
   // Get a Perturbation by non-const reference.
   template <int D, class T> inline 
   typename T::Perturbation& Simulator<D,T>::perturbation()
   {
      UTIL_CHECK(perturbationPtr_);
      return *perturbationPtr_;
   }
 
   // Does this Simulator have an associated Ramp?
   template <int D, class T> inline
   bool Simulator<D,T>::hasRamp() const
   {  return (bool)rampPtr_; }
 
   // Get a Ramp by const reference.
   template <int D, class T> inline 
   typename T::Ramp const & Simulator<D,T>::ramp() const
   {
      UTIL_CHECK(rampPtr_);
      return *rampPtr_;
   }
 
   // Get a Ramp by non-const reference.
   template <int D, class T> inline 
   typename T::Ramp& Simulator<D,T>::ramp()
   {
      UTIL_CHECK(rampPtr_);
      return *rampPtr_;
   }
 
   // Get the stored internal state by reference.
   template <int D, class T> inline 
   typename T::SimState& Simulator<D,T>::state()
   {  return state_; }
 
   // Projected Chi Matrix
 
   // Return an array of eigenvalues of the projected chi matrix.
   template <int D, class T> inline 
   DArray<double> const & Simulator<D,T>::chiEvals() const
   {  return chiEvals_; }
 
   // Return a single eigenvalue of the projected chi matrix.
   template <int D, class T> inline 
   double Simulator<D,T>::chiEval(int a) const
   {  return chiEvals_[a]; }
 
   // Return a matrix of eigenvectors of the projected chi matrix.
   template <int D, class T> inline 
   DMatrix<double> const & Simulator<D,T>::chiEvecs() const
   {  return chiEvecs_; }
 
   // Return an element of an eigenvector of the projected chi matrix.
   template <int D, class T> inline 
   double Simulator<D,T>::chiEvecs(int a, int i) const
   {  return chiEvecs_(a, i); }
 
   // Return an array of values of vector S.
   template <int D, class T> inline 
   DArray<double> const & Simulator<D,T>::sc() const
   {  return sc_; }
 
   // Return one component of vector S.
   template <int D, class T> inline 
   double Simulator<D,T>::sc(int a) const
   {  return sc_[a]; }
 
   // Hamiltonian and its components
 
   // Has the Hamiltonian been computed for the current w fields ?
   template <int D, class T> inline 
   bool Simulator<D,T>::hasHamiltonian() const
   {  return hasHamiltonian_; }
 
   // Get the precomputed total Hamiltonian.
   template <int D, class T> inline 
   double Simulator<D,T>::hamiltonian() const
   {
      UTIL_CHECK(hasHamiltonian_);
      return hamiltonian_;
   }
 
   // Get the ideal gas component of the precomputed Hamiltonian.
   template <int D, class T> inline 
   double Simulator<D,T>::idealHamiltonian() const
   {
      UTIL_CHECK(hasHamiltonian_);
      return idealHamiltonian_;
   }
 
   // Get the harmonic field component of the precomputed Hamiltonian.
   template <int D, class T> inline 
   double Simulator<D,T>::fieldHamiltonian() const
   {
      UTIL_CHECK(hasHamiltonian_);
      return fieldHamiltonian_;
   }
 
   // Get the perturbation component of the precomputed Hamiltonian.
   template <int D, class T> inline 
   double Simulator<D,T>::perturbationHamiltonian() const
   {
      UTIL_CHECK(hasHamiltonian_);
      return perturbationHamiltonian_;
   }
 
   // Fields
 
   // Have eigenvector components of the current w fields been computed?
   template <int D, class T> inline 
   bool Simulator<D,T>::hasWc() const
   {  return hasWc_; }
 
   // Return all eigencomponents of the w fields.
   template <int D, class T> inline 
   DArray<typename T::RField> const & Simulator<D,T>::wc() const
   {  return wc_; }
 
   // Return a single eigenvector component of the w fields.
   template <int D, class T> inline 
   typename T::RField const & Simulator<D,T>::wc(int a) const
   {  return wc_[a]; }
 
   // Have eigenvector components of the current c fields been computed?
   template <int D, class T> inline 
   bool Simulator<D,T>::hasCc() const
   {  return hasCc_; }
 
   // Return all eigenvector components of the current c fields.
   template <int D, class T> inline 
   DArray<typename T::RField> const & Simulator<D,T>::cc() const
   {  return cc_; }
 
   // Return a single eigenvector component of the current c fields.
   template <int D, class T> inline 
   typename T::RField const & Simulator<D,T>::cc(int a) const
   {  return cc_[a]; }
 
   // Have eigenvector components of the current d fields been computed?
   template <int D, class T>
   inline bool Simulator<D,T>::hasDc() const
   {  return hasDc_; }
 
   // Return all eigenvector components of the current d fields.
   template <int D, class T> inline 
   DArray<typename T::RField> const & Simulator<D,T>::dc() const
   {  return dc_; }
 
   // Return a single eigenvector component of the current d fields.
   template <int D, class T> inline 
   typename T::RField const & Simulator<D,T>::dc(int a) const
   {  return dc_[a]; }
 
   // Return the current converged simulation step index.
   template <int D, class T>
   inline long Simulator<D,T>::iStep()
   {  return iStep_; }
 
   // Return the current total simulation step index.
   template <int D, class T>
   inline long Simulator<D,T>::iTotalStep()
   {  return iTotalStep_; }
 
}
}
#endif