FluidMixture.h

/*-------------------------------------------------------------------
Copyright 2011 Ravishankar Sundararaman, Kendra Letchworth Weaver

This file is part of JDFTx.

JDFTx is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

JDFTx is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with JDFTx.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------*/

#ifndef JDFTX_FLUID_FLUIDMIXTURE_H
#define JDFTX_FLUID_FLUIDMIXTURE_H

#include <fluid/FluidComponent.h>
#include <fluid/Fmix.h>
#include <core/Units.h>
#include <core/Minimize.h>
#include <core/EnergyComponents.h>

class FluidMixture : public Minimizable<ScalarFieldArray>
{
public:
        const GridInfo& gInfo;
        const double T; 
        bool verboseLog; 
        vector3<> Eexternal; 

        FluidMixture(const GridInfo& gInfo, const double T=298*Kelvin);
        virtual ~FluidMixture();

        void initialize(double P=1.01325*Bar, double epsBulkOverride=0., double epsInfOverride=0.);

        double getBoilingPressure(double NliqGuess, double NvapGuess, std::vector<double>* Nvap=0) const;
        
        unsigned get_nIndep() const { return nIndepIdgas + (polarizable ? 3 : 0); }  
        unsigned get_nDensities() const { return nDensities; } 

        const std::vector<const FluidComponent*>& getComponents() const; 
        ScalarFieldArray state;

        ScalarFieldTilde rhoExternal;
        bool useMFKernel; 
        double Qtol; 

        void initState(double scale = 0.0, double Elo=-DBL_MAX, double Ehi=+DBL_MAX);

        void loadState(const char* filename);

        void saveState(const char* filename) const;

        struct Outputs
        {       
                ScalarFieldArray* N; 
                vector3<>* electricP; 
                ScalarFieldTilde* Phi_rhoExternal; 
                ScalarFieldArray* psiEff; 
                EnergyComponents* Phi; 
                
                Outputs( ScalarFieldArray* N=0, vector3<>* electricP=0,
                        ScalarFieldTilde* Phi_rhoExternal=0, ScalarFieldArray* psiEff=0, EnergyComponents* Phi=0);
        };

        double operator()(const ScalarFieldArray& indep, ScalarFieldArray& Phi_indep, Outputs outputs=Outputs()) const;

        double getFreeEnergy(Outputs outputs=Outputs()) const;

        void step(const ScalarFieldArray& dir, double alpha);

        double compute(ScalarFieldArray* grad, ScalarFieldArray* Kgrad);
        
        double sync(double x) const; 
        
private:
        unsigned nIndepIdgas; 
        unsigned nDensities; 
        bool polarizable;  
        double p;
        double Crot, Cpol; 
        
        std::vector<const FluidComponent*> component; 
        std::vector<const Fmix*> fmixArr; 

        std::vector<double> Kindep; 
        double Keps; 
        
        void addComponent(FluidComponent*); 
        friend struct FluidComponent;

        void addFmix(const Fmix* fmix); 
        friend class Fmix;

        double computeUniformEx(const std::vector< double >& Nmol, std::vector< double >& Phi_Nmol) const;

        double compute_p(double Ntot) const;
        
        friend struct BoilingPressureSolver;
        static ScalarFieldTilde coulomb(const ScalarFieldTilde& rho) { return (-4*M_PI) * Linv(O(rho)); }
};

#endif // JDFTX_FLUID_FLUIDMIXTURE_H