Coulomb.h

/*-------------------------------------------------------------------
Copyright 2012 Ravishankar Sundararaman

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_CORE_COULOMB_H
#define JDFTX_CORE_COULOMB_H

#include <core/GridInfo.h>
#include <core/string.h>
#include <memory>
#include <set>


struct CoulombParams
{       
        enum Geometry
        {       Periodic, 
                Slab, 
                Wire, 
                Cylindrical, 
                Isolated, 
                Spherical 
        };
        Geometry geometry; 
        int iDir; 
        double Rc; 
        
        double ionMargin; 
        
        bool embed; 
        vector3<> embedCenter; 
        bool embedFluidMode; 
        
        vector3<> Efield; 
        
        //Parameters for computing exchange integrals:
        enum ExchangeRegularization
        {       None, 
                AuxiliaryFunction, 
                ProbeChargeEwald, 
                SphericalTruncated, 
                WignerSeitzTruncated 
        };
        ExchangeRegularization exchangeRegularization; 
        std::set<double> omegaSet; 
        std::shared_ptr<struct Supercell> supercell; 
        
        CoulombParams();
        
        std::shared_ptr<class Coulomb> createCoulomb(const GridInfo& gInfo) const;
        
        vector3<bool> isTruncated() const;
};


struct Atom
{       double Z; 
        vector3<> pos; 
        vector3<> force; 
        int atomicNumber; 
        int sp; 
        
        Atom(double Z, vector3<> pos, vector3<> force=vector3<>(0,0,0), int atomicNumber=0, int sp=0)
        : Z(Z), pos(pos), force(force), atomicNumber(atomicNumber), sp(sp)
        {
        }
};


class Ewald
{
public:
        virtual double energyAndGrad(std::vector<Atom>& atoms) const=0;
};


class Coulomb
{
public:
        enum PointChargeMode
        {       PointChargeNone, 
                PointChargeLeft, 
                PointChargeRight 
        };
        
        ScalarFieldTilde operator()(ScalarFieldTilde&&, PointChargeMode pointChargeMode=PointChargeNone) const;
        
        ScalarFieldTilde operator()(const ScalarFieldTilde&, PointChargeMode pointChargeMode=PointChargeNone) const;
        
        double energyAndGrad(std::vector<Atom>& atoms) const; 

        ScalarField getEfieldPotential() const;
        
        complexScalarFieldTilde operator()(complexScalarFieldTilde&&, vector3<> kDiff, double omega) const;

        complexScalarFieldTilde operator()(const complexScalarFieldTilde&, vector3<> kDiff, double omega) const;

private:
        const GridInfo& gInfoOrig; 
protected:
        const CoulombParams& params;
        const GridInfo& gInfo; 
        std::shared_ptr<Ewald> ewald;
        std::map<double, std::shared_ptr<struct ExchangeEval>> exchangeEval;
        friend struct ExchangeEval;
        
        Coulomb(const GridInfo& gInfoOrig, const CoulombParams& params);
        virtual ~Coulomb();
        
        void initExchangeEval();
        
        virtual ScalarFieldTilde apply(ScalarFieldTilde&&) const=0;
        
        virtual std::shared_ptr<Ewald> createEwald(matrix3<> R, size_t nAtoms) const=0;

private:
        //Data for mesh-embedded truncation:
        GridInfo gInfoEmbed; 
        vector3<int> ivCenter; 
        vector3<> xCenter; 
        vector3<> embedScale; 
        int* embedIndex; 
        std::vector< std::pair<int,int*> > symmIndex; 
        class WignerSeitz* wsOrig; 
        double ionWidth; 
        RealKernel* ionKernel;
        ScalarFieldTilde embedExpand(const ScalarFieldTilde& in) const; 
        complexScalarFieldTilde embedExpand(complexScalarFieldTilde&& in) const; 
        ScalarFieldTilde embedShrink(const ScalarFieldTilde& in) const; 
        complexScalarFieldTilde embedShrink(complexScalarFieldTilde&& in) const; 
        friend struct FluidSolver;
        friend struct SlabEpsilon;
        friend struct ChargedDefect;
};

#endif // JDFTX_CORE_COULOMB_H