/*
 * 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.
 *
 *     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 <libint2/initialize.h>
#define BOOST_TEST_MAIN

#define BOOST_TEST_MODULE aomatrix_test

// Third party includes
#include <boost/test/unit_test.hpp>

// Local VOTCA includes
#include "votca/xtp/aomatrix.h"
#include "votca/xtp/orbitals.h"
#include <votca/tools/eigenio_matrixmarket.h>

using namespace votca::xtp;
using namespace votca;
using namespace std;

BOOST_AUTO_TEST_SUITE(aomatrix_test)

QMMolecule Methane() {

  QMMolecule mol(" ", 0);
  mol.LoadFromFile(std::string(XTP_TEST_DATA_FOLDER) +
                   "/aomatrix/molecule.xyz");
  return mol;
}

BOOST_AUTO_TEST_CASE(aomatrices_test) {
  libint2::initialize();
  QMMolecule mol = Methane();
  BasisSet basis;
  basis.Load(std::string(XTP_TEST_DATA_FOLDER) + "/aomatrix/3-21G.xml");
  AOBasis aobasis;
  aobasis.Fill(basis, mol);
  AOOverlap overlap;
  overlap.Fill(aobasis);
  Eigen::MatrixXd overlap_ref = votca::tools::EigenIO_MatrixMarket::ReadMatrix(
      std::string(XTP_TEST_DATA_FOLDER) + "/aomatrix/overlap_ref.mm");

  bool check_overlap = overlap.Matrix().isApprox(overlap_ref, 0.0001);
  BOOST_CHECK_EQUAL(check_overlap, 1);
  if (!check_overlap) {
    cout << "ref" << endl;
    cout << overlap_ref << endl;
    cout << "result" << endl;
    cout << overlap.Matrix() << endl;
  }

  AOKinetic kinetic;
  kinetic.Fill(aobasis);
  Eigen::MatrixXd kinetic_ref = votca::tools::EigenIO_MatrixMarket::ReadMatrix(
      std::string(XTP_TEST_DATA_FOLDER) + "/aomatrix/kinetic_ref.mm");

  bool check_kinetic = kinetic.Matrix().isApprox(kinetic_ref, 0.00001);
  BOOST_CHECK_EQUAL(check_kinetic, 1);
  if (!check_kinetic) {
    cout << "ref" << endl;
    cout << kinetic_ref << endl;
    cout << "result" << endl;
    cout << kinetic.Matrix() << endl;
  }

  AOCoulomb coulomb;
  coulomb.Fill(aobasis);
  Eigen::MatrixXd coulomb_ref = votca::tools::EigenIO_MatrixMarket::ReadMatrix(
      std::string(XTP_TEST_DATA_FOLDER) + "/aomatrix/coulomb_ref.mm");
  bool check_coulomb = coulomb.Matrix().isApprox(coulomb_ref, 0.00001);
  BOOST_CHECK_EQUAL(check_coulomb, 1);
  if (!check_coulomb) {
    cout << "ref" << endl;
    cout << coulomb_ref << endl;
    cout << "result" << endl;
    cout << coulomb.Matrix() << endl;
  }

  Eigen::MatrixXd ps_invSqrt = coulomb.Pseudo_InvSqrt(1e-7);
  Eigen::MatrixXd coulombinvsqrt_ref =
      votca::tools::EigenIO_MatrixMarket::ReadMatrix(
          std::string(XTP_TEST_DATA_FOLDER) +
          "/aomatrix/coulombinvsqrt_ref.mm");

  bool check_coulombinvsqrt = ps_invSqrt.isApprox(coulombinvsqrt_ref, 0.00001);
  BOOST_CHECK_EQUAL(check_coulombinvsqrt, 1);
  if (!check_coulombinvsqrt) {
    cout << "ref" << endl;
    cout << coulombinvsqrt_ref << endl;
    cout << "result" << endl;
    cout << ps_invSqrt << endl;
  }

  Eigen::MatrixXd ps_invSqrtgw = coulomb.Pseudo_InvSqrt_GWBSE(overlap, 1e-7);
  Eigen::MatrixXd coulombinvsqrtgw_ref =
      votca::tools::EigenIO_MatrixMarket::ReadMatrix(
          std::string(XTP_TEST_DATA_FOLDER) +
          "/aomatrix/coulombinvsqrtgw_ref.mm");

  bool check_coulombinvsqrtgw =
      ps_invSqrtgw.isApprox(coulombinvsqrtgw_ref, 0.00001);

  BOOST_CHECK_EQUAL(check_coulombinvsqrtgw, 1);
  if (!check_coulombinvsqrtgw) {
    cout << "ref" << endl;
    cout << coulombinvsqrtgw_ref << endl;
    cout << "result" << endl;
    cout << ps_invSqrtgw << endl;
  }
  libint2::finalize();
}

BOOST_AUTO_TEST_CASE(aomatrices_contracted_test) {
  libint2::initialize();
  QMMolecule mol("C", 0);
  mol.LoadFromFile(std::string(XTP_TEST_DATA_FOLDER) + "/aomatrix/C.xyz");
  BasisSet basis;
  basis.Load(std::string(XTP_TEST_DATA_FOLDER) + "/aomatrix/contracted.xml");
  AOBasis aobasis;
  aobasis.Fill(basis, mol);
  AOOverlap overlap;
  overlap.Fill(aobasis);
  Eigen::MatrixXd overlap_ref = votca::tools::EigenIO_MatrixMarket::ReadMatrix(
      std::string(XTP_TEST_DATA_FOLDER) +
      "/aomatrix/overlap_ref_contracted.mm");

  bool check_overlap = overlap.Matrix().isApprox(overlap_ref, 0.0001);
  if (!check_overlap) {
    std::cout << std::endl;
    std::cout << "Ref" << std::endl;
    std::cout << overlap_ref << std::endl;
    std::cout << "Result" << std::endl;
    std::cout << overlap.Matrix() << std::endl;
  }
  BOOST_CHECK_EQUAL(check_overlap, 1);
  libint2::finalize();
}

BOOST_AUTO_TEST_CASE(aocoulomb_inv_test) {
  libint2::initialize();
  QMMolecule mol = Methane();
  BasisSet basis;
  basis.Load(std::string(XTP_TEST_DATA_FOLDER) + "/aomatrix/3-21G.xml");
  AOBasis aobasis;
  aobasis.Fill(basis, mol);

  AOCoulomb cou;
  cou.Fill(aobasis);

  Eigen::MatrixXd PseudoInvSqrt = cou.Pseudo_InvSqrt(1e-7);

  Eigen::MatrixXd Reformed = PseudoInvSqrt * PseudoInvSqrt * cou.Matrix();

  bool check_inv = Reformed.isApprox(Eigen::MatrixXd::Identity(17, 17), 0.0001);
  if (!check_inv) {
    std::cout << "reformed" << endl;
    std::cout << Reformed << endl;
  }
  BOOST_CHECK_EQUAL(check_inv, 1);
  libint2::finalize();
}

/*BOOST_AUTO_TEST_CASE(large_l_test) {

  QMMolecule mol("C", 0);
  mol.LoadFromFile(std::string(XTP_TEST_DATA_FOLDER) + "/aomatrix/C2.xyz");

  BasisSet basisset;
  basisset.Load(std::string(XTP_TEST_DATA_FOLDER) + "/aomatrix/G.xml");

  BasisSet auxbasisset;
  auxbasisset.Load(std::string(XTP_TEST_DATA_FOLDER) + "/aomatrix/I.xml");
  AOBasis dftbasis;
  dftbasis.Fill(basisset, mol);

  AOBasis auxbasis;
  auxbasis.Fill(auxbasisset, mol);

  AOOverlap overlap;
  overlap.Fill(auxbasis);

  Eigen::MatrixXd overlap_ref = votca::tools::EigenIO_MatrixMarket::ReadMatrix(
      std::string(XTP_TEST_DATA_FOLDER) + "/aomatrix/overlap_ref_gi.mm");

  bool check_overlap = overlap.Matrix().isApprox(overlap_ref, 0.00001);

  BOOST_CHECK_EQUAL(check_overlap, 1);
  if (!check_overlap) {
    cout << "ref" << endl;
    cout << overlap_ref << endl;
    cout << "result" << endl;
    cout << overlap.Matrix() << endl;
  }

  AOCoulomb coulomb;
  coulomb.Fill(auxbasis);
  Eigen::MatrixXd coulomb_ref = votca::tools::EigenIO_MatrixMarket::ReadMatrix(
      std::string(XTP_TEST_DATA_FOLDER) + "/aomatrix/coulomb_ref_gi.mm");

  bool check_coulomb = coulomb.Matrix().isApprox(coulomb_ref, 0.00001);
  BOOST_CHECK_EQUAL(check_coulomb, 1);
  if (!check_coulomb) {
    cout << "ref" << endl;
    cout << coulomb_ref << endl;
    cout << "result" << endl;
    cout << coulomb.Matrix() << endl;
  }

  AOKinetic kinetic;
  kinetic.Fill(dftbasis);

  Eigen::MatrixXd kinetic_ref = votca::tools::EigenIO_MatrixMarket::ReadMatrix(
      std::string(XTP_TEST_DATA_FOLDER) + "/aomatrix/kinetic_ref_gi.mm");

  bool check_kinetic = kinetic.Matrix().isApprox(kinetic_ref, 0.00001);

  BOOST_CHECK_EQUAL(check_kinetic, 1);
  if (!check_kinetic) {
    cout << "ref" << endl;
    cout << kinetic_ref << endl;
    cout << "result" << endl;
    cout << kinetic.Matrix() << endl;
  }
} */

BOOST_AUTO_TEST_SUITE_END()