source: src/Parser/MpqcParser.cpp@ ff4fff9

CombiningParticlePotentialParsing
Last change on this file since ff4fff9 was 2affd1, checked in by Frederik Heber <heber@…>, 9 years ago

Removed molecules_deprecated from World and unnecessary includes of MoleculeListClass and all insert/erase.

  • this goes along the lines of removing "remove me when we don't need MoleculeCistClass anymore".
  • Property mode set to 100644
File size: 18.0 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
5 * Copyright (C) 2013 Frederik Heber. All rights reserved.
6 *
7 *
8 * This file is part of MoleCuilder.
9 *
10 * MoleCuilder is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 2 of the License, or
13 * (at your option) any later version.
14 *
15 * MoleCuilder is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
22 */
23
24/*
25 * MpqcParser.cpp
26 *
27 * Created on: 12.06.2010
28 * Author: heber
29 */
30
31// include config.h
32#ifdef HAVE_CONFIG_H
33#include <config.h>
34#endif
35
36#include <iomanip>
37#include <iostream>
38#include <boost/foreach.hpp>
39#include <boost/tokenizer.hpp>
40#include <sstream>
41#include <string>
42
43#include "CodePatterns/MemDebug.hpp"
44
45#include "MpqcParser.hpp"
46#include "MpqcParser_Parameters.hpp"
47
48#include "Atom/atom.hpp"
49#include "CodePatterns/Log.hpp"
50#include "CodePatterns/toString.hpp"
51#include "config.hpp"
52#include "Element/element.hpp"
53#include "Element/periodentafel.hpp"
54#include "LinearAlgebra/Vector.hpp"
55#include "molecule.hpp"
56#include "Parser/Exceptions.hpp"
57#include "World.hpp"
58
59// declare specialized static variables
60const std::string FormatParserTrait<mpqc>::name = "mpqc";
61const std::string FormatParserTrait<mpqc>::suffix = "in";
62const ParserTypes FormatParserTrait<mpqc>::type = mpqc;
63
64// a converter we often need
65ConvertTo<bool> FormatParser<mpqc>::Converter;
66
67/** Constructor of MpqcParser.
68 *
69 */
70FormatParser< mpqc >::FormatParser() :
71 FormatParser_common(new MpqcParser_Parameters())
72{}
73
74/** Destructor of MpqcParser.
75 *
76 */
77FormatParser< mpqc >::~FormatParser()
78{}
79
80/** Load an MPQC config file into the World.
81 * \param *file input stream
82 */
83void FormatParser< mpqc >::load(istream *file)
84{
85 bool MpqcSection = false;
86 bool MoleculeSection = false;
87 bool GeometrySection = false;
88 bool GeometrySection_n = false;
89 bool BasisSection = false;
90 bool AuxBasisSection = false;
91 char line[MAXSTRINGSIZE];
92 typedef boost::tokenizer<boost::char_separator<char> > tokenizer;
93 boost::char_separator<char> sep("[]");
94 boost::char_separator<char> angularsep("<>");
95 boost::char_separator<char> equalitysep(" =");
96 boost::char_separator<char> whitesep(" \t");
97 ConvertTo<double> toDouble;
98 int old_n = -1; // note down the last parsed "n" to ascertain it's ascending
99
100 molecule *newmol = World::getInstance().createMolecule();
101 newmol->ActiveFlag = true;
102 while (file->good()) {
103 file->getline(line, MAXSTRINGSIZE-1);
104 std::string linestring(line);
105 if (((linestring.find("atoms") == string::npos)
106 || (linestring.find("geometry") == string::npos))
107 && (linestring.find("}") != string::npos)) {
108 GeometrySection = false;
109 }
110 if ((linestring.find(")") != string::npos)) { // ends a section which do not overlap
111 MpqcSection = false;
112 MoleculeSection = false;
113 BasisSection = false;
114 AuxBasisSection = false;
115 }
116 if (MoleculeSection) {
117 if (GeometrySection) { // we have an atom
118// LOG(2, "DEBUG: Full line is '" << linestring << "'.");
119 // separate off [..] part
120 tokenizer tokens(linestring, sep);
121 // split part prior to [..] into tokens
122 std::string prior_part(*tokens.begin());
123 tokenizer prefixtokens(prior_part, whitesep);
124 tokenizer::iterator tok_iter = prefixtokens.begin();
125// LOG(2, "DEBUG: Current tok_iter is " << *tok_iter << ".");
126 if (GeometrySection_n) {
127 ASSERT(tok_iter != prefixtokens.end(),
128 "FormatParser< mpqc >::load() - missing n entry for MoleculeSection in line "
129 +linestring+"!");
130 // if additional n is given, parse and check but discard eventually
131 int n;
132 std::stringstream whitespacefilter(*tok_iter++);
133 whitespacefilter >> ws >> n;
134 if ((old_n != -1) && ((n-1) != old_n))
135 ELOG(2, "n index is not simply ascending by one but "
136 << n-old_n << ", specific ids are lost!");
137 if (old_n >= n) {
138 ELOG(1, "n index is not simply ascending, coordinates might get mixed!");
139 throw ParserException();
140 }
141 old_n = n;
142 }
143 ASSERT(tok_iter != prefixtokens.end(),
144 "FormatParser< mpqc >::load() - missing atoms entry for MoleculeSection in line "
145 +linestring+"!");
146// LOG(2, "DEBUG: Current tok_iter is " << *tok_iter << ".");
147 std::string element;
148 {
149 std::stringstream whitespacefilter(*tok_iter++);
150 whitespacefilter >> ws >> element;
151 }
152 // split [..] part and parse
153 tok_iter = (++tokens.begin());
154 ASSERT (tok_iter != tokens.end(),
155 "FormatParser< mpqc >::load() - missing geometry entry for MoleculeSection in line "
156 +linestring+"!");
157// LOG(2, "DEBUG: Current tok_iter is " << *tok_iter << ".");
158 std::string vector(*tok_iter);
159 tokenizer vectorcomponents(vector, whitesep);
160 Vector X;
161// if (vectorcomponents.size() != NDIM)
162// throw ParserException();
163 tok_iter = vectorcomponents.begin();
164 for (int i=0; i<NDIM; ++i) {
165// LOG(2, "DEBUG: Current tok_iter is " << *tok_iter << ".");
166 X[i] = toDouble(*tok_iter++);
167 }
168 // create atom
169 atom *newAtom = World::getInstance().createAtom();
170 newAtom->setType(World::getInstance().getPeriode()->FindElement(element));
171 newAtom->setPosition(X);
172 newmol->AddAtom(newAtom);
173 LOG(1, "Adding atom " << *newAtom);
174 }
175 }
176 if (MpqcSection) {
177 if (linestring.find("mole<") != string::npos) { // get theory
178 tokenizer tokens(linestring, angularsep);
179 tokenizer::iterator tok_iter = tokens.begin();
180 ++tok_iter;
181 ASSERT(tok_iter != tokens.end(),
182 "FormatParser< mpqc >::load() - missing token in brackets<> for mole< in line "+linestring+"!");
183 std::string value(*tok_iter);
184 std::stringstream linestream("theory = "+value);
185 linestream >> getParams();
186 } else if (linestring.find("integrals<") != string::npos) { // get theory
187 tokenizer tokens(linestring, angularsep);
188 tokenizer::iterator tok_iter = tokens.begin();
189 ++tok_iter;
190 ASSERT(tok_iter != tokens.end(),
191 "FormatParser< mpqc >::load() - missing token in brackets<> for integrals< in line "+linestring+"!");
192 std::string value(*tok_iter);
193 std::stringstream linestream("integration = "+value);
194 linestream >> getParams();
195 } else if ((linestring.find("molecule") == string::npos) && (linestring.find("basis") == string::npos)){
196 // molecule and basis must not be parsed in this section
197 tokenizer tokens(linestring, equalitysep);
198 tokenizer::iterator tok_iter = tokens.begin();
199 ASSERT(tok_iter != tokens.end(),
200 "FormatParser< mpqc >::load() - missing token before '=' for MpqcSection in line "+linestring+"!");
201 std::stringstream whitespacefilter(*tok_iter);
202 std::string key;
203 whitespacefilter >> ws >> key;
204 if (getParams().haveParameter(key)) {
205 std::stringstream linestream(linestring);
206 linestream >> getParams();
207 } else { // unknown keys are simply ignored as long as parser is incomplete
208 LOG(2, "INFO: '"+key+"' is unknown and ignored.");
209 }
210 }
211 }
212 if (BasisSection) {
213 tokenizer tokens(linestring, equalitysep);
214 tokenizer::iterator tok_iter = tokens.begin();
215 ASSERT(tok_iter != tokens.end(),
216 "FormatParser< mpqc >::load() - missing token for BasisSection in line "+linestring+"!");
217 std::string key(*tok_iter++);
218 ASSERT(tok_iter != tokens.end(),
219 "FormatParser< mpqc >::load() - missing value for BasisSection after key "+key+" in line "+linestring+"!");
220 std::string value(*tok_iter);
221 tok_iter++;
222 // TODO: use exception instead of ASSERT
223 ASSERT(tok_iter == tokens.end(),
224 "FormatParser< mpqc >::load() - more than (key = value) on line "+linestring+".");
225 if (key == "name") {
226 std::stringstream linestream("basis = "+value);
227 linestream >> getParams();
228 }
229 }
230 if (AuxBasisSection) {
231 tokenizer tokens(linestring, equalitysep);
232 tokenizer::iterator tok_iter = tokens.begin();
233 ASSERT(tok_iter != tokens.end(),
234 "FormatParser< mpqc >::load() - missing token for AuxBasisSection in line "+linestring+"!");
235 std::string key(*tok_iter++);
236 ASSERT(tok_iter != tokens.end(),
237 "FormatParser< mpqc >::load() - missing value for BasisSection after key "+key+" in line "+linestring+"!");
238 std::string value(*tok_iter);
239 tok_iter++;
240 // TODO: use exception instead of ASSERT
241 ASSERT(tok_iter == tokens.end(),
242 "FormatParser< mpqc >::load() - more than (key = value) on line "+linestring+".");
243 if (key == "name") {
244 std::stringstream linestream("aux_basis = "+value);
245 linestream >> getParams();
246 }
247 }
248 // set some scan flags
249 if (linestring.find("mpqc:") != string::npos) {
250 MpqcSection = true;
251 }
252 if (linestring.find("molecule<Molecule>:") != string::npos) {
253 MoleculeSection = true;
254 }
255 if ((linestring.find("atoms") != string::npos)
256 && (linestring.find("geometry") != string::npos)) {
257 GeometrySection = true;
258 if (linestring.find("n") != string::npos) // neither atoms nor geometry contains a letter n
259 GeometrySection_n = true;
260 }
261 if ((linestring.find("basis<GaussianBasisSet>:") != string::npos) && ((linestring.find("abasis<") == string::npos))) {
262 BasisSection = true;
263 }
264 if (linestring.find("abasis<") != string::npos) {
265 AuxBasisSection = true;
266 }
267 }
268 // refresh atom::nr and atom::name
269 newmol->getAtomCount();
270}
271
272void FormatParser< mpqc >::OutputMPQCLine(ostream * const out, const atom &_atom, const Vector *center) const
273{
274 Vector recentered(_atom.getPosition());
275 recentered -= *center;
276 *out << "\t\t" << _atom.getType()->getSymbol() << " [ ";
277 {
278 std::stringstream posstream;
279 posstream << std::setprecision(12) << recentered[0] << "\t" << recentered[1] << "\t" << recentered[2];
280 *out << posstream.str();
281 }
282 *out << " ]" << std::endl;
283};
284
285
286/** Saves all atoms and data into a MPQC config file.
287 * \param *file output stream
288 * \param atoms atoms to store
289 */
290void FormatParser< mpqc >::save(ostream *file, const std::vector<const atom *> &atoms)
291{
292 Vector center;
293// vector<const atom *> allatoms = const_cast<const World &>(World::getInstance()).
294// getAllAtoms();
295
296 // calculate center
297// for (std::vector<const atom *>::const_iterator runner = atoms.begin();runner != atoms.end(); ++runner)
298// center += (*runner)->getPosition();
299// center.Scale(1./(double)atoms.size());
300 center.Zero();
301
302 // first without hessian
303 if (file->fail()) {
304 ELOG(1, "Cannot open mpqc output file.");
305 } else {
306 *file << "% Created by MoleCuilder" << endl;
307 *file << "mpqc: (" << endl;
308 *file << "\tcheckpoint = no" << endl;
309 *file << "\tsavestate = " << getParams().getParameter(MpqcParser_Parameters::savestateParam) << endl;
310 *file << "\tdo_gradient = " << getParams().getParameter(MpqcParser_Parameters::do_gradientParam) << endl;
311 if (Converter(getParams().getParameter(MpqcParser_Parameters::hessianParam))) {
312 *file << "\tfreq<MolecularFrequencies>: (" << endl;
313 *file << "\t\tmolecule=$:molecule" << endl;
314 *file << "\t)" << endl;
315 }
316 const std::string theory = getParams().getParameter(MpqcParser_Parameters::theoryParam);
317 if (theory == getParams().getTheoryName(MpqcParser_Parameters::CLHF)) {
318 *file << "\tmole<" << getParams().getParameter(MpqcParser_Parameters::theoryParam) << ">: (" << endl;
319 *file << "\t\tmolecule = $:molecule" << endl;
320 *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::basisParam) << " = $:basis" << endl;
321 *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::maxiterParam)
322 << " = " << getParams().getParameter(MpqcParser_Parameters::maxiterParam)<< endl;
323 *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::memoryParam)
324 << " = " << getParams().getParameter(MpqcParser_Parameters::memoryParam) << endl;
325 *file << "\t)" << endl;
326 } else if (theory == getParams().getTheoryName(MpqcParser_Parameters::CLKS)) {
327 *file << "\tmole<" << getParams().getParameter(MpqcParser_Parameters::theoryParam) << ">: (" << endl;
328 *file << "\t\tfunctional<StdDenFunctional>:(name=B3LYP)" << endl;
329 *file << "\t\tmolecule = $:molecule" << endl;
330 *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::basisParam) << " = $:basis" << endl;
331 *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::maxiterParam)
332 << " = " << getParams().getParameter(MpqcParser_Parameters::maxiterParam)<< endl;
333 *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::memoryParam)
334 << " = " << getParams().getParameter(MpqcParser_Parameters::memoryParam) << endl;
335 *file << "\t)" << endl;
336 } else if (theory == getParams().getTheoryName(MpqcParser_Parameters::MBPT2)) {
337 *file << "\tmole<" << getParams().getParameter(MpqcParser_Parameters::theoryParam) << ">: (" << endl;
338 *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::basisParam) << " = $:basis" << endl;
339 *file << "\t\tmolecule = $:molecule" << endl;
340 *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::memoryParam)
341 << " = " << getParams().getParameter(MpqcParser_Parameters::memoryParam) << endl;
342 *file << "\t\treference<CLHF>: (" << endl;
343 *file << "\t\t\t" << getParams().getParameterName(MpqcParser_Parameters::maxiterParam)
344 << " = " << getParams().getParameter(MpqcParser_Parameters::maxiterParam)<< endl;
345 *file << "\t\t\t" << getParams().getParameterName(MpqcParser_Parameters::basisParam) << " = $:basis" << endl;
346 *file << "\t\t\tmolecule = $:molecule" << endl;
347 *file << "\t\t\t" << getParams().getParameterName(MpqcParser_Parameters::memoryParam)
348 << " = " << getParams().getParameter(MpqcParser_Parameters::memoryParam) << endl;
349 *file << "\t\t)" << endl;
350 *file << "\t)" << endl;
351 } else if (theory == getParams().getTheoryName(MpqcParser_Parameters::MBPT2_R12)) {
352 *file << "\tmole<" << getParams().getParameter(MpqcParser_Parameters::theoryParam) << ">: (" << endl;
353 *file << "\t\tmolecule = $:molecule" << endl;
354 *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::basisParam) << " = $:basis" << endl;
355 *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::aux_basisParam) << " = $:abasis" << endl;
356 *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::stdapproxParam)
357 << " = \"" << getParams().getParameter(MpqcParser_Parameters::stdapproxParam) << "\"" << endl;
358 *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::nfzcParam)
359 << " = " << getParams().getParameter(MpqcParser_Parameters::nfzcParam) << endl;
360 *file << "\t\t" << getParams().getParameterName(MpqcParser_Parameters::memoryParam)
361 << " = " << getParams().getParameter(MpqcParser_Parameters::memoryParam) << endl;
362 *file << "\t\tintegrals<IntegralCints>:()" << endl;
363 *file << "\t\treference<CLHF>: (" << endl;
364 *file << "\t\t\tmolecule = $:molecule" << endl;
365 *file << "\t\t\t" << getParams().getParameterName(MpqcParser_Parameters::basisParam) << " = $:basis" << endl;
366 *file << "\t\t\tmaxiter = " << getParams().getParameter(MpqcParser_Parameters::maxiterParam) << endl;
367 *file << "\t\t\tmemory = " << getParams().getParameter(MpqcParser_Parameters::memoryParam) << endl;
368 *file << "\t\t\tintegrals<" << getParams().getParameter(MpqcParser_Parameters::integrationParam) << ">:()" << endl;
369 *file << "\t\t)" << endl;
370 *file << "\t)" << endl;
371 } else {
372 ELOG(0, "Unknown level of theory requested for MPQC output file.");
373 }
374 const std::string jobtype = getParams().getParameter(MpqcParser_Parameters::jobtypeParam);
375 if (jobtype == getParams().getJobtypeName(MpqcParser_Parameters::Optimization)) {
376 *file << "\t% optimizer object for the molecular geometry" << endl;
377 *file << "\topt<QNewtonOpt>: (" << endl;
378 *file << "\t\tfunction = $..:mole" << endl;
379 *file << "\t\tupdate<BFGSUpdate>: ()" << endl;
380 *file << "\t\tconvergence<MolEnergyConvergence>: (" << endl;
381 *file << "\t\t\tcartesian = yes" << endl;
382 *file << "\t\t\tenergy = $..:..:mole" << endl;
383 *file << "\t\t)" << endl;
384 *file << "\t)" << endl;
385 }
386 *file << ")" << endl;
387 *file << "molecule<Molecule>: (" << endl;
388 *file << "\tunit = " << (World::getInstance().getConfig()->GetIsAngstroem() ? "angstrom" : "bohr" ) << endl;
389 *file << "\t{ atoms geometry } = {" << endl;
390 // output of atoms
391 BOOST_FOREACH(const atom *_atom, atoms) {
392 OutputMPQCLine(file, *_atom, &center);
393 }
394 *file << "\t}" << endl;
395 *file << ")" << endl;
396 *file << "basis<GaussianBasisSet>: (" << endl;
397 *file << "\tname = \"" << getParams().getParameter(MpqcParser_Parameters::basisParam) << "\"" << endl;
398 *file << "\tmolecule = $:molecule" << endl;
399 *file << ")" << endl;
400 if (theory == getParams().getTheoryName(MpqcParser_Parameters::MBPT2_R12)) {
401 *file << "% auxiliary basis set specification" << endl;
402 *file << "\tabasis<GaussianBasisSet>: (" << endl;
403 *file << "\tname = \"" << getParams().getParameter(MpqcParser_Parameters::aux_basisParam) << "\"" << endl;
404 *file << "\tmolecule = $:molecule" << endl;
405 *file << ")" << endl;
406 }
407 }
408}
409
410
Note: See TracBrowser for help on using the repository browser.