McNVTChemicalPotential.h

#ifndef MCMD_MC_NVT_CHEMICAL_POTENTIAL_H
#define MCMD_MC_NVT_CHEMICAL_POTENTIAL_H

/*
* Simpatico - Simulation Package for Polymeric and Molecular Liquids
*
* Copyright 2010 - 2017, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/

#include <mcMd/analyzers/SystemAnalyzer.h>  // base class template
#include <mcMd/mcSimulation/McSystem.h>     // base template parameter
#include <simp/ensembles/EnergyEnsemble.h>  // inline function
#include <util/accumulators/Average.h>      // member
#include <util/accumulators/Distribution.h> // member
#include <util/random/Random.h>             // member
#include <util/archives/Serializable.h>     // typedef

#include <cstdio>

namespace McMd
{

   using namespace Util;
   using namespace Simp;

   class Atom;

   class McNVTChemicalPotential : public SystemAnalyzer<McSystem>
   {

   public:

      McNVTChemicalPotential(McSystem& system);

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

      virtual void setup();
   
      virtual void sample(long iStep);
   
      virtual void output();

      virtual void save(Serializable::OArchive& ar);

      virtual void load(Serializable::IArchive& ar);

      template <class Archive>
      void serialize(Archive& ar, const unsigned int version);

   private:

      McSystem& system();

      Simulation& simulation();

      Boundary& boundary();

      EnergyEnsemble& energyEnsemble();

      Random& random();

      double boltzmann(double energy);

      void addEndAtom(Molecule* molPtr, int atomId, double &rosenbluth,
                      double &energy);

      void uniformCone(const double length, const double angle, 
                       const Vector n, Vector &p);

      static const int MaxTrial_ = 200;

      McSystem  *systemPtr_;

      Simulation  *simulationPtr_;

      Boundary  *boundaryPtr_;

      EnergyEnsemble  *energyEnsemblePtr_;

      Random  *randomPtr_;
   
      std::ofstream outputFile_;

      Average  accumulator_;

      std::ofstream EoutputFile_;

      Distribution  Eaccumulator_;

      double Emin_;

      double Emax_;

      double EnBin_;

      std::ofstream EcoutputFile_;

      Distribution  Ecaccumulator_;

      double Ecmin_;

      double Ecmax_;

      double EcnBin_;

      std::ofstream EmoutputFile_;

      Distribution  Emaccumulator_;

      double Emmin_;

      double Emmax_;

      double EmnBin_;

      std::ofstream BRoutputFile_;

      Distribution  BRaccumulator_;

      double BRmin_;

      double BRmax_;

      double BRnBin_;

      int  nTrial_; 

      int  speciesId_; 

      int  nMoleculeTrial_; 
   
      int nSamplePerBlock_;
   
      bool isInitialized_;

      //  friend class CbEndBaseTest;

   };

   // Inline methods

   /*
   * Get parent McSystem.
   */
   inline McSystem& McNVTChemicalPotential::system()
   {  return *systemPtr_; }

   /*
   * Get Simulation object of parent McSystem.
   */
   inline Simulation& McNVTChemicalPotential::simulation()
   {  return *simulationPtr_; }

   /*
   * Get Boundary object of parent McSystem.
   */
   inline Boundary& McNVTChemicalPotential::boundary()
   {  return *boundaryPtr_; }

   /*
   * Get EnergyEnsemble object of parent McSystem.
   */
   inline EnergyEnsemble& McNVTChemicalPotential::energyEnsemble()
   {  return *energyEnsemblePtr_; }

   /*
   * Get random object of parent McSystem.
   */
   inline Random& McNVTChemicalPotential::random()
   {  return *randomPtr_; }

   /*
   * Boltzmann weight associated with an energy difference.
   */
   inline double McNVTChemicalPotential::boltzmann(double energy)
   {  return exp(-energyEnsemblePtr_->beta()*energy); }

   /*
   * Serialize to/from an archive. 
   */
   template <class Archive>
   void McNVTChemicalPotential::serialize(Archive& ar, 
                                          const unsigned int version)
   {
      if (!isInitialized_) {
         UTIL_THROW("Error: Object not initialized.");
      }

      ar & accumulator_;
      ar & nSamplePerBlock_;
   }

}
#endif