cpc/fts/Simulator.h

#ifndef CPC_SIMULATOR_H
#define CPC_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 <prdc/cpu/CField.h>               // member (template arg)
#include <util/random/Random.h>            // member
#include <util/containers/DArray.h>        // member (template)
#include <util/containers/DMatrix.h>       // member (template)
 
 
namespace Pscf {
namespace Cpc {
 
   // Forward declarations 
   template <int D> class System;
   template <int D> class Step;
   // template <int D> class TrajectoryReader;
 
   using namespace Util;
   using namespace Prdc;
   using namespace Prdc::Cpu;

   template <int D>
   class Simulator : public ParamComposite
   {
 
   public:
 
      // Suppress automatically generated functions
      Simulator() = delete;
      Simulator(Simulator<D>& ) = delete;

      Simulator(System<D>& system);

      ~Simulator();

      void allocate();

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


      virtual void simulate(int nStep);
 
      #if 0
      virtual void analyze(int min, int max,
                           std::string classname,
                           std::string filename);
      #endif

      void clearData();


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

      virtual void clearTimers();

      long iStep();
      
      long iTotalStep();


      void analyzeInteraction();

      DMatrix<double> const & u() const;

      DArray<double> const & evals() const;

      double eval(int i) const;

      bool isPositiveEval(int i) const;

      DMatrix<double> const & evecs() const;

      double evecs(int a, int i) const;
 
      #if 0
      DArray<double> const & sc() const;

      double sc(int a) const;
      #endif


      void computeHamiltonian();

      std::complex<double> hamiltonian() const;

      std::complex<double> idealHamiltonian() const;

      std::complex<double> fieldHamiltonian() const;

      bool hasHamiltonian() const;


      void computeWc();

      DArray< CField<D> > const & wc() const;

      CField<D> const & wc(int a) const;

      bool hasWc() const;


      void computeCc();

      DArray< CField<D> > const & cc() const;

      CField<D> const & cc(int a) const;

      bool hasCc() const;


      void computeDc();

      DArray< CField<D> > const & dc() const;

      CField<D> const & dc(int i) const;

      bool hasDc() const;
      

      System<D>& system();

      Random& random();
      
      bool hasStep() const;

      Step<D>& step();
 
      #if 0
      AnalyzerManager<D>& analyzerManager();

      Factory<TrajectoryReader<D>>& trajectoryReaderFactory();
      #endif

 
   private:
 
      #if 0
      AnalyzerManager<D> analyzerManager_;
      #endif

      Random random_;
 
      // Field components in eigenvector basis

      DArray< CField<D> > wc_;

      DArray< CField<D> > cc_;

      DArray< CField<D> > dc_;
     
      // Interaction matrix and eigen-properties

      DMatrix<double> u_;

      DMatrix<double> evecs_;

      DArray<double>  evals_;

      DArray<bool> isPositiveEval_;
 
      // Hamiltonian and components
 
      std::complex<double> hamiltonian_;

      std::complex<double> idealHamiltonian_;

      std::complex<double> fieldHamiltonian_;
 
      // Pointers to associated and child objects

      System<D>* systemPtr_;

      Step<D>* stepPtr_;

      Factory< Step<D> >* stepFactoryPtr_;
 
      #if 0
      Factory< TrajectoryReader<D> >* trajectoryReaderFactoryPtr_;
      #endif
 
      // Counters and random seed

      long iStep_;
      
      long iTotalStep_;

      long seed_;
 
      // Boolean status indicators

      bool hasHamiltonian_;

      bool hasWc_;

      bool hasCc_;

      bool hasDc_;

      bool isAllocated_;
 
      // Private member functions

      void readRandomSeed(std::istream& in);

      void setup(int nStep);
 
   };
 
   // Inline functions
 
   // Management of owned and associated objects
 
   // Get the parent System by reference.
   template <int D>
   inline System<D>& Simulator<D>::system()
   {
      UTIL_ASSERT(systemPtr_);  
      return *systemPtr_; 
   }
 
   // Get the random number generator by reference.
   template <int D>
   inline Random& Simulator<D>::random()
   {  return random_; }
 
   // Get the Brownian dynamics stepper.
   template <int D>
   inline bool Simulator<D>::hasStep() const
   {  return (bool)stepPtr_; }
 
   // Get the Brownian dynamics stepper.
   template <int D>
   inline Step<D>& Simulator<D>::step()
   {
      UTIL_CHECK(hasStep());  
      return *stepPtr_; 
   }
 
