/*************************************************************************************
* *
* ***** 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 : atom.hpp *
* Author : Tokarev. K *
* Last modification : 2008/08/05 *
* Description : This module contains functions used to convert data from XML *
* nodes to program easy-to-use form *
* *
*************************************************************************************/
#include "input.hpp"
#include <vector>
using namespace std;
extern const double Cels, atm, cal, Hartree, HartreeJ, pi, Na;
string InpEnergy, InpLength, InpP, InpT;
void ReadUnits (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);
void Move2pi (double &x);
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";
// Section <units> -> <inp>
if (Input.has_units_inp())
ReadUnits(Input.firstchild("input")->firstchild("units")->firstchild("inp"));
// 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);
}
/*/ Test symmetry
fout << "Recommended sigma = " << Mol->GetSigma() << endl;
fout << endl;
//Mol->SetSigma(sigma);*/
/*// Update molecule's properties
if(!Mol->ReCompute())
return 0; // Exit on errors (unknown atoms occured)
Mol->LogAtoms(fout);
Mol->LogInfo(fout);*/
// Section [options]
//fin >> buffer;
if ((Mol->GetNAtoms() > 2) && (Input.has_rotors()))
if(!ReadRot(Mol, Input.firstchild("input")->firstchild("rotors")))
return 0;
/* // Read "sigma = X"
fin >> buffer >> buffer >> N;
if((!N) || (N>24))
{
// fout << "## Invalid rotation factor! sigma = " << N << "\n Exiting" << endl;
#ifndef SILENT
cerr << "## Invalid rotation factor! sigma = " << N << "\n Exiting" << endl;
return 0;
#else
return <number_of_error>
#endif
}
if (N != Mol->GetSigma())
cerr << "# WARNING! Sigma value " << N << " doesn't match recommended value " << Mol->GetSigma() << endl;
Mol->SetSigma(N);*/
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;
}
bool ReadGeometry (Molecule * Mol, XMLNode * node)
{
int NAtoms;
char buf[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 >> buf;
if(iss.eof())
{
cerr << "\n## File is broken or incomplete in group <geometry>, line " << k+2 <<"\n Expected: X" << endl;
return false;
}
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;
try
{
(Mol->GetAtom(k))->SetAtom(x, y, z, buf);
}
catch (GenericError e)
{
cerr << e.get_message() << endl;
return false;
}
}
// Read point group if present
buffer = "";
buffer = node->get_attr("sym");
if (buffer!="")
//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);
if (*end != 0)
{
cerr << "Error in group <elevels>: illegal value of energy="
<< buffer.c_str() << " for level#"<< i+1 << endl;
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;
// ofstream fout;
// fout.open(LogName, ios::app);
// First of all, 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;
XMLNode * pes;
for (it=nl.begin(); j<NRot; j++, it++)
{
//Rotor * &R = Mol->GetRotor(j);
vector<int> ChopGroup;
// Read rotation type
buffer = "";
buffer = (*it)->get_attr("barrier");
if (buffer != "")
{
Mol->GetRotor(j)->SetRotType(HINDERED_HARMONIC);
height = strtod(buffer.c_str(), NULL);
if (InpEnergy == "H")
height *= HartreeJ;
else if (InpEnergy == "kJ")
height *= 1e-3/Na;
else if (InpEnergy == "kcal")
height *= 1e-3*cal/Na;
else
{
cerr << "Error in group <units> : invalid unit " << InpEnergy.c_str() << " for energy!\n";
cerr << "Valid units are: kJ kcal H \n";
return false;
}
//(R->GetBarrier()).SetHeight(height);
(Mol->GetRotor(j)->GetBarrier()).SetPES(height, Mol->GetRotor(j)->GetSigma());
}
else
{
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!
{
Mol->GetRotor(j)->SetRotType(HINDERED_DISCREET);
if(!ReadBarrier(Mol->GetRotor(j)->GetBarrier(),pes))
return false;
}
}
// Read rotating group
//fin >> NAtGroup >> buffer; // <Nuber of atoms> :
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 = (*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);
}
/*/ Special data for different rotation types
int N;
double barrier, x, y, y0;
double * RotParams = R->GetRotParams();
switch (R->GetRotType())
{
case FREE:
RotParams = NULL;
break;
case HINDERED_HARMONIC:
RotParams = new double[2];
// Read "barrier = "
fin >> buffer >> buffer >> barrier;
RotParams[0] = barrier*HartreeJ;
RotParams[1] = R->GetSigma();
fout << "Rotation Parameters: " << RotParams[0] << "\t" << RotParams[1] << endl;
break;
default:
#ifndef SILENT
cerr << "Oops! Something went wrong!\n";
fout.close();
return false;
#else
return <number_of_error>
#endif
}
R->SetRotParams(RotParams); */
}
//fout.close();
return true;
}
bool ReadBarrier (PES_1D & Bar, XMLNode * node)
{
// Read type type of PES
string buffer = "";
buffer = node->get_attr("type");
if ((buffer=="")||(buffer=="discreet"))
{
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 >> 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
{
Bar.SetPES(Energy, FOURIER, 3);
}
catch (GenericError e)
{
cerr << e.get_message() << endl;
return false;
}
/*case HINDERED_DISCREET:
// Read "xy <N>"
fin >> buffer >> N;
RotParams = new double[2*N+1];
RotParams[0] = N;
// xy array
y0 = 0;
for (int m=0; m<N; m++)
{
fin >> x >> y;
if (y < y0)
y0 = y;
RotParams[2*m+1] = x*pi/180;
RotParams[2*m+2] = y*HartreeJ;
}
fout << "Energy minimum: ";
fout << y0 << " Hartree" << endl;
for (int m=0; m<N; m++)
{
RotParams[2*m+2] -= y0*HartreeJ;
fout << RotParams[2*m+1] << "\t" << RotParams[2*m+2] << endl;
}
break;*/
return true;
}
if (buffer=="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 (buffer=="cos")
{
return true;
}
return false;
}
bool ReadFreq(Molecule * Mol, XMLNodeList nl)
{
XMLNodeListIterator it;
vector<double> Freq;
int j;
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))
{
cout << "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))
{
cout << "## 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;
}
}
return true;
}
bool ReadFreqList (vector<double> & Freq, istream & iss, string style)
{
if ((style=="") || (style=="std"))
{
int i=1;
do
{
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++;
}while (!iss.eof());
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;
}
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;
}
iss >> f;
if ((n <= 0))
{
cerr << "Error in group <freq>: illegal number value " << n << " at line " << i << endl;
return false;
}
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;
}
}
/*
bool ReadFreq(Molecule * Mol, int & Range, XMLInputFile & Input, ofstream & fout)
{
// Reads description of rotations in the input file
// Returns FALSE on fatal error
char freqstyle[10], buffer[80];
double freqscale = 0;
double f = 0, f2 = 0;
int NFreq = 0;
// Read "freqstyle = X"
fin >> buffer >> buffer >> freqstyle;
if(strcmp(freqstyle, "molpro")*strcmp(freqstyle, "std")*strcmp(freqstyle, "anharmonic"))
{
cout << "## Invalid frequencies style! freqstyle = " << freqstyle << "\n Exiting" << endl;
return false;
}
// Read "range = X"
if (!strcmp(freqstyle, "anharmonic"))
{
fin >> buffer >> buffer >> Range;
if(Range <= 0)
{
cout << "## Invalid anharmonic range! range = " << Range << "\n Exiting" << endl;
return false;
}
}
// Read "freqscale = X"
fin >> buffer >> buffer >> freqscale;
if((freqscale <= 0) || (freqscale > 10))
{
cout << "## Invalid frequencies scale! freqscale = " << freqscale << "\n Exiting" << endl;
return false;
}
if (freqscale != 1.0)
fout << "Frequencies will be used with scale factor " << freqscale << endl;
fout << "Using frequencies list:" << endl;
// Read frequencies
fin >> buffer;
int NAtoms = Mol->GetNAtoms();
//Freq = new double[3*NAtoms];
while (fin)
{
f = 0;
if (!strcmp(freqstyle, "molpro"))
{
fin >> buffer;
if (strcmp(buffer, "[comments]")) // Stop if [comments] start
fin >> f; // In the input: "Number Frequency"
else
f = 0;
if (buffer[0] == '#')
continue; // Skip frequency
Mol->SetFreq(NFreq, f*freqscale);
}
if (!strcmp(freqstyle, "anharmonic"))
{
fin >> buffer;
if (strcmp(buffer, "[comments]")) // Stop if [comments] start
fin >> f >> f2; // In the input: "Number Harmonic_Frequency Anharmonic_Frequency"
else
{
f = 0;
f2 = 0;
}
if (buffer[0] == '#')
continue; // Skip frequency
Mol->SetFreq(NFreq, f*freqscale);
Mol->SetFreq_Anh(NFreq, f2*freqscale); /// ??? is it sensible to scale anharmonics?
}
if(!f)
continue; // Skip zero frequencies
fout << NFreq+1 << "\t" << Mol->GetFreq(NFreq) << "\t\t" << Mol->GetFreq_Anh(NFreq) << endl;
NFreq++;
}
if (!strcmp(freqstyle, "anharmonic"))
Mol->SetNFreq(NFreq, true);
else
Mol->SetNFreq(NFreq, false);
return true;
}*/
#ifndef SILENT
bool GetFileName (string & FileName)
{
// Returns TRUE if successfully opened
// Returns FALSE on exit
ifstream file;
cout << "Input file name:\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 << "\n File reading error... Input path again:\n>";
file.close();
}
file.close();
return true;
}
#endif
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 << "Error 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 << "Error 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 (InpT == "C")
{
TStart += Cels;
TStop += Cels;
}
else if (InpT != "K")
{
cerr << "Error in group <units>: invalid unit " << InpEnergy.c_str() << " for temperature!\n";
cerr << "Valid units are: K C \n";
return false;
}
if (dT==0)
{
cerr << "Error in group <job>: zero temperature step!\n";
return false;
}
if (dT>0)
if (TStop < TStart)
{
cerr << "Error in group <job>: TStart must be greater then TStop!\n";
return false;
}
if (dT<0)
if (TStop > TStart)
{
cerr << "Error in group <job>: TStart must be lesser then TStop!\n";
return false;
}
return true;
}
}
void Move2pi (double &x)
{
while (true)
{
if (x>=360)
x-=360;
else
break;
}
}