SpeciesInfo.h

/*-------------------------------------------------------------------
Copyright 2011 Ravishankar Sundararaman, Deniz Gunceler
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_SPECIESINFO_H
#define JDFTX_ELECTRONIC_SPECIESINFO_H

#include <electronic/common.h>
#include <electronic/RadialFunction.h>
#include <electronic/VanDerWaals.h>
#include <electronic/matrix.h>
#include <core/ScalarFieldArray.h>
#include <core/vector3.h>
#include <core/string.h>


class SpeciesInfo
{
public:

        double Z; 
        int atomicNumber; 
        string name; 
        string potfilename, pulayfilename;
        bool fromWildcard; 
        
        std::vector<vector3<> > atpos; 
        std::vector<vector3<> > velocities; 
        #ifdef GPU_ENABLED
        vector3<> *atposGpu; 
        #endif
        void sync_atpos(); 
        vector3<>* atposPref; 
        
        double dE_dnG; 
        double mass; 
        double coreRadius; 
        double ZfullCore; 
        
        struct Constraint
        {       double moveScale; 
                vector3<> d; 
                enum ConstraintType {None, Linear, Planar, HyperPlane} type; 
                
                vector3<> operator()(const vector3<>& grad);
                bool isEquivalent(const Constraint& otherConstraint, const matrix3<>& transform) const;
                int getDimension() const;
                void print(FILE* fp, const Everything&) const;
        };
        std::vector<Constraint> constraints; 
        
        std::vector< vector3<> > initialMagneticMoments; 
        double initialOxidationState; 
        
        SpeciesInfo();
        ~SpeciesInfo();
        void setup(const Everything&);
        void print(FILE* fp) const;
        void populationAnalysis(const std::vector<matrix>& RhoAll) const; //print population analysis given the density matrix in the Lowdin basis
        bool isRelativistic() const { return psi2j.size(); } 
        
        enum PseudopotentialFormat
        {       Fhi, 
                Uspp, 
                UPF 
        };
        PseudopotentialFormat getPSPFormat(){return pspFormat;}

        std::shared_ptr<ColumnBundle> getV(const ColumnBundle& Cq, matrix* M=0) const; //get projectors with qnum and basis matching Cq  (optionally cached, and optionally retrieve full M repeated over atoms)

        double EnlAndGrad(const QuantumNumber& qnum, const diagMatrix& Fq, const matrix& VdagCq, matrix& HVdagCq) const;
        
        void augmentOverlap(const ColumnBundle& Cq, ColumnBundle& OCq, matrix* VdagCq=0) const;
        
        void augmentDensityInit();
        void augmentDensitySpherical(const QuantumNumber& qnum, const diagMatrix& Fq, const matrix& VdagCq);
        void augmentDensityGrid(ScalarFieldArray& n) const;
        
        void augmentDensityGridGrad(const ScalarFieldArray& E_n, std::vector<vector3<> >* forces=0);
        void augmentDensitySphericalGrad(const QuantumNumber& qnum, const diagMatrix& Fq, const matrix& VdagCq, matrix& HVdagCq) const;
        
        //DFT+U functions: handle IonInfo::rhoAtom_*() for this species
        //The rhoAtom pointers point to the start of those relevant to this species (and ends at that pointer + rhoAtom_nMatrices())
        size_t rhoAtom_nMatrices() const;
        void rhoAtom_initZero(matrix* rhoAtomPtr) const;
        void rhoAtom_calc(const std::vector<diagMatrix>& F, const std::vector<ColumnBundle>& C, matrix* rhoAtomPtr) const;
        double rhoAtom_computeU(const matrix* rhoAtomPtr, matrix* U_rhoAtomPtr) const;
        void rhoAtom_grad(const ColumnBundle& Cq, const matrix* U_rhoAtomPtr, ColumnBundle& HCq) const;
        void rhoAtom_forces(const std::vector<diagMatrix>& F, const std::vector<ColumnBundle>& C, const matrix* U_rhoAtomPtr, std::vector<vector3<> >& forces) const;
        void rhoAtom_getV(const ColumnBundle& Cq, const matrix* U_rhoAtomPtr, ColumnBundle& psi, matrix& M) const; //get DFT+U Hamiltonian in the same format as the nonlocal pseudopotential (psi = atomic orbitals, M = matrix in that order)

        //Atomic orbital related functions:
        void accumulateAtomicDensity(ScalarFieldTildeArray& nTilde) const; 
        void accumulateAtomicPotential(ScalarFieldTilde& dTilde) const; 
        void setAtomicOrbitals(ColumnBundle& Y, bool applyO, int colOffset=0) const; 
        void setAtomicOrbitals(ColumnBundle& Y, bool applyO, unsigned n, int l, int colOffset=0, int atomColStride=0) const;  
        int nAtomicOrbitals() const; 
        int lMaxAtomicOrbitals() const; 
        int nAtomicOrbitals(int l) const; 
        int atomicOrbitalOffset(unsigned iAtom, unsigned n, int l, int m, int s) const; 

        void updateLocal(ScalarFieldTilde& Vlocps, ScalarFieldTilde& rhoIon, ScalarFieldTilde& nChargeball,
                ScalarFieldTilde& nCore, ScalarFieldTilde& tauCore) const; 
        
        std::vector< vector3<> > getLocalForces(const ScalarFieldTilde& ccgrad_Vlocps, const ScalarFieldTilde& ccgrad_rhoIon,
                const ScalarFieldTilde& ccgrad_nChargeball, const ScalarFieldTilde& ccgrad_nCore, const ScalarFieldTilde& ccgrad_tauCore) const;

        void accumNonlocalForces(const ColumnBundle& Cq, const matrix& VdagC, const matrix& E_VdagC, const matrix& grad_CdagOCq, std::vector<vector3<> >& forces) const;
        
        static matrix getYlmToSpinAngleMatrix(int l, int j2); 
        static matrix getYlmOverlapMatrix(int l, int j2); 
private:
        matrix3<> Rprev; void updateLatticeDependent(); 

        RadialFunctionG VlocRadial; 
        RadialFunctionG nCoreRadial; 
        RadialFunctionG tauCoreRadial; 

        std::vector< std::vector<RadialFunctionG> > VnlRadial; 
        std::vector<matrix> Mnl; 
        matrix MnlAll; 
        
        std::vector<matrix> Qint; 
        matrix QintAll; 
        
        std::map<std::pair<vector3<>,const Basis*>, std::shared_ptr<ColumnBundle> > cachedV; //cached projectors (identified by k-point and basis pointer)
        
        struct QijIndex
        {       int l1, p1; 
                int l2, p2; 
                int l; 
                int index; 
                bool operator<(const QijIndex&) const; 
        private:
                void sortIndices(); 
        };
        std::map<QijIndex,RadialFunctionG> Qradial; 
        matrix QradialMat; 
        matrix nAug; 
        matrix E_nAug; 
        uint64_t* nagIndex; size_t* nagIndexPtr; 

        std::vector<std::vector<RadialFunctionG> > psiRadial; 
        std::vector<std::vector<RadialFunctionG> >* OpsiRadial; 
        std::vector<std::vector<double> > atomEigs; 
        
        std::vector< std::vector<int> > Vnl2j, psi2j; 
        matrix fljAll; 
        
        struct PlusU
        {       int n, l; 
                double UminusJ; 
                
                std::vector<double> Vext; 
        };
        std::vector<PlusU> plusU; 
        
        std::shared_ptr<VanDerWaals::AtomParams> vdwOverride;

        PseudopotentialFormat pspFormat;
        
        // gaussian chargeball used to prevent dielectric spikes
        // DEPRECATED: As far as possible, use partial core correction instead
        double Z_chargeball; 
        double width_chargeball; 

        const Everything* e;
        
        double tauCore_rCut; bool tauCorePlot; 
        
        void setCore(RadialFunctionR&); 
        void readFhi(istream&); // Implemented in SpeciesInfo_readFhi.cpp
        void readUspp(istream&); //Implemented in SpeciesInfo_readUspp.cpp
        void readUPF(istream&); //Implemented in SpeciesInfo_readUPF.cpp
        void setupPulay();
        
        //Following implemented in SpeciesInfo_atomFillings.cpp
        void estimateAtomEigs(); 
        void getAtom_nRadial(int spin, double magneticMoment, RadialFunctionG& nRadial, bool forceNeutral) const; 
        void getAtomPotential(RadialFunctionG& dRadial) const; 
        
        friend struct CommandIonSpecies;
        friend struct CommandSetVDW;
        friend class VanDerWaals;
        friend struct CommandAddU;
        friend struct CommandChargeball;
        friend struct CommandTauCore;
        friend struct CommandWavefunction;
        friend class WannierMinimizer;
        friend class IonicMinimizer;
        friend class Phonon;
        friend class Dump;
};


static EnumStringMap<SpeciesInfo::Constraint::ConstraintType> constraintTypeMap
(       SpeciesInfo::Constraint::None, "None",
        SpeciesInfo::Constraint::Linear, "Linear",
        SpeciesInfo::Constraint::Planar, "Planar",
        SpeciesInfo::Constraint::HyperPlane, "HyperPlane"
);

#endif // JDFTX_ELECTRONIC_SPECIESINFO_H