/* 
 *            Copyright 2009-2018 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

#ifndef VOTCA_XTP_AOSHELL_H
#define	VOTCA_XTP_AOSHELL_H


#include <boost/math/constants/constants.hpp>
#include <votca/xtp/eigen.h>
#include <votca/tools/constants.h>

#include <votca/xtp/basisset.h>
#include "qmatom.h"

namespace votca { namespace xtp {

class AOBasis;
class AOShell;  

class AOGaussianPrimitive 
{
    friend class AOShell;
public:

    double getPowfactor()const {return _powfactor;}
    int    getPower()const{return _power;}
    double getDecay()const {return _decay;}
    const std::vector<double>& getContraction()const {return _contraction;}
    const AOShell& getShell() const{return _aoshell;}
private:
     
    int _power; // used in pseudopotenials only
    double _decay;
    std::vector<double> _contraction;
    const AOShell& _aoshell;
    double _powfactor;//used in evalspace to speed up DFT
    // private constructor, only a shell can create a primitive
    AOGaussianPrimitive( const GaussianPrimitive& gaussian, const AOShell& aoshell )
    : _power(gaussian._power),
    _decay(gaussian._decay),
    _contraction(gaussian._contraction),
    _aoshell(aoshell) {_powfactor=std::pow(2.0 * _decay / boost::math::constants::pi<double>(), 0.75) ; }

    AOGaussianPrimitive( const AOGaussianPrimitive& gaussian, const AOShell &aoshell )
    : _power(gaussian._power),
    _decay(gaussian._decay),
    _contraction(gaussian._contraction),
    _aoshell(aoshell),
    _powfactor(gaussian._powfactor){;}

};      
    
/*
 * shells in a Gaussian-basis expansion
 */
class AOShell 
{
    friend class AOBasis;
public:

    AOShell(const AOShell& shell){
          
   _type= shell._type;
    _Lmax= shell._Lmax;
    _scale= shell._scale;
    _numFunc= shell._numFunc;
    _mindecay= shell._mindecay;
    _startIndex= shell._startIndex;
    _offset= shell._offset;
    _pos= shell._pos;
    _atomindex=shell._atomindex;
    _nonlocal=shell._nonlocal;
    _gaussians.reserve(shell._gaussians.size());
    for(const auto& gaus:shell._gaussians){
        _gaussians.push_back(AOGaussianPrimitive(gaus,*this));
    }
        
    }

    const std::string& getType() const{ return _type; }
    int    getNumFunc() const{ return _numFunc ;}
    int    getStartIndex() const{ return _startIndex ;}
    int    getOffset() const{ return _offset ;}
    int    getAtomIndex() const{ return _atomindex;}
    
    int getLmax() const{ return _Lmax;}
    
    bool isCombined()const{
        return _type.length()>1;
    }
    
    bool isNonLocal(  ) const{ return _nonlocal;}
    
    const Eigen::Vector3d& getPos() const{ return _pos; }
    double getScale() const{ return _scale; }
    
    int getSize() const{ return _gaussians.size(); }
    
    void CalcMinDecay(){
     _mindecay=std::numeric_limits<double>::max();
     for(auto& gaussian:_gaussians){
         if(gaussian.getDecay()<_mindecay){
             _mindecay=gaussian.getDecay();
         }
     }
     return;
    }
    
    double getMinDecay() const{return _mindecay;}
    
    
  void EvalAOspace(Eigen::VectorBlock<Eigen::VectorXd>&  AOvalues, const Eigen::Vector3d& grid_pos ) const;
  void EvalAOspace(Eigen::VectorBlock<Eigen::VectorXd>&  AOvalues,Eigen::Block< Eigen::MatrixX3d >& AODervalues, const Eigen::Vector3d& grid_pos ) const;

    // iterator over pairs (decay constant; contraction coefficient)
    typedef std::vector< AOGaussianPrimitive >::const_iterator GaussianIterator;
    GaussianIterator begin() const{ return _gaussians.begin(); }
    GaussianIterator end()const{ return _gaussians.end(); }
   
    // adds a Gaussian 
    void  addGaussian( const GaussianPrimitive& gaussian ){
        AOGaussianPrimitive aogaussian = AOGaussianPrimitive(gaussian, *this);
        _gaussians.emplace_back( aogaussian );
        return;
    }                                                                  

    void normalizeContraction();
    
    friend std::ostream &operator<<(std::ostream &out, const AOShell& shell);
    
private:   

    // only class aobasis can construct shells    
    AOShell( const Shell& shell, const QMAtom & atom, int startIndex)
            : _type(shell.getType()),_Lmax(shell.getLmax()),
                    _scale(shell.getScale()), _numFunc(shell.getnumofFunc()),
                    _startIndex(startIndex), _offset(shell.getOffset()), _pos(atom.getPos()) , 
                    _atomindex(atom.getAtomID()) { ; }
    // for ECPs
    AOShell( const Shell& shell, const QMAtom & atom, int startIndex, bool nonlocal)
            : _type(shell.getType()),_Lmax(shell.getLmax()),
                    _scale(shell.getScale()), _numFunc(shell.getnumofFunc()),
                    _startIndex(startIndex), _offset(shell.getOffset()), _pos(atom.getPos()) , 
                    _atomindex(atom.getAtomID()),_nonlocal(nonlocal) { ; }
            

    
    // shell type (S, P, D))
    std::string _type;
    int _Lmax;
    // scaling factor
    double _scale;
    // number of functions in shell
    int _numFunc;
    double _mindecay;
    int _startIndex;
    int _offset;
    Eigen::Vector3d _pos;
    int _atomindex;
    //used for ecp calculations
    bool _nonlocal;
     
    // vector of pairs of decay constants and contraction coefficients
    std::vector< AOGaussianPrimitive > _gaussians;
    
};

    
}}

#endif	// VOTCA_XTP_AOSHELL_H