source: src/Parser/MpqcParser.cpp@ df5b8c

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * Copyright (C)  2013 Frederik Heber. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

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

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

#include <iomanip>
#include <iostream>
#include <boost/foreach.hpp>
#include <boost/tokenizer.hpp>
#include <sstream>
#include <string>

#include "CodePatterns/MemDebug.hpp"

#include "MpqcParser.hpp"
#include "MpqcParser_Parameters.hpp"

#include "Atom/atom.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/toString.hpp"
#include "config.hpp"
#include "Element/element.hpp"
#include "Element/periodentafel.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "molecule.hpp"
#include "MoleculeListClass.hpp"
#include "Parser/Exceptions.hpp"
#include "World.hpp"

// declare specialized static variables
const std::string FormatParserTrait<mpqc>::name = "mpqc";
const std::string FormatParserTrait<mpqc>::suffix = "in";
const ParserTypes FormatParserTrait<mpqc>::type = mpqc;

// a converter we often need
ConvertTo<bool> FormatParser<mpqc>::Converter;

/** Constructor of MpqcParser.
 *
 */
FormatParser< mpqc >::FormatParser() :
  FormatParser_common(new MpqcParser_Parameters())
{}

/** Destructor of MpqcParser.
 *
 */
FormatParser< mpqc >::~FormatParser()
{}

/** Load an MPQC config file into the World.
 * \param *file input stream
 */
