RadialSchrodinger.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_ELECTRONIC_RADIALSCHRODINGER_H
#define JDFTX_ELECTRONIC_RADIALSCHRODINGER_H

#include <electronic/RadialFunction.h>
#include <electronic/matrix.h>
#include <core/Minimize.h>
#include <core/scalar.h>
#include <cstring>
#include <vector>
#include <map>

class RadialSchrodinger
{
public:
        RadialSchrodinger(const std::vector<double>& rArr, const std::vector<double>& drArr,
                const std::vector<double>& V, double Z, size_t iMatch=0.);
        
        std::vector< std::vector<double> > getFillings(double nElectrons,
                std::vector< std::vector<double> >* Eptr=0);
        
        struct Outputs
        {
                std::vector<double>* n; 
                std::vector<double>* Dn; 
                std::vector<double>* tau; 
                std::vector< std::vector< std::vector<complex> > >* z; 
                std::vector< std::vector<double> >* E; 
                
                Outputs(std::vector<double>* n=0, std::vector<double>* Dn=0, std::vector<double>* tau=0,
                        std::vector< std::vector< std::vector<complex> > >* z=0, std::vector< std::vector<double> >* E=0)
                : n(n), Dn(Dn), tau(tau), z(z), E(E) {}
        };
        
        double compute(const std::vector< std::vector<double> >& F, Outputs outputs);
        
private:
        const std::vector<double>& rArr; 
        const std::vector<double>& drArr; 
        std::vector<double> Vsub; 
        double Z; 
        size_t iMatch; 
        
        std::vector< std::map<int,double> > nodesEmin, nodesEmax; 
        std::vector< std::map<double,double> > cachedEerr; 
        std::vector< std::map<int,double> > cachedE; 

        inline complex schEqn(double r, int iSub, complex z, int l, double E) const
        {       double rInv = r ? 1./r : 0.;
                double V = Vsub[iSub] - Z*rInv;
                return complex(z.imag(), 2*(z.real()*(V-E) - z.imag()*(l+1)*rInv));
        }
        
        int solveSchEqn(int l, double E, complex* zSave=0);
        
        void locateNodeCount(int l, int nNodes);
        
        double getEerr(int l, double E);
        
        double findE(int l, double Elo, double Ehi, double tol);

        double getEig(int l, int nNodes, complex* zEvec=0);
};


class InvertKS : public Minimizable<diagMatrix>
{
public:
        InvertKS(const RadialFunctionR& nTarget);
        RadialFunctionR getTau(); 

        //Functions for Minimizable interface:
        void step(const diagMatrix& dV, double alpha);
        double compute(diagMatrix* E_V);
        diagMatrix precondition(const diagMatrix& grad);
        
private:
        diagMatrix V; 
        double nElectrons; 
        const std::vector<double> &r, &dr, &n0; 
        std::vector<double> n; 
        std::vector< std::vector<double> > F; 
};

#endif // JDFTX_ELECTRONIC_RADIALSCHRODINGER_H