Changes in src/Parser/MpqcParser.cpp [ad011c:4cbca0]

File:

• src/Parser/MpqcParser.cpp (modified) (7 diffs)

src/Parser/MpqcParser.cpp

@@ -18 +18 @@
 #endif
 
+#include <iostream>
+#include <boost/tokenizer.hpp>
+#include <string>
+
 #include "CodePatterns/MemDebug.hpp"
 
@@ -25 +29 @@
 #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"
@@ -35 +41 @@
  *
  */
-MpqcParser::MpqcParser() : HessianPresent(false)
+MpqcParser::MpqcParser()
 {}
 
@@ -49 +55 @@
 void MpqcParser::load(istream *file)
 {
-  // TODO: MpqcParser::load implementation
-  ASSERT(false, "Not implemented yet");
+  bool MpqcSection = false;
+  bool MoleculeSection = false;
+  bool GeometrySection = 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;
+
+  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 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
+        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);
+      }
+    }
+    if (MpqcSection) {
+      tokenizer tokens(linestring, equalitysep);
+      tokenizer::iterator tok_iter = tokens.begin();
+      std::stringstream whitespacefilter(*tok_iter++);
+      std::string key(*tok_iter++);
+      if (params.haveParam(key)) {
+        std::stringstream linestream(linestring);
+        linestream >> params;
+      } else if (linestring.find("mole<") != string::npos) { // get theory
+        tokenizer tokens(linestring, angularsep);
+        tokenizer::iterator tok_iter = tokens.begin();
+        std::string value(*(++tok_iter));
+        std::stringstream linestream("theory = "+value);
+        linestream >> params;
+      } else if (linestring.find("integrals<") != string::npos) { // get theory
+        tokenizer tokens(linestring, angularsep);
+        tokenizer::iterator tok_iter = tokens.begin();
+        std::string value(*(++tok_iter));
+        std::stringstream linestream("integration = "+value);
+        linestream >> params;
+      }
+    }
+    if (BasisSection) {
+      tokenizer tokens(linestring, equalitysep);
+      tokenizer::iterator tok_iter = tokens.begin();
+      std::string key(*tok_iter++);
+      std::string value(*tok_iter);
+      // TODO: use exception instead of ASSERT
+      ASSERT(tok_iter == tokens.end(),
+          "MpqcParser::load() - more than (key = value) on line "+linestring+".");
+      if (key == "name") {
+        std::stringstream linestream("basis = "+value);
+        linestream >> params;
+      }
+    }
+    if (AuxBasisSection) {
+      tokenizer tokens(linestring, equalitysep);
+      tokenizer::iterator tok_iter = tokens.begin();
+      std::string key(*tok_iter++);
+      std::string value(*tok_iter);
+      // TODO: use exception instead of ASSERT
+      ASSERT(tok_iter == tokens.end(),
+          "MpqcParser::load() - more than (key = value) on line "+linestring+".");
+      if (key == "name") {
+        std::stringstream linestream("aux_basis = "+value);
+        linestream >> params;
+      }
+    }
+    // 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 geometry") != string::npos) {
+      GeometrySection = true;
+    }
+    if ((linestring.find("basis<GaussianBasisSet>:") != string::npos) && ((linestring.find("abasis<") == string::npos))) {
+      BasisSection = true;
+    }
+    if (linestring.find("abasis<") != string::npos) {
+      AuxBasisSection = true;
+    }
+  }
 }
 
-/**
- * Saves the \a atoms into as a MPQC file.
- *
- * \param file where to save the state
- * \param atoms atoms to store
- */
-void MpqcParser::save(ostream *file, const std::vector<atom *> &atoms)
-{
-  DoLog(0) && (Log() << Verbose(0) << "Saving changes to MPQC ." << std::endl);
-
-  if (HessianPresent)
-    saveHessian(file, atoms);
-  else
-    saveSimple(file, atoms);
-}
-
-/** Saves all atoms and data into a MPQC config file without hessian.
+/** Saves all atoms and data into a MPQC config file.
  * \param *file output stream
  * \param atoms atoms to store
  */
