#include "rotor.hpp"
#include "control.hpp"
#include "mathmodule.hpp"
#include "matrix.h"
#include "symmetry.hpp"
#include <cmath>
#include <iomanip>
#include <fstream>
#include <cstring>
using namespace std;
//extern const double Na, pi, c, h, R;
const double A2m = 1e-10; // Angstrom in meters
// Local declarations
double GenericRotor (double T, double alfa, PES_1D & Barrier);
// Local definitions
double GenericRotor (double T, double alfa, PES_1D & Barrier)
{
double E = Barrier.GetEnergy(alfa);
// energy per mol
return exp(-E/(R*T));
}
double RotationIntegral(double T, PES_1D & Barrier)
{
double I = Integrate(GenericRotor, 0, 2*pi, T, Barrier);
return I;
}
// Definitions of Rotor
Rotor::Rotor(): itsBarrier(0,2*pi)
{
//cerr << "\nRotor constructor...";
itsAlpha = 0;
itsSupport = 0;
itsHead = 0;
itsI0 = 0;
itsIeff = 0;
itsAtomNumber = 0;
itsRotType = HINDERED_POLYNOM;
itsRotationTheory = CLASSIC;
}
Rotor::~Rotor()
{
//cerr << "\nRotor destructor...";
delete [] itsAtomNumber;
}
void Rotor::SetAtoms (int NAtoms)
{
itsAtomNumber = new int[NAtoms];
Structure::SetAtoms(NAtoms);
}
void Rotor::SetAtomNumbers (int * atoms)
{
for (int i=0; i<itsNAtoms; i++)
itsAtomNumber[i] = atoms[i];
}
int Rotor::GetAtomNumber (int n) const
{
return itsAtomNumber[n];
}
void Rotor::SetSupport (Atom * theAtom)
{
itsSupport = theAtom;
}
Atom * Rotor::GetSupport () const
{
return itsSupport;
}
void Rotor::SetHead (int n)
{
itsHead = itsAtom[n];
}
void Rotor::SetHead (Atom * theAtom)
{
itsHead = theAtom;
}
Atom * Rotor::GetHead () const
{
return itsHead;
}
void Rotor::SetSigma (int N)
{
if (N>1)
itsTopType = SYMMETRIC;
Structure::SetSigma(N);
itsBarrier.SetSigma(itsSigma);
}
void Rotor::SetI0 (double I0)
{
itsI0 = I0;
}
double Rotor::GetI0 () const
{
return itsI0;
}
void Rotor::SetIeff (double Ieff)
{
itsIeff = Ieff;
}
double Rotor::GetIeff () const
{
return itsIeff;
}
void Rotor::SetRotType (ROTATION type)
{
itsRotType = type;
}
ROTATION Rotor::GetRotType () const
{
return itsRotType;
}
Rotation_Theory Rotor::GetTheory () const
{
return itsRotationTheory;
}
void Rotor::SetTheory (Rotation_Theory theory)
{
itsRotationTheory = theory;
}
PES & Rotor::GetBarrier ()
{
return itsBarrier;
}
double Rotor::GetInn (Atom * A, Atom * B) const
{
// Returns moment for axis going through A and B
// This function doesn't use itsInertialTensor
double x1, y1, z1, x2, y2, z2, I;
A->GetCoordinates(x1, y1, z1);
B->GetCoordinates(x2, y2, z2);
I = 0;
for (int i=0; i<itsNAtoms; i++)
{
double xi, yi, zi, mi, hi2;
itsAtom[i]->GetCoordinates(xi,yi,zi);
mi = itsAtom[i]->GetWeight()*(1e-3/Na);
hi2 = (xi-x1)*(xi-x1) + (yi-y1)*(yi-y1) + (zi-z1)*(zi-z1) -
((x2-x1)*(xi-x1)+(y2-y1)*(yi-y1)+(z2-z1)*(zi-z1))*((x2-x1)*
(xi-x1)+(y2-y1)*(yi-y1)+(z2-z1)*(zi-z1)) /
((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
I += mi*hi2*A2m*A2m; // Convert A to m !!!
}
return I;
}
double Rotor::A() const
{
return GetInn(itsHead, itsSupport);
}
double Rotor::B() const
{
double b = 0;
for (int i=0; i<itsNAtoms; i++)
{
double xi, yi, zi, mi;
itsAtom[i]->GetCoordinates(xi,yi,zi);
// Convert coordinates!
mi = itsAtom[i]->GetWeight()*(1e-3/Na);
b += mi*xi*zi*A2m*A2m;
}
return b;
}
double Rotor::C() const
{
double c = 0;
for (int i=0; i<itsNAtoms; i++)
{
double xi, yi, zi, mi;
itsAtom[i]->GetCoordinates(xi,yi,zi);
// Convert coordinates!
mi = itsAtom[i]->GetWeight()*(1e-3/Na);
c += mi*yi*zi*A2m*A2m;
}
return c;
}
double Rotor::U() const
{
double u = 0;
for (int i=0; i<itsNAtoms; i++)
{
double xi, yi, zi, mi;
itsAtom[i]->GetCoordinates(xi,yi,zi);
// Convert coordinates!
mi = itsAtom[i]->GetWeight()*(1e-3/Na);
u += mi*xi*A2m;
}
return u;
}
void Rotor::Rotate(double alpha)
{
// Rotate around axis of Rotor
double x1, y1, z1, x2, y2, z2, x, y, z, l;
long double alpha_change = alpha - itsAlpha;
Matrix R(3,3);
ColumnVector V1(3);
ColumnVector V2(3);
itsSupport->GetCoordinates(x1, y1, z1);
itsHead->GetCoordinates(x2, y2, z2);
l = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
x = (x2-x1)/l;
y = (y2-y1)/l;
z = (z2-z1)/l;
// Matrix of rotation
R(1,1) = cos(alpha_change)+(1-cos(alpha_change))*x*x;
R(1,2) = (1-cos(alpha_change))*x*y-sin(alpha_change)*z;
R(1,3) = (1-cos(alpha_change))*x*z+sin(alpha_change)*y;
R(2,1) = (1-cos(alpha_change))*y*x+sin(alpha_change)*z;
R(2,2) = cos(alpha_change)+(1-cos(alpha_change))*y*y;
R(2,3) = (1-cos(alpha_change))*y*z-sin(alpha_change)*x;
R(3,1) = (1-cos(alpha_change))*z*x-sin(alpha_change)*y;
R(3,2) = (1-cos(alpha_change))*z*y+sin(alpha_change)*x;
R(3,3) = cos(alpha_change)+(1-cos(alpha_change))*z*z;
for (int i=0; i<itsNAtoms; i++)
{
double xi, yi, zi;
itsAtom[i]->GetCoordinates(xi,yi,zi);
V1 << xi << yi << zi;
// Transformation of coordinates
V2 = R*V1;
itsAtom[i]->SetCoordinates(V2(1),V2(2),V2(3));
}
itsAlpha = alpha;
}
void Rotor::SetPointGroup (const char * pg)
{
Structure::SetPointGroup(pg);
CalculateSigma();
}
void Rotor::TestSymmetry (const char * LogName, bool exact)
{
char pg[15], pg_brutal[15];
double support[3];
char C1[] = "C1";
#ifdef POSIX
char tmpname[] = "/tmp/openthermo";
#else
char tmpname[] = "tmp";
#endif
SymmetryState s;
ofstream tmp(tmpname);
ofstream fout;
fout.open(LogName, ios::app);
// Preparing input file for Symmetry program
ExportCoordinates(tmp);
tmp.close();
/*// Find number of Head
for (int i=0; i<itsNAtoms; i++)
if (itsAtom[i] == itsHead)
{
head = i;
break;
}
// -1 - exception*/
// Running Symmetry
fout << "\n***** Symmetry (Brute force symmetry analyzer) Author: S. Patchkovskii, 1996,2003 *****" << endl;
fout << "\nStarting Symmetry for rotor...\n" << endl;
fout << "Search thresholds\n\tPrimary: " << DefaultSymmetryPrimary
<< "\n\tFinal: " << DefaultSymmetryFinal
<< "\n\tMaxit: " << DefaultSymmetryMaxit << endl;
fout.close();
itsSupport->GetCoordinates(support[0], support[1], support[2]);
s = symmetry_of_axis(pg, support, tmpname, LogName);
fout.open(LogName, ios::app);
switch (s)
{
case OK:
if (exact != true)
{
fout << "\nRestarting Symmetry to find higher axis...\n" << endl;
fout << "Search thresholds\n\tPrimary: " << Control::instance()->SymmetryPrimary()
<< "\n\tFinal: " << Control::instance()->SymmetryFinal()
<< "\n\tMaxit: " << Control::instance()->SymmetryMaxit() << endl;
fout.close();
s = symmetry_of_axis(pg_brutal, support, tmpname, LogName, -2,
Control::instance()->SymmetryPrimary(),
Control::instance()->SymmetryFinal(),
Control::instance()->SymmetryMaxit());
fout.open(LogName, ios::app);
switch (s)
{
case OK:
SetPointGroup(pg_brutal);
fout << "\nSetting " << pg_brutal << " point group" << endl;
if(strcmp(pg, pg_brutal))
{
fout << "WARNING! " << pg_brutal << " isn't exact point group!" << endl; //\n You can choose
cerr << "# WARNING! Setting " << pg_brutal << ", but it isn't exact point group of rotor" << endl; //\n You can choose
}
goto cleantmp;
case UNKNOWN_GROUP:
fout << "\nSymmetry program can't determine point group, but rotor might have axis of higher order" << endl;
// falling down
default:
SetPointGroup(pg);
fout << "\nSetting " << pg << " point group" << endl;
goto cleantmp;
}
}
else
{
SetPointGroup(pg);
fout << "\nSetting " << pg << " point group" << endl;
goto cleantmp;
}
case UNKNOWN_GROUP:
// exception!
fout << "\nSymmetry program can't determine point group, but rotor might have axis of higher order" << endl;
goto cleantmp;
default:
fout << "Can't recognize point group, C1 is used as default" << endl; //\nUse keyword \'sym\' to specify point group" << endl;
cerr << "# WARNING! Can't recognize point group, C1 is used default" << endl;
SetPointGroup(C1);
goto cleantmp;
}
cleantmp:
fout << "**********" << endl;
tmp.open(tmpname);
tmp.close();
}
void Rotor::LogAtoms (ostream & flog) const
{
flog << "\n\t\t********** Rotor XYZ coordinates **********\n";
Structure::LogAtoms(flog);
flog << "Head atom:\n" << *itsHead;
flog << "Support atom:\n" << *itsSupport;
flog << endl;
}
void Rotor::LogInfo (ostream & flog) const
{
vector<double> coeff, min, max;
flog << "\n***** Rotor Summary *****" << endl;
flog << "Point group: " << itsPointGroup << endl;
switch (itsTopType)
{
case LINEAR:
flog << "Linear top" << endl;
break;
case SYMMETRIC:
flog << "Symmetric top" << endl;
break;
case ASYMMETRIC:
flog << "Asymmetric top" << endl;
break;
default:
flog << "Unknown top\n ##ERROR! Something went wrong!" << endl;
break;
}
flog << "Rotation number: " << itsSigma;
flog << "\nIeff = " << itsIeff*1e7 << " [g*cm^2]" << endl;
switch (itsRotType)
{
case FREE:
flog << "Free rotation" << endl;
break;
case HINDERED_POLYNOM:
break;
case HINDERED_HARMONIC:
case HINDERED_DISCREET:
flog << "Hindered rotation" << endl;
flog << "\nFourier coefficients:\nl\t"<< setw(12) << setiosflags(ios::right)
<< "cos" << "\t" << setw(12) << setiosflags(ios::right) << "sin\n";
coeff = itsBarrier.GetCoeff();
flog << 0 << "\t" << setw(12) << setiosflags(ios::showpoint) << setiosflags(ios::right) << coeff.at(0) << endl;
for (unsigned i=1; i<=(coeff.size()-1)/2; i++)
{
flog << i << "\t" << setw(12) << setiosflags(ios::showpoint) << setiosflags(ios::right) << coeff.at(2*i-1)
<< "\t" << setw(12) << setiosflags(ios::showpoint) << setiosflags(ios::right) << coeff.at(2*i) << endl;
}
flog << endl;
/*flog << "\nApproximation results:\n";
for (int x=0; x<360; x+=15)
{
flog << "\t" << x << "\t";
flog.precision(12);
flog << itsBarrier.GetEnergy(x*pi/180)/1000 << endl;
}*/
if (itsRotationTheory == QUANTUM)
{
flog << "Energy levels\n \t" << setw(12) << setiosflags(ios::right)
<< "[kJ/mol]\t" << setw(12) << setiosflags(ios::right) << "[1/cm]\n";
vector<double> l = itsBarrier.GetLevels();
unsigned int i=0;
do
{
flog << i << "\t" << setw(12) << setiosflags(ios::showpoint)
<< setiosflags(ios::right) << l.at(i)*Na/1000 << "\t" << setw(12)
<< setiosflags(ios::showpoint) << setiosflags(ios::right)
<< l.at(i)/(h*c) << "\n";
i++;
if ((i==21) && (i<l.size()-11))
{
i=l.size()-10;
flog << "...\t...\n";
}
}while (/*(l.at(i) <= itsBarrier.GetHeight()) &&*/ (i<l.size()-1));
flog << endl;
flog << "Overtone frequencies [1/cm]:\n";
for (int i=1; i<20; i++)
flog << i-1 << "->" << i << "\t" << (l.at(i)-l.at(i-1))/(h*c) << endl;
}
if (itsRotType == HINDERED_DISCREET)
{
flog << "\nApproximation results:\n";
for (int x=0; x<360; x+=15)
{
flog << "\t" << x << "\t";
flog.precision(12);
flog << itsBarrier.GetEnergy(x*pi/180)/1000 << endl;
}
/*min = itsBarrier.GetMin();
flog << "Harmonic frequencies in minima:\n";
for (unsigned i=0; i<min.size(); i++)
flog << "Angle: " << min.at(i)/pi*180 << "\tVibr. frequency: "
<< sqrt(itsBarrier.d2Energy(min.at(i))/(Na*itsIeff))/(2*pi*c) << " cm^-1" << endl;*/
min = itsBarrier.GetMin();
max = itsBarrier.GetMax();
flog << "\n\tx\tEnergy[kJ/mol]";
flog << "\nMinima:\n";
for (unsigned int i=0; i<min.size(); i++)
flog << i+1 << "\t" << min.at(i)/pi*180 << "\t" << itsBarrier.GetEnergy(min.at(i))/1e3 << endl;
flog << "\nMaxima:\n";
for (unsigned int i=0; i<min.size(); i++)
flog << max.at(i)/pi*180 << "\t" << itsBarrier.GetEnergy(max.at(i))/1e3 << endl;
}
break;
default:
flog << "Unknown rotation\n ##ERROR! Something went wrong!" << endl;
}
flog << endl;
}
void Rotor::CalculateSigma ()
{
Structure::CalculateSigma();
if (itsSigma == 1)
itsTopType = ASYMMETRIC;
else
itsTopType = SYMMETRIC;
itsBarrier.SetSigma(itsSigma);
}
void Rotor::GetProjectionOfCenter(double &x, double &y, double &z) const
{
double Xa, Ya, Za, Xb, Yb, Zb, Xc, Yc, Zc, alpha;
itsSupport->GetCoordinates(Xa,Ya,Za);
itsHead->GetCoordinates(Xb,Yb,Zb);
GetMassCenter(Xc, Yc, Zc);
alpha = ((Xb-Xa)*Xc+(Yb-Ya)*Yc+(Zb-Za)*Zc)/((Xb-Xa)*(Xb-Xa)+(Yb-Ya)*(Yb-Ya)+(Zb-Za)*(Zb-Za));
x = (Xb-Xa)*alpha;
y = (Yb-Ya)*alpha;
z = (Zb-Za)*alpha;
}
bool Rotor::isBalanced() const
{
double x1,y1,z1,x2,y2,z2;
double d2;
double dmax = Control::instance()->BalancedRotor();
GetMassCenter(x1,y1,z1);
GetProjectionOfCenter(x2,y2,z2);
d2 = (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1);
if (d2 <= dmax*dmax)
return true;
else
return false;
}