source: src/Parser/Psi3Parser.cpp@ 440ac3

/*
 * 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/>.
 */

/*
 * Psi3Parser.cpp
 *
 *  Created on: Oct 04, 2011
 *      Author: heber
 */

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

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

#include "CodePatterns/MemDebug.hpp"

#include "Psi3Parser.hpp"
#include "Psi3Parser_Parameters.hpp"

#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/toString.hpp"
#include "Element/element.hpp"
#include "Element/periodentafel.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "molecule.hpp"
#include "World.hpp"

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

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

/** Constructor of Psi3Parser.
 *
 */
FormatParser< psi3 >::FormatParser() :
  FormatParser_common(new Psi3Parser_Parameters())
{}

/** Destructor of Psi3Parser.
 *
 */
FormatParser< psi3 >::~FormatParser()
{}

/** Load an PSI3 config file into the World.
 * \param *file input stream
 */
void FormatParser< psi3 >::load(istream *file)
{
  bool Psi3Section = false;
  bool GeometrySection = 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;

  molecule *newmol = World::getInstance().createMolecule();
  newmol->ActiveFlag = true;
  while (file->good()) {
    file->getline(line, MAXSTRINGSIZE-1);
    std::string linestring(line);
    LOG(3, "INFO: Current line is: " << line);
    if ((linestring.find(")") != string::npos) && (linestring.find("(") == string::npos)) {
      LOG(3, "INFO: Line contains ')' and no '(' (end of section): " << line);
      // ends a section which do not overlap
      if (GeometrySection)
        GeometrySection = false;
      else
        Psi3Section = false;
    }
    if (GeometrySection) { // we have an atom
      tokenizer tokens(linestring, sep);
//      if (tokens.size() != 2)
//        throw Psi3ParseException;
      tokenizer::iterator tok_iter = tokens.begin();
      ASSERT(tok_iter != tokens.end(),
          "FormatParser< psi3 >::load() - missing token for MoleculeSection in line "+linestring+"!");
      std::stringstream whitespacefilter(*++tok_iter);
      std::string element;
      whitespacefilter >> ws >> element;
      LOG(2, "INFO: element of atom is " << element);
      ASSERT(tok_iter != tokens.end(),
          "FormatParser< psi3 >::load() - missing token for MoleculeSection in line "+linestring+"!");
      std::string vector = *tok_iter;
      tokenizer vectorcomponents(vector, whitesep);
      Vector X;
//      if (vectorcomponents.size() != NDIM)
//        throw Psi3ParseException;
      tok_iter = vectorcomponents.begin();
      ++tok_iter;
      for (int i=0; i<NDIM; ++i) {
        LOG(4, "INFO: Current value is " << *tok_iter << ".");
        X[i] = toDouble(*tok_iter++);
      }
      LOG(2, "INFO: position of atom is " << X);
      // 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 ((Psi3Section) && (!GeometrySection)) {
      if (linestring.find("=") != string::npos) { // get param value
        tokenizer tokens(linestring, equalitysep);
        tokenizer::iterator tok_iter = tokens.begin();
        ASSERT(tok_iter != tokens.end(),
            "FormatParser< psi3 >::load() - missing token before '=' for Psi3Section in line "+linestring+"!");
        std::stringstream whitespacefilter(*tok_iter);
        std::string key;
        whitespacefilter >> ws >> key;
        //LOG(2, "INFO: key to check is: " << key);
        if (getParams().haveParameter(key)) {
          //LOG(2, "INFO: Line submitted to parameter is: " << linestring);
          std::stringstream linestream(linestring);
          linestream >> getParams();
        } else { // unknown keys are simply ignored as long as parser is incomplete
          LOG(3, "INFO: '"+key+"' is unknown and ignored.");
        }
      }
    }
    if ((linestring.find("geometry") != string::npos) && (linestring.find("(") != string::npos)) {
      LOG(3, "INFO: Line contains geometry and '(': " << line);
      GeometrySection = true;
    }
    if ((linestring.find("psi:") != string::npos) && (linestring.find("(") != string::npos)) {
      LOG(3, "INFO: Line contains psi: and '(': " << line);
      Psi3Section = true;
    }
  }
  // refresh atom::nr and atom::name
  newmol->getAtomCount();
}

void FormatParser< psi3 >::OutputPsi3Line(ostream * const out, const atom &_atom, const Vector *center) const
{
  Vector recentered(_atom.getPosition());
  recentered -= *center;
  *out << "\t( " << _atom.getType()->getSymbol() << "\t" << recentered[0] << "\t" << recentered[1] << "\t" << recentered[2] << " )" << endl;
};


/** Saves all atoms and data into a PSI3 config file.
 * \param *file output stream
 * \param atoms atoms to store
 */
