/*************************************************************************************
* *
* ***** OpenThermo ***** *
* Calculation of thermodynamic functions from molecular data *
* Copyright 2008 Konstantin Tokarev <annulen@users.sourceforge.net> *
* and others *
* *
*************************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License. See COPYING for *
* more details *
* *
*************************************************************************************
* Module name : input.cpp *
* Author : Tokarev. K *
* Last modification : 2009/04/05 *
* Description : This module contains functions used to convert data from XML *
* nodes to program easy-to-use form *
* *
*************************************************************************************/
#include "control.hpp"
#include "input.hpp"
#include "thermo.h"
#include <vector>
using namespace std;
// Globals
extern const double Cels, atm, cal, Hartree, HartreeJ, pi, Na;
namespace
{
string InpEnergy, InpLength, InpP, InpT;
Rotation_Theory DefaultRotTheory;
}
// File reading stages
void ReadUnits (XMLNode * node);
void ReadJob (XMLNode * node);
bool ReadGeometry (Molecule * Mol, XMLNode * node);
bool ReadLevels (Molecule * Mol, XMLNode * node);
bool ReadRot(Molecule * Mol, XMLNode * node);
bool ReadBarrier (PES_1D & Bar, XMLNode * node);
bool ReadFreq(Molecule * Mol, XMLNodeList nl);
bool ReadFreqList (vector<double> & Freq, istream & iss, string style);
// Utils
inline void ReadString (string in, string & out);
inline void ReadBool (string in, bool & out);
inline void ReadTSet (string in, TemperatureSet & out);
bool ConvertEnergy (double & e);
void Move2pi (double &x);
// Definitions
bool ReadMolData(Molecule * Mol, XMLInputFile & Input)
{
// char buf[80];
// string buffer("");
// ofstream fout;
// fout.open(LogName, ios::app);
// Standard units for input values
::InpEnergy = "kj";
::InpLength = "a";
::InpP = "pa";
::InpT = "k";
::DefaultRotTheory = CLASSIC;
// Section <units> -> <inp>
if (Input.has_units_inp())
ReadUnits(Input.firstchild("input")->firstchild("units")->firstchild("inp"));
// Section <job>
if (Input.has_job())
ReadJob(Input.firstchild("input")->firstchild("job"));
// Section <geometry>
if (!ReadGeometry(Mol,Input.firstchild("input")->firstchild("geometry")))
return false;
// Section [electronic_levels]
if (Input.has_electronic_levels())
{
if (!ReadLevels(Mol,Input.firstchild("input")->firstchild("elevels")))
return false;
}
else
{
Mol->SetNLevels(1);
Mol->SetLevel(0, 1, 0);
}
if ((Mol->GetNAtoms() > 2) && (Input.has_rotors()))
if(!ReadRot(Mol, Input.firstchild("input")->firstchild("rotors")))
return 0;
if(!ReadFreq(Mol, Input.firstchild("input")->children("freq")))
return false;
return true;
}
void ReadUnits (XMLNode * node)
{
/* string buffer = "";
buffer = node->get_attr("energy");
if (buffer != "")
::InpEnergy = buffer;
buffer = "";
buffer = node->get_attr("length");
if (buffer != "")
::InpLength = buffer;
buffer = "";
buffer = node->get_attr("p");
if (buffer != "")
::InpP = buffer;
buffer = "";
buffer = node->get_attr("t");
if (buffer != "")
::InpT = buffer;*/
ReadString(node->get_attr("energy"),::InpEnergy);
ReadString(node->get_attr("length"),::InpLength);
ReadString(node->get_attr("p"),::InpP);
ReadString(node->get_attr("t"),::InpT);
}
void ReadJob (XMLNode * node)
{
try
{
ReadTSet(node->get_attr("tset"), Control::itsT_Set);
}
catch (GenericError e)
{
throw GenericError("Error in tset parameter in <job> group: " + string(e.get_message()));
return;
}
try
{
ReadBool((node)->get_attr("spin"),Control::itsSpin);
}
catch (GenericError e)
{
throw GenericError("Error in spin parameter in <job> group: " + string(e.get_message()));
return;
}
}
bool ReadGeometry (Molecule * Mol, XMLNode * node)
{
int NAtoms;
char name[80];
string buffer = GetCdata(node);
#ifdef STRSTREAM
istrstream iss(buffer.c_str(), buffer.length())
#else /* STRSTREAM */
istringstream iss(buffer,istringstream::in);
#endif /* STRSTREAM */
iss >> NAtoms;
if (NAtoms <=0)
{
cerr << "Error in group <geometry>: illegal number of atoms " << NAtoms << endl;
return false;
}
Mol->SetAtoms(NAtoms);
for (int k=0; k<NAtoms; k++)
{
double x, y, z;
if(iss.eof())
{
cerr << "\n## File is broken or incomplete in group <geometry>, line " << k+2 <<"\n Expected: atomic symbol" << endl;
return false;
}
iss >> name;
if(iss.eof())
{
cerr << "\n## File is broken or incomplete in group <geometry>, line " << k+2 <<"\n Expected: X" << endl;
return false;
}
name[0]=toupper(name[0]); // Element's name starts with capital letter
iss >> x;
if(iss.eof())
{
cerr << "\n## File is broken or incomplete in group <geometry>, line " << k+2 <<"\n Expected: Y" << endl;
return false;
}
iss >> y;
if(iss.eof())
{
cerr << "\n## File is broken or incomplete in group <geometry>, line " << k+2 <<"\n Expected: Z" << endl;
return false;
}
iss >> z;
if (::InpLength == "nm")
{
x *= 10;
y *= 10;
z *= 10;
}
else if (::InpLength == "cm")
{
x *= 1e8;
y *= 1e8;
z *= 1e8;
}
else if (::InpLength == "m")
{
x *= 1e10;
y *= 1e10;
z *= 1e10;
}
else if (::InpLength !="a")
{
cerr << "Error in group <units> : invalid unit " << ::InpLength.c_str() << " for length!\n";
cerr << "Valid units are: A nm cm m\n";
return false;
}
try
{
(Mol->GetAtom(k))->SetAtom(x, y, z, name);
}
catch (GenericError e)
{
cerr << e.get_message() << endl;
return false;
}
}
// Read point group if present
buffer = "";
buffer = node->get_attr("sym");
if (buffer!="")
{
buffer.at(0) = toupper(buffer.at(0));
//try
Mol->SetPointGroup(buffer.c_str());
// catch
}
return true;
}
bool ReadLevels (Molecule * Mol, XMLNode * node)
{
int N;
char * end = 0;
N=CountNodes(node,"level");
Mol->SetNLevels(N);
// Read levels
int d;
double E;
XMLNodeListIterator it;
XMLNodeList nl;
nl = node->children("level");
int i=0;
for (it=nl.begin(); i<N; i++, it++)
{
string buffer = "";
buffer = (*it)->get_attr("degen");
if (buffer == "")
d = 1;
else
{
d = strtol(buffer.c_str(),&end,0);
if ((*end != 0)||(d<=0))
{
cerr << "Error in group <elevels>: illegal value of degen="
<< buffer.c_str() << " for level#"<< i+1 << endl;
return false;
}
}
buffer = "";
buffer = (*it)->get_attr("energy");
if (buffer == "")
{
cerr << "Error in group <elevels>: energy value for level#"
<< i+1 << " not found!\n";
return false;
}
end = 0;
E = strtod(buffer.c_str(),&end);
// read unit if given
if (*end != 0)
{
while (end[0] == ' ')
end++;
buffer = ::InpEnergy;
::InpEnergy = end;
if(!ConvertEnergy(E))
return false;
::InpEnergy = buffer;
}
else
{
if(!ConvertEnergy(E))
return false;
}
Mol->SetLevel(i, d, E);
}
return true;
}
bool ReadRot(Molecule * Mol, XMLNode * node)
{
// Reads description of rotations in the node <rotors> and calculates Ieff
// Returns FALSE on error in input
string buffer("");
int NRot, sigma, k;
double height;
XMLNodeList nl;
XMLNodeListIterator it;
// First of all, load global settings for rotors
// Separate or S Determinant
try
{
ReadBool((node)->get_attr("sdet"),Control::itsS_Det);
}
catch (GenericError e)
{
throw GenericError("Error in sdet parameter in <rotors> group: " + string(e.get_message()));
return false;
}
// Global theory
buffer = "";
buffer = (node)->get_attr("theory");
if (buffer != "")
{
if (Control::S_Det() == false)
{
if (buffer=="classic")
::DefaultRotTheory = CLASSIC;
else if (buffer=="quantum")
::DefaultRotTheory = QUANTUM;
else
throw GenericError("Error in <rotors> group: unknown theory "
+ buffer + "\nValid keywords are: classic quantum");
}
else
{
if (buffer == "quantum")
throw GenericError("S Determinant calculations are currently implemented only for classic theory");
//::DefaultRotTheory = CLASSIC;
}
}
// count rotors
NRot = CountNodes(node, "rotor");
if (NRot == 0)
{
cerr << "## Error! Zero rotors number!" << endl; // Should not be reached
return false;
}
Mol->SetNRot(NRot);
nl = node->children("rotor");
int j=0;
for (it=nl.begin(); j<NRot; j++, it++)
{
//Rotor * &R = Mol->GetRotor(j);
vector<int> ChopGroup;
buffer = "";
buffer = (*it)->get_attr("theory");
if (buffer=="")
Mol->GetRotor(j)->SetTheory(::DefaultRotTheory);
else if (buffer=="classic")
Mol->GetRotor(j)->SetTheory(CLASSIC);
else if (buffer=="quantum")
Mol->GetRotor(j)->SetTheory(QUANTUM);
else
throw GenericError("Error in <pes> : Unknown theory="+buffer);
if (Control::S_Det())
{
if (buffer == "quantum")
throw GenericError("S Determinant calculations are currently implemented only for classic theory");
Mol->GetRotor(j)->SetTheory(CLASSIC);
}
// Determine rotation type and read rotation barrier if present
buffer = "";
buffer = (*it)->get_attr("barrier");
if (buffer != "")
{
char * end = 0;
if (Control::S_Det())
throw GenericError("S Determinant calculations are currently implemented only for free rotation");
Mol->GetRotor(j)->SetRotType(HINDERED_HARMONIC);
height = strtod(buffer.c_str(), &end);
// read unit if given
if (*end != 0)
{
while (end[0] == ' ')
end++;
buffer = ::InpEnergy;
::InpEnergy = end;
if(!ConvertEnergy(height))
return false;
::InpEnergy = buffer;
}
else
{
if (!ConvertEnergy(height))
return false;
}
//(R->GetBarrier()).SetHeight(height);
(Mol->GetRotor(j)->GetBarrier()).SetPES(height, Mol->GetRotor(j)->GetSigma());
}
else
{
XMLNode * pes = 0;
try
{
pes = (*it)->firstchild("pes");
}
catch (xmlerror e)
{
if (e.get_error() == xml_name_not_found)
Mol->GetRotor(j)->SetRotType(FREE);
else
return false;
}
if(Mol->GetRotor(j)->GetRotType() != FREE) // Free isn't default!
{
if (Control::S_Det())
throw GenericError("S Determinant calculations are currently implemented only for free rotation");
Mol->GetRotor(j)->SetRotType(HINDERED_DISCREET);
try
{
if(!ReadBarrier(Mol->GetRotor(j)->GetBarrier(),pes))
return false;
}
catch (GenericError e)
{
cerr << "Error in <rotor> #" << j+1 << ":\n\t" << e.get_message() << endl;
return false;
}
}
}
// Read rotating group
buffer = "";
buffer = (*it)->get_attr("atoms");
#ifdef STRSTREAM
istrstream iss1(buffer.c_str(), buffer.length())
#else /* STRSTREAM */
istringstream iss1(buffer,istringstream::in);
#endif /* STRSTREAM */
ChopGroup.reserve(Mol->GetNAtoms());
while (!iss1.eof())
{
iss1 >> k;
ChopGroup.push_back(k);
if(k<=0)
{
cerr << "Error in <rotor> #" << j+1 << ": invalid number of atom " << k << endl;
return false;
}
}
// Read "bond : N1 N2"
buffer = "";
buffer = (*it)->get_attr("bond");
#ifdef STRSTREAM
istrstream iss2(buffer.c_str(), buffer.length())
#else /* STRSTREAM */
istringstream iss2(buffer,istringstream::in);
#endif /* STRSTREAM */
int bond[2];
if(iss2.eof())
{
cerr << "Error in <rotor> #" << j+1 <<": bonding atoms expected" << endl;
return false;
}
iss2 >> bond[0];
if(iss2.eof())
{
cerr << "Error in <rotor> #" << j+1 <<": second bonding atom expected" << endl;
return false;
}
iss2 >> bond[1];
if ((bond[0]>Mol->GetNAtoms()) || (bond[0]<1))
{
cerr << "Error in <rotor> #" << j+1 <<": incorrect value " << bond[0] << " for bonding atom\n";
return false;
}
if ((bond[1]<1) || (bond[1]<1))
{
cerr << "Error in <rotor> #" << j+1 <<": incorrect value " << bond[1] << " for bonding atom\n";
return false;
}
try
{
Mol->SetRotor(j, ChopGroup, bond);
}
catch (GenericError e)
{
cerr << e.get_message() << endl;
return false;
}
// Read sigma if present
buffer = "";
buffer = (*it)->get_attr("sigma");
if (buffer != "")
{
char * end = 0;
sigma = strtol(buffer.c_str(),&end,0);
if ((*end != 0)||(sigma<=0)||(sigma>24))
{
cerr << "Error in <rotor> #" << j+1 << ": illegal value of rotation number sigma="
<< buffer.c_str() << endl;
}
/*
if ((sigma != R->GetSigma()) && R->GetSigma())
cerr << "# WARNING! Sigma value " << sigma << " for rotor #" << j+1 << " doesn't match recommended value " << R->GetSigma() << endl;*/
Mol->GetRotor(j)->SetSigma(sigma);
}
}
//fout.close();
return true;
}
bool ReadBarrier (PES_1D & Bar, XMLNode * node)
{
// Read type type of PES
string type = "", buffer = "", degree = "";
int MaxDeg = 0;
type = node->get_attr("type");
degree = node->get_attr("degree");
if ((type=="")||(type=="discreet"))
{
if (degree != "")
{
char * end = 0;
MaxDeg = strtol(degree.c_str(),&end,0);
if ((*end != 0)||(MaxDeg<=0))
{
throw GenericError("Error in <pes> : illegal value of degree="+ degree);
//cerr << "Error in <pes> : illegal value of maximal interpolation degree = "
// << degree.c_str() << endl;
}
}
map<double,double> Energy;
buffer = GetCdata(node);
#ifdef STRSTREAM
istrstream iss(buffer.c_str(), buffer.length())
#else /* STRSTREAM */
istringstream iss(buffer,istringstream::in);
#endif /* STRSTREAM */
do
{
double x, y;
iss >> x;
if (iss.eof())
break;
iss >> y;
Move2pi(x);
x = x*pi/180;
// Units conversion
// Hartree -> J/mol
Energy[x] = y*Hartree;
}while(!iss.eof());
// init: energy values + approx. type
try
{
if (MaxDeg != 0)
Bar.SetPES(Energy, FOURIER, MaxDeg);
else
Bar.SetPES(Energy, FOURIER);
}
catch (GenericError e)
{
cerr << e.get_message() << endl;
return false;
}
return true;
}
if (type=="polynom")
{
/* case HINDERED_POLYNOM:
// Read "polynom <N> :"
fin >> buffer >> N >> buffer;
RotParams = new double[N+2];
RotParams[0] = N;
fout << "Rotation Parameters: " << RotParams[0] << "\t";
// minE = y0
// Coefficients
for (int n=0; n<N; n++)
{
fin >> y;
RotParams[n+1] = y*HartreeJ;
fout << RotParams[n+1] << "\t";
}
fin >> y;
RotParams[N+1] = y*HartreeJ;
fin >> buffer >> buffer >> x ;
y0 = 0;
for (int n=0; n<=N; n++)
y0 += RotParams[n+1]*pow(x, N-n);
fout << buffer << "Energy minimum: " << y0 << " J" << endl;
RotParams[N+1] -= y0;
fout << RotParams[N+1] << endl;
break;*/
return true;
}
if (type=="cos")
{
return true;
}
return false;
}
bool ReadFreq(Molecule * Mol, XMLNodeList nl)
{
XMLNodeListIterator it;
vector<double> Freq;
int j=0;
for (it=nl.begin(), j=1; it!=nl.end(); it++, j++)
{
string buffer = GetCdata(*it);
#ifdef STRSTREAM
istrstream iss(buffer.c_str(), buffer.length())
#else /* STRSTREAM */
istringstream iss(buffer,istringstream::in);
#endif /* STRSTREAM */
string type = "", style = "";
type = (*it)->get_attr("type");
style = (*it)->get_attr("style");
if (style=="")
{
//if ((type=="") || (type=="harm"))
style = "std";
//else
// style = "selected";
}
if (!ReadFreqList(Freq, iss, style))
{
cerr << "Group <freq> #" << j << " incorrect!\n";
return false;
}
if ((type=="") || (type=="harm"))
{
string scale = "";
double f = -1;
scale = (*it)->get_attr("scale");
if (scale != "")
{
f = strtod(scale.c_str(), NULL);
if((f <= 0) || (f > 10))
{
cerr << "## Invalid frequencies scale factor! freqscale = " << f << "\n Exiting" << endl;
return false;
}
for (unsigned int i=0; i<Freq.size(); i++)
Freq.at(i) *= f;
//ScaleFreq(Freq, f);
}
Mol->SetFreq(Freq);
}
else if (type=="anharm")
{
try
{
Mol->SetAnh(Freq);
}
catch (GenericError e)
{
cerr << e.get_message() << endl;
return false;
}
}
else if (type=="exp")
{
try
{
Mol->SetFreqExp(Freq);
}
catch (GenericError e)
{
cerr << e.get_message() << endl;
return false;
}
}
else
{
cerr << "_oops!\n";
return false;
}
Freq.clear();
}
return true;
}
bool ReadFreqList (vector<double> & Freq, istream & iss, string style)
{
if ((style=="") || (style=="std"))
{
int i=0;
while (!iss.eof())
{
double f=0;
iss >> f;
if (f < 0)
{
cerr << "Error in group <freq>: illegal frequency value " << f << " at line " << i << endl;
return false;
}
if (f == 0)
continue;
Freq.push_back(f);
i++;
}
return true;
}
else if (style=="molpro")
{
int i=1;
do
{
int n=0;
double f=0;
iss >> n;
if(iss.eof())
{
cerr << "Error in group <freq>: broken line " << i << endl;
return false;
}
if (n == 0)
continue;
iss >> f;
if (f <= 0)
{
cerr << "Error in group <freq>: illegal frequency value " << f << " at line " << i << endl;
return false;
}
Freq.push_back(f);
i++;
}while (!iss.eof());
return true;
}
else if (style=="selected") // ??????
{
map<int,double> FreqMap;
int i=1, maxn=0;
do
{
double f=0;
int n=0;
iss >> n;
if(iss.eof())
{
cerr << "Error in group <freq>: broken line " << i << endl;
return false;
}
if ((n < 0))
{
cerr << "Error in group <freq>: illegal number value " << n << " at line " << i << endl;
return false;
}
if (n == 0)
continue;
iss >> f;
if (f <= 0)
{
cerr << "Error in group <freq>: illegal frequency value " << f << " at line " << i << endl;
return false;
}
FreqMap[n] = f;
if (n>maxn)
maxn = n;
i++;
}while (!iss.eof());
// ????????????
/*for (int j=0; j<maxn; j++)
Freq.push_back(FreqMap[j]);*/
return true;
}
else
{
cerr << "Error in group <freq> : invalid frequencies style \"" << style.c_str() << "\"\n";
return false;
}
}
#ifndef SILENT
bool GetFileName (string & FileName)
{
// Returns TRUE if successfully opened
// Returns FALSE on exit
ifstream file;
cout << "Input file name with path:\tInput '#' to exit\n>";
while(true)
{
const char exit = '#';
cin >> FileName;
cin.ignore(1, '\n');
if (FileName[0] == exit) return false;
file.open(FileName.c_str()); // Temporary object
if (file.is_open())
break;
cout << "Can't open" << FileName.c_str() << endl;
FileName += ".inp";
file.open(FileName.c_str()); // Temporary object
if (file.is_open())
break;
cout << "Can't open " << FileName.c_str() << endl;
cout << "\n Input file name with path again:\n>";
file.close();
}
file.close();
return true;
}
#endif
void Preprocess (string FileName, string & PFileName)
{
#ifdef POSIX
PFileName = "/tmp/openthermo";
#else
PFileName = "tmp";
#endif
ifstream fin(FileName.c_str());
ofstream fout(PFileName.c_str());
char Line[250];
char lowerLine[250];
while (!fin.eof())
{
fin.getline(Line, 250);
//lowerLine = Line;
unsigned int i;
for (i=0; (i<strlen(Line))&&(i<250); i++)
lowerLine[i] = tolower(Line[i]);
lowerLine[i] = '\0';
//lowerLine.tolower();
if ((Line[0] != '#') && (Line[0] != '%'))
fout << lowerLine << "\n";
else
fout << "\n";
}
fin.close();
fout.close();
}
bool ReadP(XMLInputFile& i, double& p)
{
if (!i.is_valid() || !i.has_job())
return false;
if (i.has_units_inp())
ReadUnits(i.firstchild("input")->firstchild("units")->firstchild("inp"));
string s("");
s = i.firstchild("input")->firstchild("job")->get_attr("p");
if (s=="")
return false;
else
{
#ifdef STRSTREAM
istrstream iss(s.c_str(), s.length())
#else /* STRSTREAM */
istringstream iss(s,istringstream::in);
#endif /* STRSTREAM */
iss >> p;
if (::InpP == "atm")
p *= atm;
else if (::InpP != "pa")
{
cerr << "\nError in group <units>: invalid unit " << InpEnergy.c_str() << " for pressure!\n";
cerr << "Valid units are: Pa atm \n";
return false;
}
if (p>0)
return true;
else
{
cerr << "\nError in group <job>: pressure must be positive!\n";
return false;
}
}
}
bool ReadT(XMLInputFile& i, double& TStart, double& dT, double& TStop)
{
if (!i.is_valid() || !i.has_job())
return false;
if (i.has_units_inp())
ReadUnits(i.firstchild("input")->firstchild("units")->firstchild("inp"));
string s("");
s = i.firstchild("input")->firstchild("job")->get_attr("t");
if (s=="")
return false;
else
{
#ifdef STRSTREAM
istrstream iss(s.c_str(), s.length())
#else /* STRSTREAM */
istringstream iss(s,istringstream::in);
#endif /* STRSTREAM */
iss >> TStart >> dT >> TStop;
if (dT>0)
if (TStop < TStart)
{
cerr << "\nError in group <job>: TStart must be greater then TStop!\n";
return false;
}
if (dT<0)
if (TStop > TStart)
{
cerr << "\nError in group <job>: TStart must be lesser then TStop!\n";
return false;
}
if (dT==0)
{
if (TStop==0) // Single T run
{
dT = 1;
TStop = TStart;
}
else
{
cerr << "\nError in group <job>: zero temperature step!\n";
return false;
}
}
if (::InpT == "c")
{
TStart += Cels;
TStop += Cels;
}
else if (::InpT != "k")
{
cerr << "\nError in group <units>: invalid unit " << InpEnergy.c_str() << " for temperature!\n";
cerr << "Valid units are: K C \n";
return false;
}
return true;
}
}
inline void ReadString (string in, string & out)
{
if (in != "")
out = in;
}
inline void ReadBool (string in, bool & out)
{
if (in != "")
out = XML_True(in);
}
inline void ReadTSet (string in, TemperatureSet & out)
{
if (in != "")
out = String2TSet(in);
}
bool ConvertEnergy (double & e)
{
if (::InpEnergy == "h")
e *= Hartree;
else if (::InpEnergy == "j")
e *= 1;
else if (::InpEnergy == "kj")
e *= 1e3;
else if (::InpEnergy == "kcal")
e *= 1e-3*cal;
else
{
cerr << "Error in group <units> : invalid unit " << ::InpEnergy.c_str() << " for energy!\n";
cerr << "Valid units are: J kJ kcal H \n";
e = -1; // Not safe! exception shoul be used
return false;
}
return true;
}
void Move2pi (double &x)
{
while (true)
{
if (x>=360)
x-=360;
else if (x<0)
x+=360;
else
break;
}
}