   #if 0
   // Get the analyzer manager.
   template <int D>
   inline AnalyzerManager<D>& Simulator<D>::analyzerManager()
   {  return analyzerManager_; }
 
   // Get the TrajectoryReaderfactory
   template <int D>
   inline
   Factory<TrajectoryReader<D> >& Simulator<D>::trajectoryReaderFactory()
   {
      UTIL_ASSERT(trajectoryReaderFactoryPtr_);
      return *trajectoryReaderFactoryPtr_;
   }
   #endif
 
   // Interaction Matrix
 
   // Return the interaction matrix U.
   template <int D>
   inline DMatrix<double> const & Simulator<D>::u() const
   {  return u_; }
 
   // Return an array of eigenvalues of the interaction matrix U.
   template <int D>
   inline DArray<double> const & Simulator<D>::evals() const
   {  return evals_; }
 
   // Return a single eigenvalue of the interaction matrix U.
   template <int D>
   inline double Simulator<D>::eval(int a) const
   {  return evals_[a]; }
 
   // Return a matrix of eigenvectors of the interaction matrix U.
   template <int D>
   inline DMatrix<double> const & Simulator<D>::evecs() const
   {  return evecs_; }
 
   // Return an element of an eigenvector of the interaction matrix U.
   template <int D>
   inline double Simulator<D>::evecs(int a, int i) const
   {  return evecs_(a, i); }
 
   // Hamiltonian and its derivatives
 
   // Get the precomputed Hamiltonian.
   template <int D>
   inline std::complex<double> Simulator<D>::hamiltonian() const
   {
      UTIL_CHECK(hasHamiltonian_);
      return hamiltonian_;
   }
 
   // Get the ideal gas component of the precomputed Hamiltonian.
   template <int D>
   inline std::complex<double> Simulator<D>::idealHamiltonian() const
   {
      UTIL_CHECK(hasHamiltonian_);
      return idealHamiltonian_;
   }
 
   // Get the W field component of the precomputed Hamiltonian.
   template <int D>
   inline std::complex<double> Simulator<D>::fieldHamiltonian() const
   {
      UTIL_CHECK(hasHamiltonian_);
      return fieldHamiltonian_;
   }
 
   // Has the Hamiltonian been computed for the current w fields ?
   template <int D>
   inline bool Simulator<D>::hasHamiltonian() const
   {  return hasHamiltonian_; }
 
   // Fields
 
   // Return all eigencomponents of the w fields.
   template <int D>
   inline DArray< CField<D> > const & Simulator<D>::wc() const
   {  return wc_; }
 
   // Return a single eigenvector component of the w fields.
   template <int D>
   inline CField<D> const & Simulator<D>::wc(int a) const
   {  return wc_[a]; }
 
   // Have eigenvector components of current w fields been computed?
   template <int D>
   inline bool Simulator<D>::hasWc() const
   {  return hasWc_; }
 
   // Return all eigenvector components of the current c fields.
   template <int D>
   inline DArray< CField<D> > const & Simulator<D>::cc() const
   {  return cc_; }
 
   // Return a single eigenvector component of the current c fields.
   template <int D>
   inline CField<D> const & Simulator<D>::cc(int a) const
   {  return cc_[a]; }
 
   // Have eigenvector components of current c fields been computed?
   template <int D>
   inline bool Simulator<D>::hasCc() const
   {  return hasCc_; }
 
   // Return all eigenvector components of the current d fields.
   template <int D>
   inline DArray< CField<D> > const & Simulator<D>::dc() const
   {  return dc_; }
 
   // Return a single eigenvector component of the current d fields.
   template <int D>
   inline CField<D> const & Simulator<D>::dc(int a) const
   {  return dc_[a]; }
 
   // Have eigenvector components of current d fields been computed?
   template <int D>
   inline bool Simulator<D>::hasDc() const
   {  return hasDc_; }
 
   // Return the current converged simulation step index.
   template <int D>
   inline long Simulator<D>::iStep()
   {  return iStep_; }
   
   // Return the current simulation step index.
   template <int D>
   inline long Simulator<D>::iTotalStep()
   {  return iTotalStep_; }
 
   // Explicit instantiation declarations
   extern template class Simulator<1>;
   extern template class Simulator<2>;
   extern template class Simulator<3>;
 
} // namespace Cpc
} // namespace Pscf
#endif