void FormatParser< mpqc >::load(istream *file)
{
  bool MpqcSection = false;
  bool MoleculeSection = false;
  bool GeometrySection = false;
  bool GeometrySection_n = false;
  bool BasisSection = false;
  bool AuxBasisSection = false;
  char line[MAXSTRINGSIZE];
  typedef boost::tokenizer<boost::char_separator<char> > tokenizer;
  boost::char_separator<char> sep("[]");
  boost::char_separator<char> angularsep("<>");
  boost::char_separator<char> equalitysep(" =");
  boost::char_separator<char> whitesep(" \t");
  ConvertTo<double> toDouble;
  int old_n = -1; // note down the last parsed "n" to ascertain it's ascending

  molecule *newmol = World::getInstance().createMolecule();
  newmol->ActiveFlag = true;
  // TODO: Remove the insertion into molecule when saving does not depend on them anymore. Also, remove molecule.hpp include
  World::getInstance().getMolecules()->insert(newmol);
  while (file->good()) {
    file->getline(line, MAXSTRINGSIZE-1);
    std::string linestring(line);
    if (((linestring.find("atoms") == string::npos)
        || (linestring.find("geometry") == string::npos))
        && (linestring.find("}") != string::npos)) {
      GeometrySection = false;
    }
    if ((linestring.find(")") != string::npos)) { // ends a section which do not overlap
      MpqcSection = false;
      MoleculeSection = false;
      BasisSection = false;
      AuxBasisSection = false;
    }
    if (MoleculeSection) {
      if (GeometrySection) { // we have an atom
//        LOG(2, "DEBUG: Full line is '" << linestring << "'.");
        // separate off [..] part
        tokenizer tokens(linestring, sep);
        // split part prior to [..] into tokens
        std::string prior_part(*tokens.begin());
        tokenizer prefixtokens(prior_part, whitesep);
        tokenizer::iterator tok_iter = prefixtokens.begin();
//        LOG(2, "DEBUG: Current tok_iter is " << *tok_iter << ".");
        if (GeometrySection_n) {
          ASSERT(tok_iter != prefixtokens.end(),
              "FormatParser< mpqc >::load() - missing n entry for MoleculeSection in line "
              +linestring+"!");
          // if additional n is given, parse and check but discard eventually
          int n;
          std::stringstream whitespacefilter(*tok_iter++);
          whitespacefilter >> ws >> n;
          if ((old_n != -1) && ((n-1) != old_n))
            ELOG(2, "n index is not simply ascending by one but "
                << n-old_n << ", specific ids are lost!");
          if (old_n >= n) {
            ELOG(1, "n index is not simply ascending, coordinates might get mixed!");
            throw ParserException();
          }
          old_n = n;
        }
        ASSERT(tok_iter != prefixtokens.end(),
            "FormatParser< mpqc >::load() - missing atoms entry for MoleculeSection in line "
            +linestring+"!");
//        LOG(2, "DEBUG: Current tok_iter is " << *tok_iter << ".");
        std::string element;
        {
          std::stringstream whitespacefilter(*tok_iter++);
          whitespacefilter >> ws >> element;
        }
        // split [..] part and parse
        tok_iter = (++tokens.begin());
        ASSERT (tok_iter != tokens.end(),
            "FormatParser< mpqc >::load() - missing geometry entry for MoleculeSection in line "
            +linestring+"!");
//        LOG(2, "DEBUG: Current tok_iter is " << *tok_iter << ".");
        std::string vector(*tok_iter);
        tokenizer vectorcomponents(vector, whitesep);
        Vector X;
//        if (vectorcomponents.size() != NDIM)
//          throw ParserException();
        tok_iter = vectorcomponents.begin();
        for (int i=0; i<NDIM; ++i) {
//          LOG(2, "DEBUG: Current tok_iter is " << *tok_iter << ".");
          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);
        LOG(1, "Adding atom " << *newAtom);
      }
    }
    if (MpqcSection) {
      if (linestring.find("mole<") != string::npos) { // get theory
        tokenizer tokens(linestring, angularsep);
        tokenizer::iterator tok_iter = tokens.begin();
        ++tok_iter;
        ASSERT(tok_iter != tokens.end(),
            "FormatParser< mpqc >::load() - missing token in brackets<> for mole< in line "+linestring+"!");
        std::string value(*tok_iter);
        std::stringstream linestream("theory = "+value);
        linestream >> getParams();
      } else if (linestring.find("integrals<") != string::npos) { // get theory
        tokenizer tokens(linestring, angularsep);
        tokenizer::iterator tok_iter = tokens.begin();
        ++tok_iter;
        ASSERT(tok_iter != tokens.end(),
            "FormatParser< mpqc >::load() - missing token in brackets<> for integrals< in line "+linestring+"!");
        std::string value(*tok_iter);
        std::stringstream linestream("integration = "+value);
        linestream >> getParams();
      } else if ((linestring.find("molecule") == string::npos) && (linestring.find("basis") == string::npos)){
        // molecule and basis must not be parsed in this section
        tokenizer tokens(linestring, equalitysep);
        tokenizer::iterator tok_iter = tokens.begin();
        ASSERT(tok_iter != tokens.end(),
            "FormatParser< mpqc >::load() - missing token before '=' for MpqcSection in line "+linestring+"!");
        std::stringstream whitespacefilter(*tok_iter);
        std::string key;
        whitespacefilter >> ws >> key;
        if (getParams().haveParameter(key)) {
          std::stringstream linestream(linestring);
          linestream >> getParams();
        } else { // unknown keys are simply ignored as long as parser is incomplete
          LOG(2, "INFO: '"+key+"' is unknown and ignored.");
        }
      }
    }
    if (BasisSection) {
      tokenizer tokens(linestring, equalitysep);
      tokenizer::iterator tok_iter = tokens.begin();
      ASSERT(tok_iter != tokens.end(),
          "FormatParser< mpqc >::load() - missing token for BasisSection in line "+linestring+"!");
      std::string key(*tok_iter++);
      ASSERT(tok_iter != tokens.end(),
          "FormatParser< mpqc >::load() - missing value for BasisSection after key "+key+" in line "+linestring+"!");
      std::string value(*tok_iter);
      tok_iter++;
      // TODO: use exception instead of ASSERT
      ASSERT(tok_iter == tokens.end(),
          "FormatParser< mpqc >::load() - more than (key = value) on line "+linestring+".");
      if (key == "name") {
        std::stringstream linestream("basis = "+value);
        linestream >> getParams();
      }
    }
    if (AuxBasisSection) {
      tokenizer tokens(linestring, equalitysep);
      tokenizer::iterator tok_iter = tokens.begin();
      ASSERT(tok_iter != tokens.end(),
          "FormatParser< mpqc >::load() - missing token for AuxBasisSection in line "+linestring+"!");
      std::string key(*tok_iter++);
      ASSERT(tok_iter != tokens.end(),
          "FormatParser< mpqc >::load() - missing value for BasisSection after key "+key+" in line "+linestring+"!");
      std::string value(*tok_iter);
      tok_iter++;
      // TODO: use exception instead of ASSERT
      ASSERT(tok_iter == tokens.end(),
          "FormatParser< mpqc >::load() - more than (key = value) on line "+linestring+".");
      if (key == "name") {
        std::stringstream linestream("aux_basis = "+value);
        linestream >> getParams();
      }
    }
    // set some scan flags
    if (linestring.find("mpqc:") != string::npos) {
      MpqcSection = true;
    }
    if (linestring.find("molecule<Molecule>:") != string::npos) {
      MoleculeSection = true;
    }
    if ((linestring.find("atoms") != string::npos)
        && (linestring.find("geometry") != string::npos)) {
      GeometrySection = true;
      if (linestring.find("n") != string::npos) // neither atoms nor geometry contains a letter n
        GeometrySection_n = true;
    }
    if ((linestring.find("basis<GaussianBasisSet>:") != string::npos) && ((linestring.find("abasis<") == string::npos))) {
      BasisSection = true;
    }
    if (linestring.find("abasis<") != string::npos) {
      AuxBasisSection = true;
    }
  }
  // refresh atom::nr and atom::name
  newmol->getAtomCount();
}