-void MpqcParser::saveSimple(ostream *file, const std::vector<atom *> &atoms)
+void MpqcParser::save(ostream *file, const std::vector<atom *> &atoms)
 {
   Vector center;
@@ -81 +189 @@
   for (vector<atom *>::iterator runner = allatoms.begin();runner != allatoms.end(); ++runner)
     center += (*runner)->getPosition();
-  center.Scale(1./allatoms.size());
+  center.Scale(1./(double)allatoms.size());
 
   // first without hessian
@@ -89 +197 @@
     *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 << "\tsavestate = " << params.getString(MpqcParser_Parameters::savestateParam) << endl;
+    *file << "\tdo_gradient = " << params.getString(MpqcParser_Parameters::do_gradientParam) << endl;
+    if (params.getBool(MpqcParser_Parameters::hessianParam)) {
+      *file << "\tfreq<MolecularFrequencies>: (" << endl;
+      *file << "\t\tmolecule=$:molecule" << endl;
+      *file << "\t)" << endl;
+    }
+    switch (params.getTheory()) {
+      case MpqcParser_Parameters::CLHF:
+        *file << "\tmole<" << params.getString(MpqcParser_Parameters::theoryParam) << ">: (" << endl;
+        *file << "\t\tmolecule = $:molecule" << endl;
+        *file << "\t\tbasis = $:basis" << endl;
+        *file << "\t\tmaxiter = " << toString(params.getInt(MpqcParser_Parameters::maxiterParam))<< endl;
+        *file << "\t\tmemory = " << toString(params.getInt(MpqcParser_Parameters::memoryParam)) << endl;
+        *file << "\t)" << endl;
+        break;
+      case MpqcParser_Parameters::CLKS:
+        *file << "\tmole<" << params.getString(MpqcParser_Parameters::theoryParam) << ">: (" << endl;
+        *file << "\t\tfunctional<StdDenFunctional>:(name=B3LYP)" << endl;
+        *file << "\t\tmolecule = $:molecule" << endl;
+        *file << "\t\tbasis = $:basis" << endl;
+        *file << "\t\tmaxiter = " << toString(params.getInt(MpqcParser_Parameters::maxiterParam))<< endl;
+        *file << "\t\tmemory = " << toString(params.getInt(MpqcParser_Parameters::memoryParam)) << endl;
+        *file << "\t)" << endl;
+        break;
+      case MpqcParser_Parameters::MBPT2:
+        *file << "\tmole<" << params.getString(MpqcParser_Parameters::theoryParam) << ">: (" << endl;
+        *file << "\t\tbasis = $:basis" << endl;
+        *file << "\t\tmolecule = $:molecule" << endl;
+        *file << "\t\tmemory = " << toString(params.getInt(MpqcParser_Parameters::memoryParam)) << endl;
+        *file << "\t\treference<CLHF>: (" << endl;
+        *file << "\t\t\tmaxiter = " << toString(params.getInt(MpqcParser_Parameters::maxiterParam))<< endl;
+        *file << "\t\t\tbasis = $:basis" << endl;
+        *file << "\t\t\tmolecule = $:molecule" << endl;
+        *file << "\t\t\tmemory = " << toString(params.getInt(MpqcParser_Parameters::memoryParam)) << endl;
+        *file << "\t\t)" << endl;
+        *file << "\t)" << endl;
+        break;
+      case MpqcParser_Parameters::MBPT2_R12:
+        *file << "\tmole<" << params.getString(MpqcParser_Parameters::theoryParam) << ">: (" << endl;
+        *file << "\t\tmolecule = $:molecule" << endl;
+        *file << "\t\tbasis = $:basis" << endl;
+        *file << "\t\taux_basis = $:abasis" << endl;
+        *file << "\t\tstdapprox = \"" << params.getString(MpqcParser_Parameters::stdapproxParam) << "\"" << endl;
+        *file << "\t\tnfzc = " << toString(params.getInt(MpqcParser_Parameters::nfzcParam)) << endl;
+        *file << "\t\tmemory = " << toString(params.getInt(MpqcParser_Parameters::memoryParam)) << 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 = " << toString(params.getInt(MpqcParser_Parameters::maxiterParam)) << endl;
+        *file << "\t\t\tmemory = " << toString(params.getInt(MpqcParser_Parameters::memoryParam)) << endl;
+        *file << "\t\t\tintegrals<" << params.getString(MpqcParser_Parameters::integrationParam) << ">:()" << 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;
@@ -111 +269 @@
     *file << ")" << endl;
     *file << "basis<GaussianBasisSet>: (" << endl;
-    *file << "\tname = \"" << World::getInstance().getConfig()->basis << "\"" << endl;
+    *file << "\tname = \"" << params.getString(MpqcParser_Parameters::basisParam) << "\"" << endl;
     *file << "\tmolecule = $:molecule" << endl;
     *file << ")" << endl;
+    if (params.getTheory() == MpqcParser_Parameters::MBPT2_R12) {
+      *file << "% auxiliary basis set specification" << endl;
+      *file << "\tabasis<GaussianBasisSet>: (" << endl;
+      *file << "\tname = \"" << params.getString(MpqcParser_Parameters::aux_basisParam) << "\"" << endl;
+      *file << "\tmolecule = $:molecule" << endl;
+      *file << ")" << endl;
+    }
   }
 }
 
-/** Saves all atoms and data into a MPQC config file with hessian.
- * \param *file output stream
- * \param atoms atoms to store
- */
-void MpqcParser::saveHessian(ostream *file, const std::vector<atom *> &atoms)
+MpqcParser_Parameters & MpqcParser::getParams()
 {
-  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;
-  }
+  return params;
 }
 
-/** Sets whether hessian is desired or not
- * \param hessian statement
- */
-void MpqcParser::setHessian(bool hessian)
-{
-  HessianPresent = hessian;
-}