pscf/floryHuggins/Mixture.h

#ifndef PSCF_FLORY_HUGGINS_MIXTURE_H
#define PSCF_FLORY_HUGGINS_MIXTURE_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 <pscf/chem/Monomer.h>            // Member template argument
#include <pscf/floryHuggins/Molecule.h>    // Member template argument
#include <util/containers/DArray.h>       // Member template
#include <util/containers/DMatrix.h>      // Member template
 
namespace Pscf {
   class Interaction;
   class LuSolver;
   class MixtureBase;
}
 
namespace Pscf {
namespace FloryHuggins {
 
   using namespace Util;

   class Mixture : public ParamComposite
   {
   public:

      Mixture();

      ~Mixture();



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

      void initialize(MixtureBase const& mixture);

      void setNMolecule(int nMolecule);
      
      void setNMonomer(int nMonomer);



      void setComposition(DArray<double> const & phi);

      void computeMu(Interaction const & interaction, double xi = 0.0);

      void computePhi(Interaction const & interaction, 
                      DArray<double> const & mu, 
                      DArray<double> const & phi, 
                      double&  xi);

      void computeFreeEnergy(Interaction const & interaction);


 
      Molecule& molecule(int id);

      double mu(int id) const;

      double phi(int id) const;

      double c(int id) const;

      double fHelmholtz() const;

      double pressure() const;

      int nMolecule() const;

      int nMonomer() const;


      void validate() const;
 
   private:

      DArray<Molecule> molecules_;

      DArray<double> mu_;

      DArray<double> phi_;

      DArray<double> c_;

      DArray<double> w_;

      DArray<double> residual_;

      DArray<double> dX_;

      DMatrix<double> dWdC_;

      DMatrix<double> dWdPhi_;

      DMatrix<double> jacobian_;

      DArray<double> phiOld_;

      double fHelmholtz_;

      double pressure_;

      LuSolver* solverPtr_;

      int nMolecule_;

      int nMonomer_;

      bool hasComposition_;

      void computeC();

      void adjustXi(DArray<double> const & mu, double& xi); 

      void computeResidual(DArray<double> const & mu,
                           double& error);
 
   };
 
   // Inline member functions
 
   inline Molecule& Mixture::molecule(int id)
   {  
      UTIL_ASSERT(id >= 0);  
      UTIL_ASSERT(id < nMolecule_);  
      return molecules_[id]; 
   }
 
   inline double Mixture::mu(int id) const
   {  
      UTIL_ASSERT(id >= 0);  
      UTIL_ASSERT(id < nMolecule_);  
      return mu_[id]; 
   }
 
   inline double Mixture::phi(int id) const
   {
      UTIL_ASSERT(id >= 0);  
      UTIL_ASSERT(id < nMolecule_);  
      return phi_[id]; 
   }
 
   inline double Mixture::c(int id) const
   {  
      UTIL_ASSERT(id >= 0);  
      UTIL_ASSERT(id < nMonomer_);  
      return c_[id]; 
   }
 
   inline double Mixture::fHelmholtz() const
   {  return fHelmholtz_; }
 
   inline double Mixture::pressure() const
   {  return pressure_; }
 
   inline int Mixture::nMolecule() const
   {  return nMolecule_; }
 
   inline int Mixture::nMonomer() const
   {  return nMonomer_; }
 
} // namespace FloryHuggins
} // namespace Pscf
#endif