/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/*
 * MpqcParser.cpp
 *
 *  Created on: 12.06.2010
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <iostream>
#include <boost/tokenizer.hpp>
#include <string>

#include "CodePatterns/MemDebug.hpp"

#include "MpqcParser.hpp"

#include "atom.hpp"
#include "config.hpp"
#include "element.hpp"
#include "molecule.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/toString.hpp"
#include "CodePatterns/Verbose.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "periodentafel.hpp"
#include "World.hpp"


/** Constructor of MpqcParser.
 *
 */
MpqcParser::MpqcParser() :
  HessianPresent(false),
  theory(MBPT2)
{}

/** Destructor of MpqcParser.
 *
 */
MpqcParser::~MpqcParser()
{}

/** Load an MPQC config file into the World.
 * \param *file input stream
 */
void MpqcParser::load(istream *file)
{
  bool GeometryFollows = false;
  char line[MAXSTRINGSIZE];
  typedef boost::tokenizer<boost::char_separator<char> > tokenizer;
  boost::char_separator<char> sep("[]");
  boost::char_separator<char> whitesep(" \t");
  ConvertTo<double> toDouble;

  molecule *newmol = World::getInstance().createMolecule();
  newmol->ActiveFlag = true;
  while (file->good()) {
    file->getline(line, MAXSTRINGSIZE-1);
    std::string linestring(line);
    if ((linestring.find("atoms geometry") == string::npos) && (linestring.find("}") != string::npos)) {
      GeometryFollows = false;
    }
    if (GeometryFollows) { // we have an atom
      tokenizer tokens(linestring, sep);
//      if (tokens.size() != 2)
//        throw MpqcParseException;
      tokenizer::iterator tok_iter = tokens.begin();
      std::stringstream whitespacefilter(*tok_iter++);
      std::string element;
      whitespacefilter >> element;
      std::string vector = *tok_iter;
      tokenizer vectorcomponents(vector, whitesep);
      Vector X;
//      if (vectorcomponents.size() != NDIM)
//        throw MpqcParseException;
      tok_iter = vectorcomponents.begin();
      for (int i=0; i<NDIM; ++i) {
        X[i] = toDouble(*tok_iter++);
      }
      // create atom
      atom *newAtom = World::getInstance().createAtom();
      newAtom->setType(World::getInstance().getPeriode()->FindElement(element));
      newAtom->setPosition(X);
      newmol->AddAtom(newAtom);
      DoLog(1) && (Log() << Verbose(1) << "Adding atom " << *newAtom << std::endl);
    }
    // set some scan flags
    if (linestring.find("atoms geometry") != string::npos) {
      GeometryFollows = true;
    }
  }
}

/** Saves all atoms and data into a MPQC config file.
 * \param *file output stream
 * \param atoms atoms to store
 */
void MpqcParser::save(ostream *file, const std::vector<atom *> &atoms)
{
  Vector center;
  vector<atom *> allatoms = World::getInstance().getAllAtoms();

  // calculate center
  for (vector<atom *>::iterator runner = allatoms.begin();runner != allatoms.end(); ++runner)
    center += (*runner)->getPosition();
  center.Scale(1./allatoms.size());

  // first without hessian
  if (file->fail()) {
    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open mpqc output file." << endl);
  } else {
    *file << "% Created by MoleCuilder" << endl;
    *file << "mpqc: (" << endl;
    *file << "\tsavestate = no" << endl;
    *file << "\tdo_gradient = yes" << endl;
    if (HessianPresent) {
      *file << "\tfreq<MolecularFrequencies>: (" << endl;
      *file << "\t\tmolecule=$:molecule" << endl;
      *file << "\t)" << endl;
    }
    switch (theory) {
      case CLHF:
        *file << "\tmole<CLHF>: (" << endl;
        *file << "\t\tmolecule = $:molecule" << endl;
        *file << "\t\tbasis = $:basis" << endl;
        *file << "\t\tmemory = 16000000" << endl;
        *file << "\t)" << endl;
        break;
      case CLKS:
        *file << "\tmole<CLKS>: (" << endl;
        *file << "\t\tfunctional<StdDenFunctional>:(name=B3LYP)" << endl;
        *file << "\t\tmolecule = $:molecule" << endl;
        *file << "\t\tbasis = $:basis" << endl;
        *file << "\t\tmemory = 16000000" << endl;
        *file << "\t)" << endl;
        break;
      case MBPT2:
        *file << "\tmole<MBPT2>: (" << endl;
        *file << "\t\tbasis = $:basis" << endl;
        *file << "\t\tmolecule = $:molecule" << endl;
        *file << "\t\treference<CLHF>: (" << endl;
        *file << "\t\t\tmaxiter = 1000" << endl;
        *file << "\t\t\tbasis = $:basis" << endl;
        *file << "\t\t\tmolecule = $:molecule" << endl;
        *file << "\t\t)" << endl;
        *file << "\t)" << endl;
        break;
      case MBPT2_R12:
        *file << "\tmole<MBPT2_R12>: (" << endl;
        *file << "\t\tmolecule = $:molecule" << endl;
        *file << "\t\tbasis = $:basis" << endl;
        *file << "\t\taux_basis = $:abasis" << endl;
        *file << "\t\tstdapprox = \"A'\"" << endl;
        *file << "\t\tnfzc = 1" << endl;
        *file << "\t\tmemory = 16000000" << endl;
        *file << "\t\tintegrals<IntegralCints>:()" << endl;
        *file << "\t\treference<CLHF>: (" << endl;
        *file << "\t\t\tmolecule = $:molecule" << endl;
        *file << "\t\t\tbasis = $:basis" << endl;
        *file << "\t\t\tmaxiter = 1000" << endl;
        *file << "\t\t\tmemory = 16000000" << endl;
        *file << "\t\t\tintegrals<IntegralCints>:()" << endl;
        *file << "\t\t)" << endl;
        *file << "\t)" << endl;
        break;
      default:
        DoeLog(0) && (eLog() << Verbose(0)
            << "Unknown level of theory requested for MPQC output file." << std::endl);
        break;
    }
    *file << ")" << endl;
    *file << "molecule<Molecule>: (" << endl;
    *file << "\tunit = " << (World::getInstance().getConfig()->GetIsAngstroem() ? "angstrom" : "bohr" ) << endl;
    *file << "\t{ atoms geometry } = {" << endl;
    // output of atoms
    for (vector<atom *>::iterator AtomRunner = allatoms.begin(); AtomRunner != allatoms.end(); ++AtomRunner) {
      (*AtomRunner)->OutputMPQCLine(file, &center);
    }
    *file << "\t}" << endl;
    *file << ")" << endl;
    *file << "basis<GaussianBasisSet>: (" << endl;
    *file << "\tname = \"" << World::getInstance().getConfig()->basis << "\"" << endl;
    *file << "\tmolecule = $:molecule" << endl;
    *file << ")" << endl;
    if (theory == MBPT2_R12) {
      *file << "% auxiliary basis set specification" << endl;
      *file << "\tabasis<GaussianBasisSet>: (" << endl;
      *file << "\tname = \"aug-cc-pVDZ\"" << endl;
      *file << "\tmolecule = $:molecule" << endl;
      *file << ")" << endl;
    }
  }
}

/** Sets whether hessian is desired or not
 * \param hessian statement
 */
void MpqcParser::setHessian(bool hessian)
{
  HessianPresent = hessian;
}

/** Sets the desired level of solving theory to use
 * \param _theory shorthand of the theory
 */
void MpqcParser::setTheory(enum Theory _theory)
{
  theory = _theory;
}