#include "qc/basis.hpp"
#include "qc/hf.hpp"
#include "qc/integrals.hpp"

#include <cmath>
#include <cstdio>

int main() {
  const double bond_bohr = 1.401;
  qc::Basis basis = qc::sto3g_hydrogen_dimer(bond_bohr);
  std::vector<qc::Atom> atoms;
  atoms.push_back({1.0, {-0.5 * bond_bohr, 0.0, 0.0}});
  atoms.push_back({1.0, {0.5 * bond_bohr, 0.0, 0.0}});

  const qc::MolecularIntegrals ints = qc::compute_integrals(basis, atoms);
  const qc::RhfResult rhf = qc::rhf_closed_shell(ints, 2, 1e-10, 100, 6, &atoms);

  std::printf("H2 STO-3G RHF (R = %.3f bohr)\n", bond_bohr);
  std::printf("iterations: %d  converged: %s\n", rhf.iterations,
              rhf.converged ? "yes" : "no");
  std::printf("electronic energy: %.10f Hartree\n", rhf.electronic_energy);
  std::printf("nuclear repulsion: %.10f Hartree\n", rhf.nuclear_repulsion);
  std::printf("total energy:      %.10f Hartree\n", rhf.total_energy);
  std::printf("orbital energies (Ha):");
  for (std::size_t i = 0; i < rhf.orbital_energies.size(); ++i) {
    std::printf(" %.8f", rhf.orbital_energies[i]);
  }
  std::printf("\n");
  return 0;
}