/*************************************************************************************
* *
* ***** 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.cpp *
* Author : Tokarev. K *
* Last modification : 2008/07/23 *
* Description : This module contains definitions of class 'Atom', *
* which is used for storage of atomic symbols, coordinates, *
* weights and other properties *
* *
*************************************************************************************/
#include <fstream>
#include <iomanip>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <string>
#include "atom.hpp"
#include "control.hpp"
#include "hash.h"
#include "matrix.h"
#include "exception.hpp"
#include "xml.hpp"
using namespace std;
bool CreateTable();
bool InitializeTable();
struct Isotope
{
int Number;
double Spin;
double Weight;
double Abundance;
};
struct Element
{
int Number;
double Weight;
Hash<int,Isotope> Isotopes;
int itsMainIsotope;
Isotope MainIsotope ()
{
return Isotopes[itsMainIsotope];
}
};
bool isInitialized = false;
string ElementsFileName, IsotopesFileName;
int NElements;
Hash <string,Element> PeriodicTable;
double Spin2Double (string s)
{
char spin[8];
strncpy(spin, s.c_str(), 8);
// Last symbol maybe + or -
int l = strlen(spin);
if (spin[l-1] == '+' || spin[l-1] == '-')
spin[l-1] = '\0';
l = strlen(spin);
// Now last symbols may be /2
if (l>2)
{
if (spin[l-2] == '/')
{
spin[l-2] = '\0';
return strtol(spin, 0, 10) / 2.0;
}
else
return strtol(spin, 0, 10);
}
else
return strtol(spin, 0, 10);
}
// Definitions of Atom
Atom::Atom () : itsAtomicNumber(0), itsIsotopeNumber(0), itsNuclearSpin(0),
itsX(0), itsY(0), itsZ(0), itsWeight(0)
{
#ifdef OBJECT_COUNTERS
NumOfAtoms++;
#endif
strcpy(itsName,"Xx");
}
Atom::~Atom ()
{
#ifdef OBJECT_COUNTERS
NumOfAtoms--;
#endif
}
void Atom::SetAtom (double x, double y, double z, char * Name, int i)
{
itsX = x;
itsY = y;
itsZ = z;
#ifndef SILENT
if (strlen(Name) > ATOM_NAME_LENGTH)
cerr << "# WARNING : Strange atom name encountered! Truncating to "<< ATOM_NAME_LENGTH << " symbols..." << endl;
#endif
strncpy (itsName,Name, ATOM_NAME_LENGTH);
if(!InitializeTable())
{
throw GenericError("## ERROR : Periodic Table not found or corrupted!");
return;
}
itsAtomicNumber = PeriodicTable[Name].Number;
if (i == 0)
{
itsIsotopeNumber = 0;
// Natural mixture
itsWeight = PeriodicTable[Name].Weight;
// Spin of the most abundant isotope
itsNuclearSpin = PeriodicTable[Name].MainIsotope().Spin;
}
else if (i == -1)
{
itsIsotopeNumber = PeriodicTable[Name].itsMainIsotope;
// Weight of the most abundant isotope
itsWeight = PeriodicTable[Name].MainIsotope().Weight;
// Spin of the most abundant isotope
itsNuclearSpin = PeriodicTable[Name].MainIsotope().Spin;
}
else
{
if (PeriodicTable[Name].Isotopes.find(i) != PeriodicTable[Name].Isotopes.end())
{
itsIsotopeNumber = i;
itsWeight = PeriodicTable[Name].Isotopes[i].Weight;
itsNuclearSpin = PeriodicTable[Name].Isotopes[i].Spin;
}
else
{
cerr << "Illegal isotope number " << i << endl;
throw GenericError("## ERROR : Illegal isotope number!");
return;
}
}
//itsWeight = PeriodicTable[Name].Weight;
/*if (!itsWeight)
{
string e = "## ERROR : Element \"";
e += Name;
e += "\" wasn't fount in ";
e += ElementsFileName;
e += " or has zero mass";
throw GenericError(e);
return;
}*/
//itsAtomicNumber = PeriodicTable[Name].Number;
//itsNuclearSpin = PeriodicTable[Name].Spin;
//itsNuclearSpin = PeriodicTable[Name].MainIsotope().Spin;
}
void Atom::SetCoordinates (double x, double y, double z)
{
itsX = x;
itsY = y;
itsZ = z;
}
void Atom::GetCoordinates (double & x, double & y, double & z) const
{
x = itsX;
y = itsY;
z = itsZ;
}
void Atom::GetCoordinates (ColumnVector &R) const
{
R.ReSize(3);
R << itsX << itsY << itsZ;
}
/*double Atom::GetWeight () const
{
return itsWeight;
}
double Atom::GetWeightSI () const
{
return itsWeight*(1e-3/Na);
}
int Atom::GetAtomicNumber () const
{
return itsAtomicNumber;
}
const char * Atom::GetName () const
{
return itsName;
}*/
double Atom::GetNuclearSpin () const
{
return itsNuclearSpin;
}
ostream & operator<< (ostream & theStream, Atom & theAtom)
{
/*#ifndef __WATCOMC__
theStream << setw(ATOM_NAME_LENGTH) << left << theAtom.itsName << "\t";
theStream << setprecision(8);
theStream << setw(12) << showpoint << right << theAtom.itsX << "\t";
theStream << setw(12) << theAtom.itsY << "\t" ;
theStream << setw(12) << theAtom.itsZ << "\t";
theStream << setprecision(9);
theStream << setw(13) << resetiosflags(ios::showpoint) << left << theAtom.itsWeight << "\t";
theStream << setw(3) << setprecision(1) << theAtom.itsNuclearSpin << endl;
#else*/
theStream << setw(ATOM_NAME_LENGTH) << setiosflags(ios::left) << theAtom.itsName << "\t";
theStream << setprecision(8);
theStream << setw(12) << setiosflags(ios::showpoint) << setiosflags(ios::right) << theAtom.itsX << "\t";
theStream << setw(12) << theAtom.itsY << "\t" ;
theStream << setw(12) << theAtom.itsZ << "\t";
theStream << setprecision(9);
theStream << setw(13) << resetiosflags(ios::showpoint) << setiosflags(ios::left) << theAtom.itsWeight << "\t";
theStream << setw(3) << setprecision(1) << theAtom.itsNuclearSpin << endl;
//#endif
theStream << setprecision(8);
return theStream;
}
// Local definitions
bool InitializeTable ()
{
// Returns TRUE if periodic table successfully read from "thermo.dat"
// Returns FALSE on error
if(isInitialized)
return true;
ElementsFileName = string(BODR_PATH) + "elements.xml";
IsotopesFileName = string(BODR_PATH) + "isotopes.xml";
//ifstream fin(ElementsFileName.c_str());
XMLContextPtr ctx1(new XMLContext);
XMLDocument Elements(ctx1);
XMLContextPtr ctx2(new XMLContext);
XMLDocument Isotopes(ctx2);
XMLNodeList nl, nl_scalar, nl_elements, nl_isotopes;
XMLNodeListIterator it1, it2, it3, it4;
string buffer;
Elements.load_file(ElementsFileName);
Isotopes.load_file(IsotopesFileName);
nl_elements = Elements.firstchild("list")->children("atom");
nl_isotopes = Isotopes.firstchild("cml")->children("isotopeList");
for (it1 = nl_elements.begin(); it1 != nl_elements.end(); it1++)
{
string Name = (*it1)->get_attr("id");
PeriodicTable[Name];
nl_scalar = (*it1)->children("scalar");
for (it2 = nl_scalar.begin(); it2 != nl_scalar.end(); it2++)
{
buffer = (*it2)->get_attr("dictRef");
if (buffer == "bo:atomicNumber")
PeriodicTable[Name].Number = strtol(GetCdata(*it2).c_str(), 0, 10);
else if(buffer == "bo:mass")
PeriodicTable[Name].Weight = strtod(GetCdata(*it2).c_str(), 0);
}
for (it2 = nl_isotopes.begin(); it2 != nl_isotopes.end(); it2++)
{
buffer = (*it2)->get_attr("id");
if (buffer == Name)
{
nl = (*it2)->children("isotope");
for (it3 = nl.begin(); it3 != nl.end(); it3++)
{
Isotope newIsotope;
int n = strtol((*it3)->get_attr("number").c_str(), 0, 10);
newIsotope.Number = n;
newIsotope.Weight = 0;
newIsotope.Abundance = 0;
newIsotope.Spin = 0;
nl_scalar = (*it3)->children("scalar");
if (Control::instance()->Spin())
{
for (it4 = nl_scalar.begin(); it4 != nl_scalar.end(); it4++)
{
buffer = (*it4)->get_attr("dictRef");
if (buffer == "bo:exactMass")
newIsotope.Weight = strtod(GetCdata(*it4).c_str(), 0);
else if (buffer == "bo:relativeAbundance")
newIsotope.Abundance = strtod(GetCdata(*it4).c_str(), 0);
else if (buffer == "bo:spin")
newIsotope.Spin = Spin2Double(GetCdata(*it4));
}
}
else
{
for (it4 = nl_scalar.begin(); it4 != nl_scalar.end(); it4++)
{
buffer = (*it4)->get_attr("dictRef");
if (buffer == "bo:exactMass")
newIsotope.Weight = strtod(GetCdata(*it4).c_str(), 0);
else if (buffer == "bo:relativeAbundance")
newIsotope.Abundance = strtod(GetCdata(*it4).c_str(), 0);
}
}
PeriodicTable[Name].Isotopes[n] = newIsotope;
}
}
}
Hash<int,Isotope>::const_iterator it;
double a = 0;
for (it = PeriodicTable[Name].Isotopes.begin(); it != PeriodicTable[Name].Isotopes.end(); it++)
if ((*it).second.Abundance > a)
{
a = (*it).second.Abundance;
PeriodicTable[Name].itsMainIsotope = (*it).second.Number;
}
}
/*Hash<int,Isotope>::const_iterator it;
for(it = PeriodicTable["B"].Isotopes.begin(); it != PeriodicTable["B"].Isotopes.end(); it++)
cout << "B" << (*it).second.Number << " " << (*it).second.Weight << " " << (*it).second.Abundance
<< " " << (*it).second.Spin << endl;
cout << "Main: " << PeriodicTable["B"].itsMainIsotope << endl;*/
/*
ElementsFileName = "";
// For POSIX: $HOME/.elements.dat
// For WIN32 and others: elements.dat
#ifdef POSIX
ElementsFileName += getenv("HOME");
ElementsFileName += "/.";
#endif
ElementsFileName += "elements.dat";
ifstream fin(ElementsFileName.c_str());
if (!fin)
{
cerr << "File " << ElementsFileName.c_str() << " was not found. Creating new..." << endl;
fin.close();
if (!CreateTable())
return false;
fin.open(ElementsFileName.c_str());
}
if (!fin)
{
string e = "Cannot create ";
e += ElementsFileName;
fin.close();
throw GenericError(e);
return false;
}
fin >> NElements;
for (int k=0; k<NElements; k++)
{
char Name[ATOM_NAME_LENGTH];
if(!fin) return false;
fin >> Name;
Name[0]=toupper(Name[0]); // Element's name starts with capital letter
PeriodicTable[Name];
fin >> PeriodicTable[Name].Weight >> PeriodicTable[Name].Spin;
if (!Control::instance()->Spin())
PeriodicTable[Name].Spin = 0;
#ifdef NO_SPIN
PeriodicTable[Name].Spin = 0;
#endif
PeriodicTable[Name].Number = k+1;
}*/
isInitialized = true;
//fin.close();
return true;
}
/*
bool CreateTable ()
{
ofstream fout(ElementsFileName.c_str());
if(!fout.is_open())
{
string e = "## Unable to write";
e += ElementsFileName;
throw GenericError(e);
}
/* {
#ifndef SILENT
cerr << "\nUnable to write" << ElementsFileName << "\n Exiting";
return false;
#else
return <number_of_error>
#endif
}
fout.precision(10);
fout.width(16);
fout << "55\nH\t1.00794\t\t0.5\n";
fout << "He\t4.002602\t1.5\nLi\t6.941\t\t1.5\nBe\t9.012182\t1.5\nB\t10.811\t\t3.2\nC\t12.0107\t\t0\n";
fout << "N\t14.0067\t\t1\nO\t15.9994\t\t0\nF\t18.9984032\t0.5\nNe\t20.1797\t\t0\n";
fout << "Na\t22.98976928\t1.5\nMg\t24.305\t\t0\nAl\t26.9815386\t2.5\nSi\t28.0855\t\t0\n";
fout << "P\t30.973762\t0.5\nS\t32.065\t\t0\nCl\t35.453\t\t1.5\nAr\t39.948\t\t0\n";
fout << "K\t39.0983\t\t1.5\nCa\t40.078\t\t0\nSc\t44.955912\t3.5\nTi\t47.867\t\t0\nV\t50.9415\t\t3.5\n";
fout << "Cr\t51.9961\t\t0\nMn\t54.938045\t2.5\nFe\t55.845\t\t0\nCo\t58.933195\t3.5\nNi\t58.6934\t\t0\n";
fout << "Cu\t63.546\t\t1.5\nZn\t65.409\t\t0\nGa\t69.723\t\t1.5\nGe\t72.64\t\t0\nAs\t74.9216\t\t1.5\n";
fout << "Se\t78.96\t\t0\nBr\t79.904\t\t1.5\nKr\t83.798\t\t0\n";
fout << "Rb\t85.4678\t\t2.5\nSr\t87.62\t\t0\nY\t88.90585\t0.5\nZr\t91.224\t\t0\nNb\t92.90638\t4.5\n";
fout << "Mo\t95.94\t\t0\nTc\t98.000\t\t6\nRu\t101.07\t\t0\nRh\t102.9055\t0.5\nPd\t106.42\t\t0\nAg\t107.8682\t0.5\n";
fout << "Cd\t112.411\t\t0\nIn\t114.818\t\t4.5\nSn\t118.71\t\t0\nSb\t121.76\t\t2.5\nTe\t127.6\t\t0\n";
fout << "I\t126.90447\t2.5\nXe\t131.293\t\t0\nW\t183.84\t\t0";
fout.close();
return true;
}*/