void FormatParser< psi3 >::save(ostream *file, const std::vector<const atom *> &atoms)
{
  Vector center;
//  vector<atom *> allatoms = 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());

  // first without hessian
  if (file->fail()) {
    ELOG(1, "Cannot open psi3 output file.");
  } else {
    *file << "% Created by MoleCuilder" << std::endl;
    *file << "psi: (" << std::endl;
    *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::labelParam)
                  << " = \"" << getParams().getParameter(Psi3Parser_Parameters::labelParam) << "\"" << std::endl;
    *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::jobtypeParam)
                  << " = " << getParams().getParameter(Psi3Parser_Parameters::jobtypeParam) << std::endl;
    *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::wavefunctionParam)
                  << " = " << getParams().getParameter(Psi3Parser_Parameters::wavefunctionParam) << std::endl;
    *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::maxiterParam)
                  << " = " << getParams().getParameter(Psi3Parser_Parameters::maxiterParam) << std::endl;
    *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::referenceParam)
                  << " = " << getParams().getParameter(Psi3Parser_Parameters::referenceParam) << std::endl;
    *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::basisParam)
                  << " = \"" << getParams().getParameter(Psi3Parser_Parameters::basisParam) << "\"" << std::endl;
    const std::string reference = getParams().getParameter(Psi3Parser_Parameters::referenceParam);
//    if (reference == getParams().getReferenceName(Psi3Parser_Parameters::RHF)) {
//    }
//    if (reference == getParams().getReferenceName(Psi3Parser_Parameters::ROHF)) {
//    }
    if ((reference == getParams().getReferenceName(Psi3Parser_Parameters::UHF))
        || (reference == getParams().getReferenceName(Psi3Parser_Parameters::TWOCON))) {
      const unsigned int multiplicity = calculateMultiplicity(atoms);
      *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::multiplicityParam)
                    << " = " << multiplicity << std::endl;
//                    << " = " << getParams().getParameter(Psi3Parser_Parameters::multiplicityParam) << std::endl;
      *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::chargeParam)
                    << " = " << getParams().getParameter(Psi3Parser_Parameters::chargeParam) << std::endl;
    }
    if (reference == getParams().getReferenceName(Psi3Parser_Parameters::TWOCON)) {
      *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::soccParam)
                    << " = " << getParams().getParameter(Psi3Parser_Parameters::soccParam) << std::endl;
      *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::doccParam)
                    << " = " << getParams().getParameter(Psi3Parser_Parameters::doccParam) << std::endl;
      *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::subgroupParam)
                    << " = " << getParams().getParameter(Psi3Parser_Parameters::subgroupParam) << std::endl;
      *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::unique_axisParam)
                    << " = " << getParams().getParameter(Psi3Parser_Parameters::unique_axisParam) << std::endl;
    }
    if ((reference != getParams().getReferenceName(Psi3Parser_Parameters::RHF))
        && (reference != getParams().getReferenceName(Psi3Parser_Parameters::ROHF))
        && (reference != getParams().getReferenceName(Psi3Parser_Parameters::UHF))
        && (reference != getParams().getReferenceName(Psi3Parser_Parameters::TWOCON)))
    {
      ELOG(0, "Unknown level of reference requested for Psi3 output file.");
    }
    *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::freeze_coreParam)
                  << " = " << getParams().getParameter(Psi3Parser_Parameters::freeze_coreParam) << std::endl;
    *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::unitsParam)
                  << " = " << getParams().getParameter(Psi3Parser_Parameters::unitsParam) << std::endl;
    *file << "\tgeometry = (" << std::endl;
    // output of atoms
    BOOST_FOREACH(const atom *_atom, atoms) {
      OutputPsi3Line(file, *_atom, &center);
    }
    *file << "\t)" << std::endl;
    *file << "\t" << getParams().getParameterName(Psi3Parser_Parameters::originParam)
                  << " = " << getParams().getParameter(Psi3Parser_Parameters::originParam) << std::endl;
    *file << ")" << std::endl;
  }
}

unsigned int FormatParser< psi3 >::calculateMultiplicity(
    const std::vector<const atom *> &atoms) const
{
  // add up the number of electrons
  double valencies = 0.;
  BOOST_FOREACH(const atom *_atom, atoms) {
    valencies += _atom->getType()->getAtomicNumber();
  }

  // add doubly up all bond degrees (two electrons per bond)
  unsigned int degrees = 0;
  BOOST_FOREACH(const atom *_atom, atoms) {
    BOOST_FOREACH(bond::ptr _bond, _atom->getListOfBonds()) {
      degrees += 2*_bond->getDegree();
    }
  }

  // return difference+1 as multiplicity
  return (int)fabs(valencies-degrees)+1;
}