IonInfo.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_IONINFO_H
#define JDFTX_ELECTRONIC_IONINFO_H

#include <electronic/common.h>
#include <electronic/SpeciesInfo.h>
#include <electronic/IonicMinimizer.h>
#include <electronic/matrix.h>
#include <core/ScalarField.h>
#include <vector>
#include <core/Thread.h>

enum CoordsType {CoordsLattice, CoordsCartesian}; 

enum ForcesOutputCoords { ForcesCoordsPositions, ForcesCoordsLattice, ForcesCoordsCartesian, ForcesCoordsContravariant };
static EnumStringMap<ForcesOutputCoords> forcesOutputCoordsMap(
        ForcesCoordsPositions, "Positions",
        ForcesCoordsLattice, "Lattice",
        ForcesCoordsCartesian, "Cartesian",
        ForcesCoordsContravariant, "Contravariant");

enum coreOverlapCheck { additive, vector, none };

class IonInfo
{
public:
        std::vector< std::shared_ptr<SpeciesInfo> > species; 
        std::vector<string> pspFilenamePatterns; 
        
        CoordsType coordsType; 
        ForcesOutputCoords forcesOutputCoords; 
        coreOverlapCheck coreOverlapCondition; 
        bool vdWenable; 
        double vdWscale; 
        
        IonicGradient forces; 
        
        ScalarFieldTilde Vlocps; 
        ScalarFieldTilde rhoIon; 
        ScalarFieldTilde nChargeball; 
        ScalarField nCore; 
        ScalarField tauCore; 
        
        IonInfo();
        
        void setup(const Everything&);
        void printPositions(FILE*) const; 
        bool checkPositions() const; 
        double getZtot() const; 
        
        void update(Energies&); 

        double ionicEnergyAndGrad(IonicGradient& forces) const;

        double EnlAndGrad(const QuantumNumber& qnum, const diagMatrix& Fq, const std::vector<matrix>& VdagCq, std::vector<matrix>& HVdagCq) const;
        
        void augmentOverlap(const ColumnBundle& Cq, ColumnBundle& OCq, std::vector<matrix>* VdagCq=0) const;
        
        //Multi-stage density augmentation and gradient propagation (see corresponding functions in SpeciesInfo)
        void augmentDensityInit() const;
        void augmentDensitySpherical(const QuantumNumber& qnum, const diagMatrix& Fq, const std::vector<matrix>& VdagCq) const;
        void augmentDensityGrid(ScalarFieldArray& n) const;
        void augmentDensityGridGrad(const ScalarFieldArray& E_n, IonicGradient* forces=0) const;
        void augmentDensitySphericalGrad(const QuantumNumber& qnum, const diagMatrix& Fq, const std::vector<matrix>& VdagCq, std::vector<matrix>& HVdagCq) const;
        
        void project(const ColumnBundle& Cq, std::vector<matrix>& VdagCq, matrix* rotExisting=0) const; //Update pseudopotential projections (optionally retain non-zero ones with specified rotation)
        void projectGrad(const std::vector<matrix>& HVdagCq, const ColumnBundle& Cq, ColumnBundle& HCq) const; //Propagate projected gradient (HVdagCq) to full gradient (HCq)
        
        //rhoAtom is a flat array of atomic density matrices per U type, with index order (outer to inner): species, Uparam(n,l), spin, atom
        size_t rhoAtom_nMatrices() const; 
        void rhoAtom_initZero(std::vector<matrix>& rhoAtom) const; 
        void rhoAtom_calc(const std::vector<diagMatrix>& F, const std::vector<ColumnBundle>& C, std::vector<matrix>& rhoAtom) const; 
        double rhoAtom_computeU(const std::vector<matrix>& rhoAtom, std::vector<matrix>& U_rhoAtom) const; 
        void rhoAtom_grad(const ColumnBundle& Cq, const std::vector<matrix>& U_rhoAtom, ColumnBundle& HCq) const; 
        void rhoAtom_forces(const std::vector<diagMatrix>& F, const std::vector<ColumnBundle>& C, const std::vector<matrix>& U_rhoAtom, IonicGradient& forces) const; 
        void rhoAtom_getV(const ColumnBundle& Cq, const std::vector<matrix>& U_rhoAtom, std::vector<ColumnBundle>& psi, std::vector<matrix>& M) const; 
        
        ColumnBundle rHcommutator(const ColumnBundle &Y, int iDir) const; 
        
        int nAtomicOrbitals() const; 
        ColumnBundle getAtomicOrbitals(int q, bool applyO, int extraCols=0) const; 
        
        enum IonWidthMethod
        {       IonWidthEcut, 
                IonWidthFFTbox, 
                IonWidthManual 
        }
        ionWidthMethod; 
        double ionWidth; 
        bool shouldPrintForceComponents;

private:
        const Everything* e;
        
        void pairPotentialsAndGrad(Energies* ener=0, IonicGradient* forces=0) const;
};

#endif // JDFTX_ELECTRONIC_IONINFO_H