ElecVars.h

/*-------------------------------------------------------------------
Copyright 2011 Ravishankar Sundararaman
Copyright 1996-2003 Sohrab Ismail-Beigi

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_ELECTRONIC_ELECVARS_H
#define JDFTX_ELECTRONIC_ELECVARS_H

#include <electronic/common.h>
#include <fluid/FluidSolverParams.h>
#include <core/ScalarFieldArray.h>
#include <string>
#include <memory>

class ElecVars
{
public:
        //Independent variables:
        std::vector<ColumnBundle> Y; 
        std::vector<matrix> B; 
        double subspaceRotationFactor; 
        
        //Derived quantities:
        std::vector<matrix> B_evecs; 
        std::vector<diagMatrix> B_eigs; 

        std::vector<matrix> Hsub; 
        std::vector<matrix> Hsub_evecs; 
        std::vector<diagMatrix> Hsub_eigs; 
        
        std::vector<ColumnBundle> C; 
        std::vector<diagMatrix> F;  
        
        std::vector< std::vector<matrix> > VdagC; //cached pseudopotential projections (by state and then species)
        std::vector<matrix> grad_CdagOC; 
        
        //Densities and potentials:
        ScalarFieldArray n; 
        ScalarFieldArray get_nXC() const; 
        ScalarField get_nTot() const { return n.size()==1 ? n[0] : n[0]+n[1]; } 
        
        ScalarFieldArray tau; 
        
        ScalarFieldTilde d_fluid; 
        ScalarFieldTilde V_cavity; 
        
        ScalarFieldArray Vscloc; 
        ScalarFieldArray Vtau; 
        
        std::vector<matrix> rhoAtom, U_rhoAtom; 
        
        //External interactions:
        ScalarFieldArray Vexternal; 
        ScalarFieldTilde rhoExternal; 
        bool rhoExternalSelfEnergy; 
        struct BoxPotential
        {       vector3<> min, max;
                double Vin, Vout, convolve_radius;
        };
        std::vector<BoxPotential> boxPot; 

        //Fluid properties
        FluidSolverParams fluidParams;
        std::shared_ptr<struct FluidSolver> fluidSolver;
        string fluidInitialStateFilename;

        //Wavefunction initialization:
        string wfnsFilename; 
        std::shared_ptr<struct ColumnBundleReadConversion> readConversion; 
        bool isRandom; 
        bool initLCAO; 

        string eigsFilename; 
        
        //Auxiliary hamiltonian initialization
        string HauxFilename; 
        bool HauxInitialized; 
        
        double overlapCondition; 
        
        string nFilenamePattern; 
        string VFilenamePattern; 
        
        ElecVars();
        void setup(const Everything &everything);

        void EdensityAndVscloc(Energies& ener, const ExCorr* alternateExCorr=0);
        
        double elecEnergyAndGrad(Energies& ener, ElecGradient* grad=0, ElecGradient* Kgrad=0, bool calc_Hsub=false); 
        
        void setEigenvectors(int q=-1); 
        
        int nOccupiedBands(int q) const; 

        ScalarFieldArray KEdensity() const;
        
        ScalarFieldArray calcDensity() const;
        
        void orthonormalize(int q);
        
        double applyHamiltonian(int q, const diagMatrix& Fq, ColumnBundle& HCq, Energies& ener, bool need_Hsub=false);
        
        void orthonormalizeGrad(int q, const diagMatrix& Fq, const ColumnBundle& HCq, ColumnBundle& gradYq, double KErollover=1., ColumnBundle* KgradYq=0, matrix* gradBq=0, matrix* KgradBq=0);

        double bandEnergyAndGrad(int q, Energies& ener, ColumnBundle* grad=0, ColumnBundle* Kgrad=0);
        
private:
        const Everything* e;
        
        std::vector<string> VexternalFilename; 
        friend struct CommandVexternal;
        friend class InverseKohnSham; 
        
        string rhoExternalFilename; 
        friend struct CommandRhoExternal;

        int lcaoIter; 
        double lcaoTol; 
        int LCAO(); 
        friend struct CommandWavefunction;
        friend struct CommandLcaoParams;
        friend class Dump;
        
        std::vector<matrix> U; // U[q] = Y[q]^O(Y[q])
        std::vector<matrix> U_evecs; // eigenvectors of U[q] in columns
        std::vector<diagMatrix> U_eigs; // eigenvalues of U[q]
        std::vector<matrix> Umhalf; // Uhmalf[q] = invsqrt(U[q])
        std::vector<matrix> V; // V=cis(B) or dagger(B_evecs) for subspace rotations / aux hamiltonian respectively
};
#endif // JDFTX_ELECTRONIC_ELECVARS_H