/*
 *            Copyright 2009-2020 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.
 *
 */

// Third party includes
#include <boost/algorithm/string.hpp>

// Local VOTCA includes
#include "votca/tools/getline.h"
#include "votca/tools/globals.h"
#include "votca/xtp/ecpaobasis.h"
#include "votca/xtp/orbitals.h"
#include "votca/xtp/qmpackage.h"
#include "votca/xtp/qmpackagefactory.h"

namespace votca {
namespace xtp {
using std::flush;

tools::Property QMPackage::ParseCommonOptions(const tools::Property& options) {

  std::string key = "package";

  _settings.read_property(options, key);

  if (tools::VotcaShareSet()) {
    Settings qmpackage_defaults{key};
    qmpackage_defaults.load_from_xml(this->FindDefaultsFile());
    _settings.amend(qmpackage_defaults);
  } else {
    std::cout << "Warning: VOTCASHARE environment variable not defined\n";
  }
  _settings.validate();

  _charge = _settings.get<Index>("charge");
  _spin = _settings.get<Index>("spin");
  _basisset_name = _settings.get("basisset");

  if (_settings.has_key("cleanup")) {
    _cleanup = _settings.get("cleanup");
  }

  if (getPackageName() != "xtp") {
    _scratch_dir = _settings.get("scratch");
  }
  return _settings.to_property("package");
}

void QMPackage::ReorderOutput(Orbitals& orbitals) const {
  if (!orbitals.hasQMAtoms()) {
    throw std::runtime_error("Orbitals object has no QMAtoms");
  }

  AOBasis dftbasis = orbitals.SetupDftBasis();
  // necessary to update nuclear charges on qmatoms
  if (orbitals.hasECPName()) {
    ECPBasisSet ecps;
    ecps.Load(orbitals.getECPName());
    ECPAOBasis ecpbasis;
    ecpbasis.Fill(ecps, orbitals.QMAtoms());
  }

  if (orbitals.hasMOs()) {
    OrbReorder reorder(ShellReorder(), ShellMulitplier());
    reorder.reorderOrbitals(orbitals.MOs().eigenvectors(), dftbasis);
    XTP_LOG(Log::info, *_pLog) << "Reordered MOs" << flush;
  }

  return;
}

Eigen::MatrixXd QMPackage::ReorderMOsBack(const Orbitals& orbitals) const {
  if (!orbitals.hasQMAtoms()) {
    throw std::runtime_error("Orbitals object has no QMAtoms");
  }
  AOBasis dftbasis = orbitals.SetupDftBasis();
  Eigen::MatrixXd result = orbitals.MOs().eigenvectors();
  bool reverseOrder = true;
  OrbReorder reorder(ShellReorder(), ShellMulitplier(), reverseOrder);
  reorder.reorderOrbitals(result, dftbasis);
  return result;
}

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

  std::vector<QMPackage::MinimalMMCharge> multipoles_split;
  // Calculate virtual charge positions
  double a = _settings.get<double>("dipole_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;
  if (std::abs(qA) > 1e-12) {
    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 a2 = 2 * a;
    for (Index i = 0; i < 3; i++) {
      double q = es.eigenvalues()[i] / (a2 * a2);
      const Eigen::Vector3d vec1 =
          aps.getPos() + 0.5 * a2 * es.eigenvectors().col(i);
      const Eigen::Vector3d vec2 =
          aps.getPos() - 0.5 * a2 * es.eigenvectors().col(i);
      multipoles_split.push_back(MinimalMMCharge(vec1, q));
      multipoles_split.push_back(MinimalMMCharge(vec2, q));
    }
  }
  return multipoles_split;
}

std::vector<std::string> QMPackage::GetLineAndSplit(
    std::ifstream& input_file, const std::string separators) const {
  std::string line;
  tools::getline(input_file, line);
  boost::trim(line);
  tools::Tokenizer tok(line, separators);
  return tok.ToVector();
}

std::string QMPackage::FindDefaultsFile() const {
  return tools::GetVotcaShare() + "/xtp/data/qmpackage_defaults.xml";
}

}  // namespace xtp
}  // namespace votca