/* * 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 #endif #include #include #include #include "CodePatterns/MemDebug.hpp" #include "MpqcParser.hpp" #include "MpqcParser_Parameters.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() { parameters = new MpqcParser_Parameters(); } /** Destructor of MpqcParser. * */ MpqcParser::~MpqcParser() { delete parameters; } /** Load an MPQC config file into the World. * \param *file input stream */ void MpqcParser::load(istream *file) { bool MpqcSection = false; bool MoleculeSection = false; bool GeometrySection = false; bool BasisSection = false; bool AuxBasisSection = false; char line[MAXSTRINGSIZE]; typedef boost::tokenizer > tokenizer; boost::char_separator sep("[]"); boost::char_separator angularsep("<>"); boost::char_separator equalitysep(" ="); boost::char_separator whitesep(" \t"); ConvertTo 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(); ASSERT(tok_iter != tokens.end(), "MpqcParser::load() - missing token for MoleculeSection in line "+linestring+"!"); std::stringstream whitespacefilter(*tok_iter++); std::string element; whitespacefilter >> ws >> element; ASSERT(tok_iter != tokens.end(), "MpqcParser::load() - missing token for MoleculeSection in line "+linestring+"!"); 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; isetType(World::getInstance().getPeriode()->FindElement(element)); newAtom->setPosition(X); newmol->AddAtom(newAtom); DoLog(1) && (Log() << Verbose(1) << "Adding atom " << *newAtom << std::endl); } } 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(), "MpqcParser::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(), "MpqcParser::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(), "MpqcParser::load() - missing token before '=' for MpqcSection in line "+linestring+"!"); std::stringstream whitespacefilter(*tok_iter); std::string key; whitespacefilter >> ws >> key; if (getParams().haveParam(key)) { std::stringstream linestream(linestring); linestream >> getParams(); } else { // unknown keys are simply ignored as long as parser is incomplete DoLog(2) && (Log() << Verbose(2) << "INFO: '"+key+"' is unknown and ignored." << std::endl); } } } if (BasisSection) { tokenizer tokens(linestring, equalitysep); tokenizer::iterator tok_iter = tokens.begin(); ASSERT(tok_iter != tokens.end(), "MpqcParser::load() - missing token for BasisSection in line "+linestring+"!"); std::string key(*tok_iter++); ASSERT(tok_iter != tokens.end(), "MpqcParser::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(), "MpqcParser::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(), "MpqcParser::load() - missing token for AuxBasisSection in line "+linestring+"!"); std::string key(*tok_iter++); ASSERT(tok_iter != tokens.end(), "MpqcParser::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(), "MpqcParser::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:") != string::npos) { MoleculeSection = true; } if (linestring.find("atoms geometry") != string::npos) { GeometrySection = true; } if ((linestring.find("basis:") != 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(); } /** 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 &atoms) { Vector center; // vector allatoms = World::getInstance().getAllAtoms(); // calculate center for (std::vector::const_iterator runner = atoms.begin();runner != atoms.end(); ++runner) center += (*runner)->getPosition(); center.Scale(1./(double)atoms.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 = " << getParams().getString(MpqcParser_Parameters::savestateParam) << endl; *file << "\tdo_gradient = " << getParams().getString(MpqcParser_Parameters::do_gradientParam) << endl; if (getParams().getBool(MpqcParser_Parameters::hessianParam)) { *file << "\tfreq: (" << endl; *file << "\t\tmolecule=$:molecule" << endl; *file << "\t)" << endl; } switch (getParams().getTheory()) { case MpqcParser_Parameters::CLHF: *file << "\tmole<" << getParams().getString(MpqcParser_Parameters::theoryParam) << ">: (" << endl; *file << "\t\tmolecule = $:molecule" << endl; *file << "\t\tbasis = $:basis" << endl; *file << "\t\tmaxiter = " << toString(getParams().getInt(MpqcParser_Parameters::maxiterParam))<< endl; *file << "\t\tmemory = " << toString(getParams().getInt(MpqcParser_Parameters::memoryParam)) << endl; *file << "\t)" << endl; break; case MpqcParser_Parameters::CLKS: *file << "\tmole<" << getParams().getString(MpqcParser_Parameters::theoryParam) << ">: (" << endl; *file << "\t\tfunctional:(name=B3LYP)" << endl; *file << "\t\tmolecule = $:molecule" << endl; *file << "\t\tbasis = $:basis" << endl; *file << "\t\tmaxiter = " << toString(getParams().getInt(MpqcParser_Parameters::maxiterParam))<< endl; *file << "\t\tmemory = " << toString(getParams().getInt(MpqcParser_Parameters::memoryParam)) << endl; *file << "\t)" << endl; break; case MpqcParser_Parameters::MBPT2: *file << "\tmole<" << getParams().getString(MpqcParser_Parameters::theoryParam) << ">: (" << endl; *file << "\t\tbasis = $:basis" << endl; *file << "\t\tmolecule = $:molecule" << endl; *file << "\t\tmemory = " << toString(getParams().getInt(MpqcParser_Parameters::memoryParam)) << endl; *file << "\t\treference: (" << endl; *file << "\t\t\tmaxiter = " << toString(getParams().getInt(MpqcParser_Parameters::maxiterParam))<< endl; *file << "\t\t\tbasis = $:basis" << endl; *file << "\t\t\tmolecule = $:molecule" << endl; *file << "\t\t\tmemory = " << toString(getParams().getInt(MpqcParser_Parameters::memoryParam)) << endl; *file << "\t\t)" << endl; *file << "\t)" << endl; break; case MpqcParser_Parameters::MBPT2_R12: *file << "\tmole<" << getParams().getString(MpqcParser_Parameters::theoryParam) << ">: (" << endl; *file << "\t\tmolecule = $:molecule" << endl; *file << "\t\tbasis = $:basis" << endl; *file << "\t\taux_basis = $:abasis" << endl; *file << "\t\tstdapprox = \"" << getParams().getString(MpqcParser_Parameters::stdapproxParam) << "\"" << endl; *file << "\t\tnfzc = " << toString(getParams().getInt(MpqcParser_Parameters::nfzcParam)) << endl; *file << "\t\tmemory = " << toString(getParams().getInt(MpqcParser_Parameters::memoryParam)) << endl; *file << "\t\tintegrals:()" << endl; *file << "\t\treference: (" << endl; *file << "\t\t\tmolecule = $:molecule" << endl; *file << "\t\t\tbasis = $:basis" << endl; *file << "\t\t\tmaxiter = " << toString(getParams().getInt(MpqcParser_Parameters::maxiterParam)) << endl; *file << "\t\t\tmemory = " << toString(getParams().getInt(MpqcParser_Parameters::memoryParam)) << endl; *file << "\t\t\tintegrals<" << getParams().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: (" << endl; *file << "\tunit = " << (World::getInstance().getConfig()->GetIsAngstroem() ? "angstrom" : "bohr" ) << endl; *file << "\t{ atoms geometry } = {" << endl; // output of atoms for (std::vector::const_iterator AtomRunner = atoms.begin(); AtomRunner != atoms.end(); ++AtomRunner) { (*AtomRunner)->OutputMPQCLine(file, ¢er); } *file << "\t}" << endl; *file << ")" << endl; *file << "basis: (" << endl; *file << "\tname = \"" << getParams().getString(MpqcParser_Parameters::basisParam) << "\"" << endl; *file << "\tmolecule = $:molecule" << endl; *file << ")" << endl; if (getParams().getTheory() == MpqcParser_Parameters::MBPT2_R12) { *file << "% auxiliary basis set specification" << endl; *file << "\tabasis: (" << endl; *file << "\tname = \"" << getParams().getString(MpqcParser_Parameters::aux_basisParam) << "\"" << endl; *file << "\tmolecule = $:molecule" << endl; *file << ")" << endl; } } }