void FormatParser< mpqc >::OutputMPQCLine(ostream * const out, const atom &_atom, const Vector *center) const
{
  Vector recentered(_atom.getPosition());
  recentered -= *center;
  *out << "\t\t" << _atom.getType()->getSymbol() << " [ ";
  {
    std::stringstream posstream;
    posstream << std::setprecision(12) << recentered[0] << "\t" << recentered[1] << "\t" << recentered[2];
    *out << posstream.str();
  }
  *out << " ]" << std::endl;
};


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

  // calculate center
//  for (std::vector<const atom *>::const_iterator runner = atoms.begin();runner != atoms.end(); ++runner)
//    center += (*runner)->getPosition();
//  center.Scale(1./(double)atoms.size());
  center.Zero();

  // first without hessian
  if (file->fail()) {
    ELOG(1, "Cannot open mpqc output file.");
  } else {
    *file << "% Created by MoleCuilder" << endl;
    *file << "mpqc: (" << endl;
    *file << "\tcheckpoint = no" << endl;
    *file << "\tsavestate = " << getParams().getParameter(MpqcParser_Parameters::savestateParam) << endl;
    *file << "\tdo_gradient = " << getParams().getParameter(MpqcParser_Parameters::do_gradientParam) << endl;
    if (Converter(getParams().getParameter(MpqcParser_Parameters::hessianParam))) {
      *file << "\tfreq<MolecularFrequencies>: (" << endl;
      *file << "\t\tmolecule=$:molecule" << endl;
      *file << "\t)" << endl;
    }
    const std::string theory = getParams().getParameter(MpqcParser_Parameters::theoryParam);
    if (theory == getParams().getTheoryName(MpqcParser_Parameters::CLHF)) {
      *file << "\tmole<" << getParams().getParameter(MpqcParser_Parameters::theoryParam) << ">: (" << endl;
      *file << "\t\tmolecule = $:molecule" << endl;
      *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::basisParam) << " = $:basis" << endl;
      *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::maxiterParam)
          << " = " << getParams().getParameter(MpqcParser_Parameters::maxiterParam)<< endl;
      *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::memoryParam)
          << " = " << getParams().getParameter(MpqcParser_Parameters::memoryParam) << endl;
      *file << "\t)" << endl;
    } else if (theory == getParams().getTheoryName(MpqcParser_Parameters::CLKS)) {
      *file << "\tmole<" << getParams().getParameter(MpqcParser_Parameters::theoryParam) << ">: (" << endl;
      *file << "\t\tfunctional<StdDenFunctional>:(name=B3LYP)" << endl;
      *file << "\t\tmolecule = $:molecule" << endl;
      *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::basisParam) << " = $:basis" << endl;
      *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::maxiterParam)
          << " = " << getParams().getParameter(MpqcParser_Parameters::maxiterParam)<< endl;
      *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::memoryParam)
          << " = " << getParams().getParameter(MpqcParser_Parameters::memoryParam) << endl;
      *file << "\t)" << endl;
    } else if (theory == getParams().getTheoryName(MpqcParser_Parameters::MBPT2)) {
      *file << "\tmole<" << getParams().getParameter(MpqcParser_Parameters::theoryParam) << ">: (" << endl;
      *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::basisParam) << " = $:basis" << endl;
      *file << "\t\tmolecule = $:molecule" << endl;
      *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::memoryParam)
          << " = " << getParams().getParameter(MpqcParser_Parameters::memoryParam) << endl;
      *file << "\t\treference<CLHF>: (" << endl;
      *file << "\t\t\t" << getParams().getParameterName(MpqcParser_Parameters::maxiterParam)
          << " = " << getParams().getParameter(MpqcParser_Parameters::maxiterParam)<< endl;
      *file << "\t\t\t" << getParams().getParameterName(MpqcParser_Parameters::basisParam) << " = $:basis" << endl;
      *file << "\t\t\tmolecule = $:molecule" << endl;
      *file << "\t\t\t" << getParams().getParameterName(MpqcParser_Parameters::memoryParam)
          << " = " << getParams().getParameter(MpqcParser_Parameters::memoryParam) << endl;
      *file << "\t\t)" << endl;
      *file << "\t)" << endl;
    } else if (theory == getParams().getTheoryName(MpqcParser_Parameters::MBPT2_R12)) {
      *file << "\tmole<" << getParams().getParameter(MpqcParser_Parameters::theoryParam) << ">: (" << endl;
      *file << "\t\tmolecule = $:molecule" << endl;
      *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::basisParam) << " = $:basis" << endl;
      *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::aux_basisParam) << " = $:abasis" << endl;
      *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::stdapproxParam)
          << " = \"" << getParams().getParameter(MpqcParser_Parameters::stdapproxParam) << "\"" << endl;
      *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::nfzcParam)
          << " = " << getParams().getParameter(MpqcParser_Parameters::nfzcParam) << endl;
      *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::memoryParam)
          << " = " << getParams().getParameter(MpqcParser_Parameters::memoryParam) << endl;
      *file << "\t\tintegrals<IntegralCints>:()" << endl;
      *file << "\t\treference<CLHF>: (" << endl;
      *file << "\t\t\tmolecule = $:molecule" << endl;
      *file << "\t\t\t" << getParams().getParameterName(MpqcParser_Parameters::basisParam) << " = $:basis" << endl;
      *file << "\t\t\tmaxiter = " << getParams().getParameter(MpqcParser_Parameters::maxiterParam) << endl;
      *file << "\t\t\tmemory = " << getParams().getParameter(MpqcParser_Parameters::memoryParam) << endl;
      *file << "\t\t\tintegrals<" << getParams().getParameter(MpqcParser_Parameters::integrationParam) << ">:()" << endl;
      *file << "\t\t)" << endl;
      *file << "\t)" << endl;
    } else {
      ELOG(0, "Unknown level of theory requested for MPQC output file.");
    }
    const std::string jobtype = getParams().getParameter(MpqcParser_Parameters::jobtypeParam);
    if (jobtype == getParams().getJobtypeName(MpqcParser_Parameters::Optimization)) {
      *file << "\t% optimizer object for the molecular geometry" << endl;
      *file << "\topt<QNewtonOpt>: (" << endl;
      *file << "\t\tfunction = $..:mole" << endl;
      *file << "\t\tupdate<BFGSUpdate>: ()" << endl;
      *file << "\t\tconvergence<MolEnergyConvergence>: (" << endl;
      *file << "\t\t\tcartesian = yes" << endl;
      *file << "\t\t\tenergy = $..:..:mole" << 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
    BOOST_FOREACH(const atom *_atom, atoms) {
      OutputMPQCLine(file, *_atom, &center);
    }
    *file << "\t}" << endl;
    *file << ")" << endl;
    *file << "basis<GaussianBasisSet>: (" << endl;
    *file << "\tname = \"" << getParams().getParameter(MpqcParser_Parameters::basisParam) << "\"" << endl;
    *file << "\tmolecule = $:molecule" << endl;
    *file << ")" << endl;
    if (theory == getParams().getTheoryName(MpqcParser_Parameters::MBPT2_R12)) {
      *file << "% auxiliary basis set specification" << endl;
      *file << "\tabasis<GaussianBasisSet>: (" << endl;
      *file << "\tname = \"" << getParams().getParameter(MpqcParser_Parameters::aux_basisParam) << "\"" << endl;
      *file << "\tmolecule = $:molecule" << endl;
      *file << ")" << endl;
    }
  }
}