/*
 *            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.
 *
 */

#include "votca/xtp/qmpackage.h"

#include <boost/algorithm/string.hpp>

namespace votca {
    namespace xtp {
      using std::flush;
        void QMPackage::ReorderOutput(Orbitals& orbitals) {
            BasisSet dftbasisset;
            dftbasisset.LoadBasisSet(_basisset_name);
            if (!orbitals.hasQMAtoms()) {
                throw std::runtime_error("Orbitals object has no QMAtoms");
            }

            AOBasis dftbasis;
            dftbasis.AOBasisFill(dftbasisset, orbitals.QMAtoms());
            //necessary to update nuclear charges on qmatoms
            if (_write_pseudopotentials) {
                BasisSet ecps;
                ecps.LoadPseudopotentialSet(_ecp_name);
                AOBasis ecpbasis;
                ecpbasis.ECPFill(ecps, orbitals.QMAtoms());
            }

            if (orbitals.hasAOOverlap()) {
                dftbasis.ReorderMatrix(orbitals.AOOverlap(), getPackageName(), "xtp");
                XTP_LOG(logDEBUG, *_pLog) << "Reordered Overlap matrix" << flush;
            }
            if (orbitals.hasAOVxc()) {
                dftbasis.ReorderMatrix(orbitals.AOVxc(), getPackageName(), "xtp");
                XTP_LOG(logDEBUG, *_pLog) << "Reordered VXC matrix" << flush;
            }
            if (orbitals.hasMOCoefficients()) {
                dftbasis.ReorderMOs(orbitals.MOCoefficients(), getPackageName(), "xtp");
                XTP_LOG(logDEBUG, *_pLog) << "Reordered MOs" << flush;
            }

            return;
        }

        Eigen::MatrixXd QMPackage::ReorderMOsBack(const Orbitals& orbitals) const{
            BasisSet dftbasisset;
            dftbasisset.LoadBasisSet(_basisset_name);
            if (!orbitals.hasQMAtoms()) {
                throw std::runtime_error("Orbitals object has no QMAtoms");
            }
            AOBasis dftbasis;
            dftbasis.AOBasisFill(dftbasisset, orbitals.QMAtoms());
            Eigen::MatrixXd result=orbitals.MOCoefficients();
            dftbasis.ReorderMOs(result, "xtp", getPackageName());
            return result;
        }

        std::vector<QMPackage::MinimalMMCharge > QMPackage::SplitMultipoles(const PolarSite& aps) {

            std::vector< QMPackage::MinimalMMCharge > multipoles_split;
            // Calculate virtual charge positions
            double a = _dpl_spacing; // this is in a0
            double mag_d = aps.getDipole().norm();// this is in e * a0
            const Eigen::Vector3d dir_d = aps.getDipole().normalized();
            const Eigen::Vector3d A = aps.getPos() + 0.5 * a * dir_d;
            const Eigen::Vector3d B = aps.getPos() - 0.5 * a * dir_d;
            double qA = mag_d / a;
            double qB = -qA;
            multipoles_split.push_back(MinimalMMCharge(A, qA));
            multipoles_split.push_back(MinimalMMCharge(B, qB));


            if (aps.getRank() > 1) {
                const Eigen::Matrix3d components = aps.CalculateCartesianMultipole();
                Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> es;
                es.computeDirect(components);
                double a = 2 * _dpl_spacing;
                for (int i = 0; i < 3; i++) {
                    double q = es.eigenvalues()[i] / (a * a);
                    const Eigen::Vector3d vec1 = aps.getPos() + 0.5 * a * es.eigenvectors().col(i);
                    const Eigen::Vector3d vec2 = aps.getPos() - 0.5 * a * es.eigenvectors().col(i);
                    multipoles_split.push_back(MinimalMMCharge(vec1, q));
                    multipoles_split.push_back(MinimalMMCharge(vec2, q));
                }
            }
            return multipoles_split;
        }
        
      void QMPackage::setMultipoleBackground(const std::shared_ptr<MMRegion>& PolarSegments ) {
        if(PolarSegments->size()==0){
          std::cout<<"WARNING::The Multipole Background has no entries!"<<std::endl;
          return;
        }
      _PolarSegments = PolarSegments;
      _write_charges = true;
      
      WriteChargeOption();
    }
      
      std::vector<std::string> QMPackage::GetLineAndSplit(std::ifstream& input_file,const std::string separators ){
          std::string line;
          getline(input_file, line);
          boost::trim(line);
          std::vector<std::string> row;
          boost::algorithm::split(row, line, boost::is_any_of(separators), boost::algorithm::token_compress_on);
          return row;
        }


    }
}