ElecInfo.h

/*-------------------------------------------------------------------
Copyright 2011 Ravishankar Sundararaman, Kendra Letchworth Weaver
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_ELECINFO_H
#define JDFTX_ELECTRONIC_ELECINFO_H

#include <electronic/common.h>
#include <core/vector3.h>

enum SpinType {
        SpinNone, 
        SpinZ, 
        SpinVector, 
        SpinOrbit 
};

class QuantumNumber
{
public:
        vector3<> k; 
        int spin;  
        double weight; 
        
        QuantumNumber() : spin(0) {}
        int index() const { return spin<0 ? 1 : 0; } 
};

inline vector3<> getCoord(const QuantumNumber& qnum) { return qnum.k; }

class ElecInfo
{
public:
        int nBands, nStates; 
        int nDensities, spinWeight, qWeightSum; 
        int qStart, qStop; 
        bool isMine(int q) const { return qDivision.isMine(q); } 
        int whose(int q) const { return qDivision.whose(q); } 
        int qStartOther(int iProc) const { return qDivision.start(iProc); } 
        int qStopOther(int iProc) const { return qDivision.stop(iProc); } 
        
        SpinType spinType; 
        double nElectrons; 
        std::vector<QuantumNumber> qnums; 
        
        bool isNoncollinear() const { return spinType==SpinVector || spinType==SpinOrbit; }
        int spinorLength() const { return isNoncollinear() ? 2 : 1; }
        
        enum FillingsUpdate
        {       ConstantFillings, 
                FermiFillingsMix, 
                FermiFillingsAux, 
                MaximumOverlapMethod 
        }
        fillingsUpdate;
        
        double kT; 
        double mu; 
        
        int mixInterval; 
        bool subspaceRotation; 
        
        bool hasU; 

        std::vector<std::tuple<int, int, double>> customFillings;
        string initialFillingsFilename; 
        
        ElecInfo();
        void setup(const Everything &e, std::vector<diagMatrix>& F, Energies& ener); 
        void printFillings(FILE* fp) const;
        void printFermi(const char* suffix, const double* muOverride=0) const; //Fermi fillings report (compute mu from eigenvalues in eVars if muOverride not provided)
        void mixFillings(std::vector<diagMatrix>& F, Energies& ener); 
        void updateFillingsEnergies(const std::vector<diagMatrix>& F, Energies&) const; 

        //Fermi function utilities:
        inline double muEff(double mu, double Bz, int q) const { return mu + Bz*qnums[q].spin; }
        double fermi(double mu, double eps) const { return 0.5*(1.-tanh(betaBy2*(eps-mu))); } 
        double fermiPrime(double mu, double eps) const { return -0.5*betaBy2/pow(cosh(betaBy2*(eps-mu)), 2); } 
        diagMatrix fermi(double mu, const diagMatrix& eps) const; 
        diagMatrix fermiPrime(double mu, const diagMatrix& eps) const; 
        
        matrix fermiGrad(double mu, const diagMatrix& eps, const matrix& gradF) const;

        double nElectronsFermi(double mu, const std::vector<diagMatrix>& eps, double& Bz) const; 
        
        double findMu(const std::vector<diagMatrix>& eps, double nElectrons, double& Bz) const; 
        
        void kpointsPrint(FILE* fp, bool printSpin=false) const; 
        void kpointPrint(FILE* fp, int q, bool printSpin=false) const; 
        
        int findHOMO(int q) const; 

        //Parallel I/O utilities for diagMatrix/matrix array (one-per-kpoint, with nBands rows and columns unless overridden):
        void read(std::vector<diagMatrix>&, const char *fname, int nRowsOverride=0) const;
        void read(std::vector<matrix>&, const char *fname, int nRowsOverride=0, int nColsOverride=0) const;
        void write(const std::vector<diagMatrix>&, const char *fname, int nRowsOverride=0) const;
        void write(const std::vector<matrix>&, const char *fname, int nRowsOverride=0, int nColsOverride=0) const;
        void appendWrite(const std::vector<diagMatrix>&, const char *fname, int nRowsOverride=0) const;

private:
        const Everything* e;
        double betaBy2; 
        TaskDivision qDivision; 
        
        //Initial fillings:
        int nBandsOld; 
        double Qinitial, Minitial; 
        bool Mconstrain; 
        friend struct CommandElecInitialFillings;
        friend struct CommandElecInitialCharge;
        friend struct CommandElecInitialMagnetization;
        friend struct CommandInitialState;
        friend class ElecVars;
        friend struct LCAOminimizer;
        
        double magnetizationFermi(double mu, double Bz, const std::vector<diagMatrix>& eps, double& nElectrons) const; 
        
        //Fillings mix / mu-controller properties:
        double fillingMixFraction; 
        double dnPrev, muMeasuredPrev; 
        double Cmeasured, Cweight; 
        double dnMixFraction; 
        friend struct CommandElecFermiFillings;
        friend struct CommandTargetMu;
        
        //k-points:
        vector3<int> kfold; 
        friend struct CommandKpointFolding;
        friend class Everything;
        friend class Phonon;
        void kpointsFold(); 
        void kpointsReduce(); 
};

#endif // JDFTX_ELECTRONIC_ELECINFO_H