/*
 * 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 "Helpers/MemDebug.hpp"

#include "MpqcParser.hpp"

#include "atom.hpp"
#include "config.hpp"
#include "element.hpp"
#include "Helpers/Log.hpp"
#include "Helpers/Verbose.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "periodentafel.hpp"
#include "World.hpp"


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

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

/** Load an MPQC config file into the World.
 * \param *file input stream
 */
void MpqcParser::load(istream *file)
{
  // TODO: MpqcParser::load implementation
  ASSERT(false, "Not implemented yet");
}

void MpqcParser::save(ostream *file)
{
  DoLog(0) && (Log() << Verbose(0) << "Saving changes to MPQC ." << std::endl);

  if (HessianPresent)
    saveHessian(file);
  else
    saveSimple(file);
}

/** Saves all atoms and data into a MPQC config file without hessian.
 * \param *file output stream
 */
void MpqcParser::saveSimple(ostream *file)
{
  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;
    *file << "\tmole<MBPT2>: (" << endl;
    *file << "\t\tmaxiter = 200" << endl;
    *file << "\t\tbasis = $:basis" << endl;
    *file << "\t\tmolecule = $:molecule" << endl;
    *file << "\t\treference<CLHF>: (" << endl;
    *file << "\t\t\tbasis = $:basis" << endl;
    *file << "\t\t\tmolecule = $:molecule" << endl;
    *file << "\t\t)" << endl;
    *file << "\t)" << endl;
    *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;
  }
}

/** Saves all atoms and data into a MPQC config file with hessian.
 * \param *file output stream
 */
void MpqcParser::saveHessian(ostream *file)
{
  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());

  // with 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;
    *file << "\tmole<CLHF>: (" << endl;
    *file << "\t\tmaxiter = 200" << endl;
    *file << "\t\tbasis = $:basis" << endl;
    *file << "\t\tmolecule = $:molecule" << endl;
    *file << "\t)" << endl;
    *file << "\tfreq<MolecularFrequencies>: (" << endl;
    *file << "\t\tmolecule=$:molecule" << endl;
    *file << "\t)" << endl;
    *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;
  }
}

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